tkf · June 3, 2019 01:06 · finmod · May 10, 2017 · tkf · May 28, 2017
diff --git a/examples.ipynb b/examples.ipynb
diff --git a/LyapunovExponentsWithForwardDiff.jl b/LyapunovExponentsWithForwardDiff.jl
 module LyapunovExponentsWithForwardDiff

 using DifferentialEquations
 using ForwardDiff
 using ParameterizedFunctions
 using ProgressMeter
 using RecipesBase

 type LyapunovExponentsResult
    sol_p
    sol_pt
    values
 end

 Base.show(io::IO, res::LyapunovExponentsResult) =
    print(io, "Lyapunov Exponents: ", res.values[:, end])

 type PhaseTangentParam
    phase_dynamics
    Jacobian
    function PhaseTangentParam(phase_dynamics, x0)
        new(phase_dynamics,
            similar(x0, (length(x0), length(x0))))
    end
 end

 function tangent_dynamics(t, u, param, du)
    ForwardDiff.jacobian!(
        param.Jacobian,
        (y, x) -> param.phase_dynamics(t, x, y),
        (@view du[:, 1]),       # phase space derivative goes here
        (@view u[:, 1]),
    )
    A_mul_B!((@view du[:, 2:end]), param.Jacobian, (@view u[:, 2:end]))
 end

 function keepgoing!(sol)
    sol.prob.u0 = sol(sol.prob.tspan[end])
    solve(sol.prob)
 end

 """
 S_n = ((n-1)/n) S_{n-1} + r_n / n
 """
 @inline function lyap_add_R!(n, lyap, R)
    a = (n - 1) / n
    b = 1 - a
    for i = 1:length(lyap)
        lyap[i] = a * lyap[i] + b * log(R[i, i])
    end
 end

 """ A = A * diag(sgn) """
 @inline function A_mul_sign!(A, sgn)
    for i = 1:size(A)[2]
        if !sgn[i]
            A[:, i] *= -1
        end
    end
    A
 end

 """ A = diag(sgn) * A """
 @inline function sign_mul_A!(sgn, A)
    for i = 1:size(A)[1]
        if !sgn[i]
            A[i, :] *= -1
        end
    end
    A
 end


 """ sgn = sign(diag(A))  """
 @inline function sign_diag!(sgn, A)
    for i = 1:size(A)[1]
        sgn[i] = A[i, i] >= 0
    end
    sgn
 end

 """
 `lyapunov_exponents(phase_dynamics, u0, t_chunk, num_tran, num_attr; ...)`

 ### Positional Arguments

 * `phase_dynamics`: Definition of the phase space dynamics in the
  inplace `(t, u, du)` format (as in `ODEProblem`).
 * `u0`: Initial state for the phase space dynamics.
 * `t_chunk`: Length of numerical integration between orthonormalization.
 * `num_tran`: Number of iterations to through away to get rid of the
  transient dynamics.
 * `num_attr`: Number of orthonormalization steps for Lyapunov exponent
  calculation (which is presumably done inside an attractor).

 ### Keyword Arguments

 * `dim_lyap`: Number of Lyapunov exponents to be calculated.
  Default to the full system dimension.
 * `Q0`: The initial guess of the Gram-Schmidt "Lyapunov vectors".
  Default to the identity matrix.
 """
 function lyapunov_exponents(
        phase_dynamics,
        u0,
        t_chunk,
        num_tran,
        num_attr;
        dim_lyap=length(u0),
        Q0=eye(length(u0), dim_lyap))

    # ODEProblem for dynamics in phase space
    pprob = ODEProblem(phase_dynamics, u0, (0.0, t_chunk))

    psol = solve(pprob)
    @showprogress 1 "Transient dynamics..." for _ in 2:num_tran
        psol = keepgoing!(psol)
    end

    # ODEProblem for dynamics in phase & tangent spaces
    pt0 = similar(u0, (length(u0), dim_lyap + 1))
    pt0[:, 1] = psol(psol.prob.tspan[end])  # phase space initial condition
    pt0[:, 2:end] = Q0                      # tangent space ...
    ptparam = PhaseTangentParam(phase_dynamics, u0)
    tprob = ODEProblem(
        ParameterizedFunction(tangent_dynamics, ptparam),
        pt0,
        (0.0, t_chunk),
    )

    # H2 method in Geist, Parlitz, Lauterborn (1990)
    lyap = similar(u0, dim_lyap)
    lehist = similar(u0, (dim_lyap, num_attr))
    signR = Array(Bool, dim_lyap)

    local tsol
    xt = tprob.u0[:, 1]
    Q = tprob.u0[:, 2:end]
    @showprogress 1 "Computing Lyapunov exponents..." for n in 1:num_attr
        tprob.u0[:, 1] = xt
        tprob.u0[:, 2:end] = Q
        tsol = solve(tprob)
        uend = tsol(tsol.prob.tspan[end])
        xt = uend[:, 1]
        P = uend[:, 2:end]

        F = qrfact!(P)
        Q = F[:Q][:, 1:dim_lyap]
        R = F[:R]

        sign_diag!(signR, R)        # signR = diagm(sign(diag(R)))
        A_mul_sign!(Q, signR)       # Q = Q * signR
        sign_mul_A!(signR, R)       # R = signR * R

        lyap_add_R!(n, lyap, R)
        lehist[:, n] = lyap
    end

    lehist /= t_chunk
    LyapunovExponentsResult(
        psol,
        tsol,
        lehist,
    )
 end

 """ Plot `LyapunovExponentsResult` via `RecipesBase`."""
 @recipe function f(res::LyapunovExponentsResult;
                   known=nothing)
    dim_lyap = size(res.values)[1]
    layout --> (dim_lyap, 1)
    xscale --> :log10

    ylims = [[minimum(res.values[i, :]),
              maximum(res.values[i, :])]
             for i = 1:dim_lyap]
    if known != nothing
        for i = 1:dim_lyap
            ylims[i][1] = min(ylims[i][1], known[i])
            ylims[i][2] = max(ylims[i][2], known[i])
        end
    end
    for i = 1:dim_lyap
        ymin, ymax = ylims[i]
        dy = ymax - ymin
        ylims[i] = [ymin - dy * 0.05,
                    ymax + dy * 0.05]
    end

    for i in 1:dim_lyap
        @series begin
            subplot := i
            label --> ""
            ylabel := "LE$i"
            ylim --> ylims[i]
            res.values[i, :]
        end
        if known != nothing
            @series begin
                subplot := i
                linetype := :hline
                label --> ""

                # repeating ylabel/ylim; otherwise they are ignored
                ylabel := "LE$i"
                ylim --> ylims[i]

                [known[i]]
            end
        end
    end
 end

 """ Some example dynamical systems and their known Lyapunov exponents. """
 module Examples

 """
 Lorenz system.

 * https://en.wikipedia.org/wiki/Lorenz_system
 * http://sprott.physics.wisc.edu/chaos/comchaos.htm
 * E. N. Lorenz, J. Atmos. Sci. 20, 130-141 (1963)
 """
 function lorenz(t, u, du)
    du[1] = 10.0(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
 end
 lorenz_les_known = [0.9056, 0, -14.5723]

 """
 Simplest piecewise linear dissipative chaotic flow.

 * http://sprott.physics.wisc.edu/chaos/comchaos.htm
 * S. J. Linz and J. C. Sprott, Phys. Lett. A 259, 240-245 (1999)
 """
 function piecewise_linear(t, u, du)
    du[1] = u[2]
    du[2] = u[3]
    du[3] = -0.6 * u[3] - u[2] - (u[1] > 0 ? u[1] : -u[1]) + 1
 end
 piecewise_linear_les_known = [0.0362, 0, -0.6362]

 end

 end
	module LyapunovExponentsWithForwardDiff

	using DifferentialEquations
	using ForwardDiff
	using ParameterizedFunctions
	using ProgressMeter
	using RecipesBase

	type LyapunovExponentsResult
	sol_p
	sol_pt
	values
	end

	Base.show(io::IO, res::LyapunovExponentsResult) =
	print(io, "Lyapunov Exponents: ", res.values[:, end])

	type PhaseTangentParam
	phase_dynamics
	Jacobian
	function PhaseTangentParam(phase_dynamics, x0)
	new(phase_dynamics,
	similar(x0, (length(x0), length(x0))))
	end
	end

	function tangent_dynamics(t, u, param, du)
	ForwardDiff.jacobian!(
	param.Jacobian,
	(y, x) -> param.phase_dynamics(t, x, y),
	(@view du[:, 1]), # phase space derivative goes here
	(@view u[:, 1]),
	)
	A_mul_B!((@view du[:, 2:end]), param.Jacobian, (@view u[:, 2:end]))
	end

	function keepgoing!(sol)
	sol.prob.u0 = sol(sol.prob.tspan[end])
	solve(sol.prob)
	end

	"""
	S_n = ((n-1)/n) S_{n-1} + r_n / n
	"""
	@inline function lyap_add_R!(n, lyap, R)
	a = (n - 1) / n
	b = 1 - a
	for i = 1:length(lyap)
	lyap[i] = a * lyap[i] + b * log(R[i, i])
	end
	end

	""" A = A * diag(sgn) """
	@inline function A_mul_sign!(A, sgn)
	for i = 1:size(A)[2]
	if !sgn[i]
	A[:, i] *= -1
	end
	end
	A
	end

	""" A = diag(sgn) * A """
	@inline function sign_mul_A!(sgn, A)
	for i = 1:size(A)[1]
	if !sgn[i]
	A[i, :] *= -1
	end
	end
	A
	end


	""" sgn = sign(diag(A)) """
	@inline function sign_diag!(sgn, A)
	for i = 1:size(A)[1]
	sgn[i] = A[i, i] >= 0
	end
	sgn
	end

	"""
	`lyapunov_exponents(phase_dynamics, u0, t_chunk, num_tran, num_attr; ...)`

	### Positional Arguments

	* `phase_dynamics`: Definition of the phase space dynamics in the
	inplace `(t, u, du)` format (as in `ODEProblem`).
	* `u0`: Initial state for the phase space dynamics.
	* `t_chunk`: Length of numerical integration between orthonormalization.
	* `num_tran`: Number of iterations to through away to get rid of the
	transient dynamics.
	* `num_attr`: Number of orthonormalization steps for Lyapunov exponent
	calculation (which is presumably done inside an attractor).

	### Keyword Arguments

	* `dim_lyap`: Number of Lyapunov exponents to be calculated.
	Default to the full system dimension.
	* `Q0`: The initial guess of the Gram-Schmidt "Lyapunov vectors".
	Default to the identity matrix.
	"""
	function lyapunov_exponents(
	phase_dynamics,
	u0,
	t_chunk,
	num_tran,
	num_attr;
	dim_lyap=length(u0),
	Q0=eye(length(u0), dim_lyap))

	# ODEProblem for dynamics in phase space
	pprob = ODEProblem(phase_dynamics, u0, (0.0, t_chunk))

	psol = solve(pprob)
	@showprogress 1 "Transient dynamics..." for _ in 2:num_tran
	psol = keepgoing!(psol)
	end

	# ODEProblem for dynamics in phase & tangent spaces
	pt0 = similar(u0, (length(u0), dim_lyap + 1))
	pt0[:, 1] = psol(psol.prob.tspan[end]) # phase space initial condition
	pt0[:, 2:end] = Q0 # tangent space ...
	ptparam = PhaseTangentParam(phase_dynamics, u0)
	tprob = ODEProblem(
	ParameterizedFunction(tangent_dynamics, ptparam),
	pt0,
	(0.0, t_chunk),
	)

	# H2 method in Geist, Parlitz, Lauterborn (1990)
	lyap = similar(u0, dim_lyap)
	lehist = similar(u0, (dim_lyap, num_attr))
	signR = Array(Bool, dim_lyap)

	local tsol
	xt = tprob.u0[:, 1]
	Q = tprob.u0[:, 2:end]
	@showprogress 1 "Computing Lyapunov exponents..." for n in 1:num_attr
	tprob.u0[:, 1] = xt
	tprob.u0[:, 2:end] = Q
	tsol = solve(tprob)
	uend = tsol(tsol.prob.tspan[end])
	xt = uend[:, 1]
	P = uend[:, 2:end]

	F = qrfact!(P)
	Q = F[:Q][:, 1:dim_lyap]
	R = F[:R]

	sign_diag!(signR, R) # signR = diagm(sign(diag(R)))
	A_mul_sign!(Q, signR) # Q = Q * signR
	sign_mul_A!(signR, R) # R = signR * R

	lyap_add_R!(n, lyap, R)
	lehist[:, n] = lyap
	end

	lehist /= t_chunk
	LyapunovExponentsResult(
	psol,
	tsol,
	lehist,
	)
	end

	""" Plot `LyapunovExponentsResult` via `RecipesBase`."""
	@recipe function f(res::LyapunovExponentsResult;
	known=nothing)
	dim_lyap = size(res.values)[1]
	layout --> (dim_lyap, 1)
	xscale --> :log10

	ylims = [[minimum(res.values[i, :]),
	maximum(res.values[i, :])]
	for i = 1:dim_lyap]
	if known != nothing
	for i = 1:dim_lyap
	ylims[i][1] = min(ylims[i][1], known[i])
	ylims[i][2] = max(ylims[i][2], known[i])
	end
	end
	for i = 1:dim_lyap
	ymin, ymax = ylims[i]
	dy = ymax - ymin
	ylims[i] = [ymin - dy * 0.05,
	ymax + dy * 0.05]
	end

	for i in 1:dim_lyap
	@series begin
	subplot := i
	label --> ""
	ylabel := "LE$i"
	ylim --> ylims[i]
	res.values[i, :]
	end
	if known != nothing
	@series begin
	subplot := i
	linetype := :hline
	label --> ""

	# repeating ylabel/ylim; otherwise they are ignored
	ylabel := "LE$i"
	ylim --> ylims[i]

	[known[i]]
	end
	end
	end
	end

	""" Some example dynamical systems and their known Lyapunov exponents. """
	module Examples

	"""
	Lorenz system.

	* https://en.wikipedia.org/wiki/Lorenz_system
	* http://sprott.physics.wisc.edu/chaos/comchaos.htm
	* E. N. Lorenz, J. Atmos. Sci. 20, 130-141 (1963)
	"""
	function lorenz(t, u, du)
	du[1] = 10.0(u[2]-u[1])
	du[2] = u[1]*(28.0-u[3]) - u[2]
	du[3] = u[1]u[2] - (8/3)u[3]
	end
	lorenz_les_known = [0.9056, 0, -14.5723]

	"""
	Simplest piecewise linear dissipative chaotic flow.

	* http://sprott.physics.wisc.edu/chaos/comchaos.htm
	* S. J. Linz and J. C. Sprott, Phys. Lett. A 259, 240-245 (1999)
	"""
	function piecewise_linear(t, u, du)
	du[1] = u[2]
	du[2] = u[3]
	du[3] = -0.6 * u[3] - u[2] - (u[1] > 0 ? u[1] : -u[1]) + 1
	end
	piecewise_linear_les_known = [0.0362, 0, -0.6362]

	end

	end