text:
@testset "D'Esopo-Pape" begin
    @testset "Generic tests for graphs" begin
        g4 = PathDiGraph(5)
        d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
        d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))

        for g in testdigraphs(g4)
            y = @inferred(desopo_pape_shortest_paths(g, 2, d1))
            z = @inferred(desopo_pape_shortest_paths(g, 2, d2))

            @test y.parents == z.parents == [0, 0, 2, 3, 4]
            @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
        end

        @testset "Graph with cycle" begin
            gx = PathGraph(5)
            add_edge!(gx, 2, 4)
            d = ones(Int, 5, 5)
            d[2, 3] = 100
            for g in testgraphs(gx)
                z = @inferred(desopo_pape_shortest_paths(g, 1, d))
                @test z.dists == [0, 1, 3, 2, 3]
                @test z.parents == [0, 1, 4, 2, 4]
            end
        end

        m = [0 2 2 0 0; 2 0 0 0 3; 2 0 0 1 2;0 0 1 0 1;0 3 2 1 0]
        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 2, 5)
        add_edge!(G, 3, 5)
        add_edge!(G, 3, 4)
        add_edge!(G, 4, 5)

        for g in testgraphs(G)
            y = @inferred(desopo_pape_shortest_paths(g, 1, m))
            @test y.parents == [0, 1, 1, 3, 3]
            @test y.dists == [0, 2, 2, 3, 4]
        end
    end

    @testset "Graph with self loop" begin
        G = SimpleGraph(5)
        add_edge!(G, 2, 2)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 3, 3)
        add_edge!(G, 1, 5)
        add_edge!(G, 2, 4)
        add_edge!(G, 4, 5)
        m = [0 10 2 0 15; 10 9 0 1 0; 2 0 1 0 0; 0 1 0 0 2; 15 0 0 2 0]
        for g in testgraphs(G)
            z = @inferred(desopo_pape_shortest_paths(g, 1 , m))
            y = @inferred(dijkstra_shortest_paths(g, 1, m))
            @test isapprox(z.dists, y.dists)
        end
    end

    @testset "Disconnected graph" begin
        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 4, 5)
        inf = typemax(eltype(G))
        for g in testgraphs(G)
            z = @inferred(desopo_pape_shortest_paths(g, 1))
            @test z.dists == [0, 1, 1, inf, inf]
            @test z.parents == [0, 1, 1, 0, 0]
        end
    end

    @testset "Empty graph" begin
        G = SimpleGraph(3)
        inf = typemax(eltype(G))
        for g in testgraphs(G)
            z = @inferred(desopo_pape_shortest_paths(g, 1))
            @test z.dists == [0, inf, inf]
            @test z.parents == [0, 0, 0]
        end
    end

    @testset "Random Graphs" begin
        @testset "Simple graphs" begin
            for i = 1:5
                nvg = Int(floor(250*rand()))
                neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                seed = Int(floor(100*rand()))
                g = SimpleGraph(nvg, neg; seed = seed)
                z = desopo_pape_shortest_paths(g, 1)
                y = dijkstra_shortest_paths(g, 1)
                @test isapprox(z.dists, y.dists)
            end
        end

        @testset "Simple DiGraphs" begin
            for i = 1:5
                nvg = Int(floor(250*rand()))
                neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                seed = Int(floor(100*rand()))
                g = SimpleDiGraph(nvg, neg; seed = seed)
                z = desopo_pape_shortest_paths(g, 1)
                y = dijkstra_shortest_paths(g, 1)
                @test isapprox(z.dists, y.dists)
            end
        end
    end

    @testset "Different types of graph" begin
        G = CompleteGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = CompleteDiGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = CycleGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = CycleDiGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = StarGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = WheelGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = RoachGraph(9)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        G = CliqueGraph(5, 19)
        z = desopo_pape_shortest_paths(G, 1)
        y = dijkstra_shortest_paths(G, 1)
        @test isapprox(z.dists, y.dists)

        @testset "Small Graphs" begin
            for s in [:bull, :chvatal, :cubical, :desargues,
                      :diamond, :dodecahedral, :frucht, :heawood,
                      :house, :housex, :icosahedral, :krackhardtkite, :moebiuskantor,
                      :octahedral, :pappus, :petersen, :sedgewickmaze, :tutte,
                      :tetrahedral, :truncatedcube, :truncatedtetrahedron, :truncatedtetrahedron_dir]
                G = smallgraph(s)
                z = desopo_pape_shortest_paths(G, 1)
                y = dijkstra_shortest_paths(G, 1)
                @test isapprox(z.dists, y.dists)
            end
        end
    end

end


text:
include("test_schurparlett.jl")
include("test_aaa.jl")
include("test_ratkrylov.jl")

ext: jl


text:
"""
This function adds the ramping limits of generators when there are no CommitmentVariables
"""
function rampconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{ThermalRampLimitDispatch}, system_formulation::Type{S}, time_periods::Int64; args...) where {T <: PowerSystems.ThermalGen, S <: PM.AbstractPowerFormulation}

    devices = [d for d in devices if !isa(d.tech.ramplimits,Nothing)]

    if !isempty(devices)

        p_th = m[:p_th]
        time_index = m[:p_th].axes[2]
        name_index = [d.name for d in devices]

        if :initialpower in keys(args)
            initialpower = args[:initialpower]
        else
            initialpower = Dict([name=>devices[ix].tech.activepower for (ix,name) in enumerate(name_index)])
        end

        (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

        rampdown_th = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)
        rampup_th = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)

        for (ix,name) in enumerate(name_index)
            t1 = time_index[1]
            rampdown_th[name,t1] = @constraint(m,  initialpower[name] - p_th[name,t1] <= devices[ix].tech.ramplimits.down)
            rampup_th[name,t1] = @constraint(m,  p_th[name,t1] - initialpower[name] <= devices[ix].tech.ramplimits.up)
        end

        for t in time_index[2:end], (ix,name) in enumerate(name_index)
            rampdown_th[name,t] = @constraint(m,  p_th[name,t-1] - p_th[name,t] <= devices[ix].tech.ramplimits.down)
            rampup_th[name,t] = @constraint(m,  p_th[name,t] - p_th[name,t-1] <= devices[ix].tech.ramplimits.up)
        end

        JuMP.register_object(m, :rampdown_th, rampdown_th)
        JuMP.register_object(m, :rampup_th, rampup_th)

    else
        @warn "Data doesn't contain generators with ramping limits"

    end

    return m

end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function rampconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64; args...) where {T <: PowerSystems.ThermalGen, D <: AbstractThermalCommitmentForm, S <: AbstractDCPowerModel}

    devices = [d for d in devices if !isa(d.tech.ramplimits,Nothing)]

    if !isempty(devices)

        p_th = m[:p_th]
        start_th = m[:start_th]
        stop_th = m[:stop_th]

        time_index = m[:p_th].axes[2]
        name_index = [d.name for d in devices]

        if :initialpower in keys(args)
            initialpower = args[:initialpower]
        else
            initialpower = Dict([name=>devices[ix].tech.activepower for (ix,name) in enumerate(name_index)])
        end

        (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

        rampdown_th = JuMP.JuMPArray(Array{ConstraintRef}(undef,length(name_index), time_periods), name_index, time_index)
        rampup_th = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)

        for (ix,name) in enumerate(name_index)
            t1 = time_index[1]
            rampdown_th[name,t1] = @constraint(m, initialpower[name] - p_th[name,t1] <= devices[ix].tech.ramplimits.down + devices[ix].tech.activepowerlimits.max * stop_th[name,t1] )
            rampup_th[name,t1] = @constraint(m, p_th[name,t1] - initialpower[name] <= devices[ix].tech.ramplimits.up + devices[ix].tech.activepowerlimits.min * start_th[name,t1] )
        end

        for t in time_index[2:end], (ix,name) in enumerate(name_index)
            rampdown_th[name,t] = @constraint(m, p_th[name,t-1] - p_th[name,t] <= devices[ix].tech.ramplimits.down + devices[ix].tech.activepowerlimits.max * stop_th[name,t] )
            rampup_th[name,t] = @constraint(m, p_th[name,t] - p_th[name,t-1] <= devices[ix].tech.ramplimits.up + devices[ix].tech.activepowerlimits.min * start_th[name,t] )
        end

        JuMP.register_object(m, :rampdown_th, rampdown_th)
        JuMP.register_object(m, :rampup_th, rampup_th)

    else
        @warn "There are no generators with Ramping Limits Data in the System"
    end

    return m

end


text:
abstract type AbstractWakeCombinationModel end

struct LinearFreestreamSuperposition <: AbstractWakeCombinationModel
    # Lissaman 1979
    # new_deficit_sum = old_deficit_sum + wind_speed*deltav

end

struct SumOfSquaresFreestreamSuperposition <: AbstractWakeCombinationModel
    # Katic et al. 1986
    # new_deficit_sum = sqrt(old_deficit_sum**2 + (wind_speed*deltav)**2)

end

struct SumOfSquaresLocalVelocitySuperposition <: AbstractWakeCombinationModel
    # Voutsinas 1990
    # new_deficit_sum = sqrt(old_deficit_sum**2 + (turb_inflow*deltav)**2)

end

struct LinearLocalVelocitySuperposition <: AbstractWakeCombinationModel
    # Niayifar and Porte Agel 2015, 2016
    # new_deficit_sum = old_deficit_sum + turb_inflow*deltav

end

function check_negative_deficits!(new_deficit_sum, wind_speed)
    if new_deficit_sum > wind_speed
        new_deficit_sum = wind_speed - eps(typeof(new_deficit_sum))
    end
end

function wake_combination_model(deltav, wind_speed, turb_inflow, old_deficit_sum,  model::LinearFreestreamSuperposition)
    # Lissaman 1979

    new_deficit_sum = old_deficit_sum + wind_speed*deltav

    check_negative_deficits!(new_deficit_sum, wind_speed)

    return new_deficit_sum

end

function wake_combination_model(deltav, wind_speed, turb_inflow, old_deficit_sum, model::SumOfSquaresFreestreamSuperposition)
    # Katic et al. 1986

    new_deficit_sum = nansafesqrt(old_deficit_sum^2 + (wind_speed*deltav)^2)

    check_negative_deficits!(new_deficit_sum, wind_speed)

    return new_deficit_sum

end

function wake_combination_model(deltav, wind_speed, turb_inflow, old_deficit_sum, model::SumOfSquaresLocalVelocitySuperposition)
    # Voutsinas 1990

    # square root of the sum of the squares
    new_deficit_sum = nansafesqrt(old_deficit_sum^2 + (turb_inflow*deltav)^2)
    # enforce non-negative wind speed
    check_negative_deficits!(new_deficit_sum, turb_inflow)
    return new_deficit_sum

end

function wake_combination_model(deltav, wind_speed, turb_inflow, old_deficit_sum, model::LinearLocalVelocitySuperposition)
    # Niayifar and Porte Agel 2015, 2016

    new_deficit_sum = old_deficit_sum + turb_inflow*deltav

    check_negative_deficits!(new_deficit_sum, turb_inflow)

    return new_deficit_sum

end


text:
using Compat
using LinearOperators

"Abstract type for statistics returned by a solver"
@compat abstract type KrylovStats end;

"Type for statistics returned by non-Lanczos solvers"
mutable struct SimpleStats <: KrylovStats
  solved :: Bool
  inconsistent :: Bool
  residuals :: Array{Float64,1}
  Aresiduals :: Array{Float64,1}
  status :: String
end

export cgTN
"""A standard implementation of the Conjugate Gradient method.
The only non-standard point about it is that it does not check
that the operator is definite.
It is possible to check that the system is inconsistent by
monitoring ‖p‖, which would cost an extra norm computation per
iteration.

Dominique Orban, <dominique.orban@gerad.ca>
Salt Lake City, UT, March 2015.

Try to adapt a stopping criterion based on regulα for the ARCq.

First reformulate the TR case using the characteristic function to confirm it is
equivalent to the usual implementation. A must be symmetric to define a quadratic objective
q(x) = 0.5*x'*A*x - b'*x

  JPD february 09 2017, Montréal"""
cgTN(A::Array{TA, 2}, b::Array{Tb, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, verbose::Bool = false) where {TA <: Number, Tb <: Number} =
  cgTN(LinearOperator(A), b, atol=atol, rtol=rtol, itmax=itmax, verbose=verbose);

cgTN(A::SparseMatrixCSC{TA, IA}, b::Array{Tb, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, verbose::Bool = false) where {TA <: Number, Tb <: Number, IA <: Integer} =
  cgTN(LinearOperator(A), b, atol=atol, rtol=rtol, itmax=itmax, verbose=verbose);

"""The conjugate gradient method to solve the symmetric linear system Ax=b.

The method does _not_ abort if A is not definite.
"""
function cgTN(A::LinearOperator, b::Array{T, 1}; atol::Float64 = 1e-08, rtol::Float64 = 1e-06, itmax::Int = 0, verbose::Bool = false) where T <: Real
    n = size(b, 1);
    (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size");
    #isequal(triu(A)',tril(A)) || error("Must supply Hermitian matrix")

    verbose && @printf("CG: system of %d equations in %d variables\n", n, n);

    # Initial state.
    x = zeros(n)
    # x = zero(n)
    x̂ = copy(x)

    γ = dot(b, b);
    γ == 0 && return x;
    r = copy(b);
    p = copy(r);

    σ = 0.0

    iter = 0;
    itmax == 0 && (itmax = 2 * n);

    rNorm = sqrt(γ);
    rNorms = [rNorm;];
    ε = atol + rtol * rNorm;
    verbose && @printf("%5d  %8.1e ", iter, rNorm)

    solved = rNorm <= ε;
    tired = iter >= itmax;
    neg_curv = false;
    status = "unknown";

    #q = s ->  0.5 * dot(s, copy(A * s)) - dot(b, s)

    #m = s ->  q(s) + norm(s)^3/(3*regulα)
    #hO = α -> m(x+α*p)
    Ap = copy(A * p);  # Bug in LinearOperators? A side effect spoils the computation without the copy.
    # pAp = BLAS.dot(n, p, 1, Ap, 1);
    pAp = p ⋅ Ap
    if pAp<=0
        neg_curv = true
        status = "gradient negative curvature"
        stats = SimpleStats(solved, false, rNorms, [], status);
        return (p, stats);
    else
        α = γ / pAp;
        # BLAS.axpy!(n,  α,  p, 1, x, 1);  # Faster than x = x + σ * p;
        @. x = x + σ * p;
        # BLAS.axpy!(n, -α, Ap, 1, r, 1);  # Faster than r = r - α * Ap;
        @. r = r - α * Ap;
        γ_next = BLAS.dot(n, r, 1, r, 1);
        rNorm = sqrt(γ);
        push!(rNorms, rNorm);

        solved = (rNorm <= ε) ;
        tired = iter >= itmax;

        if !solved
            β = γ_next / γ;
            γ = γ_next;
            BLAS.scal!(n, β, p, 1)
            # BLAS.axpy!(n, 1.0, r, 1, p, 1);  # Faster than p = r + β * p;
            @. p = r + β * n
        end
        iter = iter + 1;
    end

    while ! (solved || tired || neg_curv)
        Ap = copy(A * p);  # Bug in LinearOperators? A side effect spoils the computation without the copy.
        pAp = BLAS.dot(n, p, 1, Ap, 1);

        α = γ / pAp;

        verbose && @printf("%8.1e  %7.1e  %7.1e\n", pAp, α, σ);

        # Move along p from x to the min if either
        # the next step leads farther than the min of the regularized model or
        # we have nonpositive curvature.
        if (pAp <= 0.0)
            α = 0.0
            neg_curv = true
            verbose &&println("negative curvature encountered ",iter," CG iterations.")
        else
            verbose && @printf("    %8.1e  %7.1e  %7.1e\n", pAp, α, σ);
            BLAS.axpy!(n,  α,  p, 1, x, 1);  # Faster than x = x + σ * p;
            # @. x = x + σ * p;
            BLAS.axpy!(n, -α, Ap, 1, r, 1);  # Faster than r = r - α * Ap;
            # @. r = r - α * Ap;
            γ_next = BLAS.dot(n, r, 1, r, 1);
            rNorm = sqrt(γ);
            push!(rNorms, rNorm);

            solved = (rNorm <= ε) ;
            tired = iter >= itmax;
        end

        if !(solved || neg_curv)
            β = γ_next / γ;
            γ = γ_next;
            BLAS.scal!(n, β, p, 1)
            BLAS.axpy!(n, 1.0, r, 1, p, 1);  # Faster than p = r + β * p;
            # @. p = r + β * p;
        end
        iter = iter + 1;
        verbose && @printf("%5d  %8.1e ", iter, rNorm);
    end
    verbose && @printf("\n");

  status = neg_curv ? "negative curvature avoided" : (tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, false, rNorms, [], status);
  return (x, stats);
end


text:
immutable Octonion{T<:Real} <: Number
  s::T
  v1::T
  v2::T
  v3::T
  v4::T
  v5::T
  v6::T
  v7::T
  norm::Bool
end

Octonion( s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real, n::Bool = false) =
  Octonion( promote( s, v1, v2, v3, v4, v5, v6, v7 )..., n )
Octonion( x::Real ) = Octonion( x, zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), abs(x) == one(x) )
Octonion( z::Complex ) = Octonion( z.re, z.im, zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), abs(z) == one(z.re) )
Octonion( q::Quaternion ) = Octonion( q.s, q.v1, q.v2, q.v3, zero(q.s), zero(q.s), zero(q.s), zero(q.s), q.norm )
Octonion( s::Real, a::Vector ) = Octonion( s, a[1], a[2], a[3], a[4], a[5], a[6], a[7] )
Octonion( a::Vector ) = Octonion( 0, a[1], a[2], a[3], a[4], a[5], a[6], a[7] )

convert{T}(::Type{Octonion{T}}, x::Real) =
  Octonion(convert(T,x), convert(T,0), convert(T,0), convert(T,0), convert(T,0), convert(T,0), convert(T,0), convert(T,0))
convert{T}(::Type{Octonion{T}}, z::Complex) =
  Octonion(convert(T,real(z)), convert(T,imag(z)), convert(T,0), convert(T,0), convert(T,0), convert(T,0), convert(T,0), convert(T,0))
convert{T}(::Type{Octonion{T}}, q::Quaternion) =
  Octonion(convert(T,real(z)), convert(T,q.v1), convert(T,q.v2), convert(T,q.v3), convert(T,0), convert(T,0), convert(T,0), convert(T,0))
convert{T<:Real}(::Type{Octonion{T}}, o::Octonion{T}) = o
convert{T}(::Type{Octonion{T}}, o::Octonion) =
  Octonion(convert(T,o.s), convert(T,o.v1), convert(T,o.v2), convert(T,o.v3), convert(T,o.v4), convert(T,o.v5), convert(T,o.v6), convert(T,o.v7), o.norm)

promote_rule{T<:Real}(::Type{Octonion{T}}, ::Type{T}) = Octonion{T}
promote_rule{T<:Real}(::Type{Octonion}, ::Type{T}) = Octonion
promote_rule{T<:Real,S<:Real}(::Type{Octonion{T}}, ::Type{S}) = Octonion{promote_type(T,S)}
promote_rule{T<:Real,S<:Real}(::Type{Complex{T}}, ::Type{Octonion{S}}) = Octonion{promote_type(T,S)}
promote_rule{T<:Real,S<:Real}(::Type{Quaternion{T}}, ::Type{Octonion{S}}) = Octonion{promote_type(T,S)}
promote_rule{T<:Real,S<:Real}(::Type{Octonion{T}}, ::Type{Octonion{S}}) = Octonion{promote_type(T,S)}

octo( p, v1, v2, v3, v4, v5, v6, v7, n = false ) = Octonion( p, v1, v2, v3, v4, v5, v6, v7, n )
octo( x ) = Octonion( x )
octo( s, a ) = Octonion( s, a )

function show(io::IO, o::Octonion)
  pm(x) = x < 0 ? " - $(-x)" : " + $x"
  print(io, o.s, pm(o.v1), "im", pm(o.v2), "jm", pm(o.v3), "km", pm(o.v4), "ilm", pm(o.v5), "jlm", pm(o.v6), "klm", pm(o.v7), "lm")
end

real( o::Octonion ) = o.s
imag( o::Octonion ) = [ o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7 ]

(/)( o::Octonion, x::Real ) = Octonion( o.s/x, o.v1/x, o.v2/x, o.v3/x, o.v4/x, o.v5/x, o.v6/x, o.v7/x )

conj( o::Octonion ) = Octonion( o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm )
abs( o::Octonion ) = sqrt(o.s*o.s + o.v1*o.v1 + o.v2*o.v2 + o.v3*o.v3 + o.v4*o.v4 + o.v5*o.v5 + o.v6*o.v6 + o.v7*o.v7 )
abs2( o::Octonion ) = o.s*o.s + o.v1*o.v1 + o.v2*o.v2  + o.v3*o.v3 + o.v4*o.v4 + o.v5*o.v5 + o.v6*o.v6 + o.v7*o.v7
inv( o::Octonion ) = o.norm ? conj(o) : conj(o)/abs2(o)

function normalize( o::Octonion )
  if ( o.norm )
    return o
  end
  o = o / abs( o )
  Octonion( o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true )
end

function normalizea( o::Octonion )
  if ( o.norm )
    return ( o, one( o.s ) )
  end
  a = abs( o )
  o = o / a
  ( Octonion( o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true ), a )
end

(-)( o::Octonion) = Octonion( -o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm )

(+)( o::Octonion, w::Octonion ) = Octonion( o.s + w.s,
                                            o.v1 + w.v1,
                                            o.v2 + w.v2,
                                            o.v3 + w.v3,
                                            o.v4 + w.v4,
                                            o.v5 + w.v5,
                                            o.v6 + w.v6,
                                            o.v7 + w.v7 )

(-)( o::Octonion, w::Octonion ) = Octonion( o.s - w.s,
                                            o.v1 - w.v1,
                                            o.v2 - w.v2,
                                            o.v3 - w.v3,
                                            o.v4 - w.v4,
                                            o.v5 - w.v5,
                                            o.v6 - w.v6,
                                            o.v7 - w.v7 )

macro ocmul( i, j )
  si = @compat Symbol( "v$i" )
  sj = @compat Symbol( "v$j" )
  esc( :( o.$si * w.$sj - o.$sj * w.$si ) )
end

(*)( o::Octonion, w::Octonion ) =
  Octonion(
            o.s*w.s - o.v1*w.v1 - o.v2*w.v2 - o.v3*w.v3 - o.v4*w.v4 - o.v5*w.v5 - o.v6*w.v6 - o.v7*w.v7,
            o.s*w.v1 + o.v1*w.s + @ocmul( 2, 3 ) + @ocmul( 6, 5 ) + @ocmul( 7, 4 ),
            o.s*w.v2 + o.v2*w.s + @ocmul( 3, 1 ) + @ocmul( 4, 6 ) + @ocmul( 7, 5 ),
            o.s*w.v3 + o.v3*w.s + @ocmul( 1, 2 ) + @ocmul( 5, 4 ) + @ocmul( 7, 6 ),
            o.s*w.v4 + o.v4*w.s + @ocmul( 1, 7 ) + @ocmul( 3, 5 ) + @ocmul( 6, 2 ),
            o.s*w.v5 + o.v5*w.s + @ocmul( 1, 6 ) + @ocmul( 2, 7 ) + @ocmul( 4, 3 ),
            o.s*w.v6 + o.v6*w.s + @ocmul( 2, 4 ) + @ocmul( 3, 7 ) + @ocmul( 5, 1 ),
            o.s*w.v7 + o.v7*w.s + @ocmul( 4, 1 ) + @ocmul( 5, 2 ) + @ocmul( 6, 3 )
          )

(/)( o::Octonion, w::Octonion ) = o * inv(w)

function exp( o::Octonion )
  s = o.s
  se = exp( s )
  scale = se
  th = abs( Octonion( imag( o)  ) )
  if th > 0
    scale *= sin( th ) / th
  end
  Octonion( se * cos( th ),
            scale * o.v1,
            scale * o.v2,
            scale * o.v3,
            scale * o.v4,
            scale * o.v5,
            scale * o.v6,
            scale * o.v7,
            abs( s ) == 0 )
end

function log( o::Octonion )
  o, a = normalizea( o )
  s = o.s
  M = abs( Octonion( imag( o ) ) )
  th = atan2( M, s )
  if M > 0
    M = th / M
    return Octonion( log( a ),
                     o.v1 * M,
                     o.v2 * M,
                     o.v3 * M,
                     o.v4 * M,
                     o.v5 * M,
                     o.v6 * M,
                     o.v7 * M )
  else
    return Octonion( log( a ), th , 0.0, 0.0 )
  end
end

(^)(o::Octonion, w::Octonion) = exp( w * log( o ) )

function sqrt( o::Octonion )
  exp( 0.5 * log( o ) )
end

octorand() = octo( randn(), randn(), randn(), randn(), randn(), randn(), randn(), randn() )


text:
using WarpedImageSeries, RegisterDeformation, CoordinateTransformations

using Test

all_nan_eq(a1, a2) = all(map((x,y)->isequal(x,y), a1, a2))

#2D
n = 100
a = rand(100,100,n);
tfms = fill(Translation(0.1,0.2), n);
ws = warpedseries(a, tfms);
@test size(a) == size(ws)
wrpd = warp(a[:,:,1], tfms[1]);
@test all_nan_eq(ws[:,:,1], wrpd[Base.front(axes(a))...])
wrpd = warp(a[:,:,2], tfms[2]);
@test all_nan_eq(ws[:,:,2], wrpd[Base.front(axes(a))...])

tfms2 = fill(Translation(1.1,1.2), n);
ws2 = warpedseries(a, tfms2);
@test size(ws2,1) == size(a,1) - 1
@test size(ws2,2) == size(a,2) - 1

sls = ws2[:,:,2:4];
for (n,i) in enumerate(2:4)
    @test all_nan_eq(ws2[:,:,i], sls[:,:,n])
end
#@time hm = ws2[:,:,1];
#@time hm = ws2[:,:,1:6];

#3D
n = 10
a = rand(10,10,10,n);
tfms = fill(Translation(0.1,0.2,0.15), n);
ws = warpedseries(a, tfms);
@test size(a) == size(ws)
wrpd = warp(a[:,:,:,1], tfms[1]);
@test all_nan_eq(ws[:,:,:,1], wrpd[Base.front(axes(a))...])
wrpd = warp(a[:,:,:,2], tfms[2]);
@test all_nan_eq(ws[:,:,:,2], wrpd[Base.front(axes(a))...])

tfms2 = fill(Translation(1.1,1.2,1.15), n);
ws2 = warpedseries(a, tfms2);
@test size(ws2,1) == size(a,1) - 1
@test size(ws2,2) == size(a,2) - 1
@test size(ws2,3) == size(a,3) - 1

sls = ws2[:,:,:,2:4];
for (n,i) in enumerate(2:4)
    @test all_nan_eq(ws2[:,:,:,i], sls[:,:,:,n])
end


text:
abstract Type2FuzzySet <: AbstractFuzzySet

#####################
abstract IntervalType2FuzzySet <: Type2FuzzySet
import Base: getindex

function getindex(mf::IntervalType2FuzzySet, range)
	out = hcat([[mf[value]::Tuple{Float64,Float64}...] for value in range]...)
	(out[1,:],out[2,:])
end

getindex(mf::IntervalType2FuzzySet, values...) = getindex(mf,[values...])

export getindex

#####################
# Interval 2 Membership Functions
immutable GaussType2 <: IntervalType2FuzzySet
	upper::GaussMF
	lower::GaussMF
	GaussType2(muU,muL,sigmaU,sigmaL) = new(GaussMF(muU,sigmaU),GaussMF(muL,sigmaL))
	GaussType2(upper,lower) = new(upper,lower)

end
export GaussType2

import Base: getindex
function getindex(mf::GaussType2, value::Real)
	m1 = (mf.lower.mu <= mf.upper.mu) ? mf.lower : mf.upper
	m2 = (mf.lower.mu > mf.upper.mu) ? mf.lower : mf.upper
	if m1.mu <= value <= (m1.mu+m2.mu)/2
		mu1 = 1
		mu2 = m2[value]
	elseif (m1.mu+m2.mu)/2 < value <= m2.mu
		mu1 = 1
		mu2 = m1[value]
	else
		mu1 = m1[value]
		mu2 = m2[value]
	end
	return (max(mu1,mu2),min(mu1,mu2))
end
export getindex

immutable TriangleType2 <: IntervalType2FuzzySet
	lower::TriangleMF
	upper::TriangleMF
	# TriangleType2(muU,muL,sigmaU,sigmaL) = new(GaussMF(muU,sigmaU),GaussMF(muL,sigmaL),muU,muL)
end
export TriangleType2

import Base: getindex
function getindex(mf::TriangleType2, value::Real)
	m1 = (mf.lower.b <= mf.upper.b) ? mf.lower : mf.upper
	m2 = (mf.lower.b > mf.upper.b) ? mf.lower : mf.upper
	if m1.b <= value <= (m1.b+m2.b)/2
		mu1 = 1
		mu2 = m2[value]
	elseif (m1.b+m2.b)/2 < value <= m2.b
		mu1 = 1
		mu2 = m1[value]
	else
		mu1 = m1[value]
		mu2 = m2[value]
	end
	return (max(mu1,mu2),min(mu1,mu2))
end
export getindex

#####################


text:
function effective_excitation_hamiltonian(trivial::Bool, ham::MPOHamiltonian, Bs, p::Float64, mpsleft::InfiniteMPS, paramsleft, mpsright::InfiniteMPS, paramsright; renorm=true)
    #does not "know" that B is left gauged, so it's possible to speed this up further
    toret = zero.(Bs)
    len = length(mpsleft);

    ids = collect(Iterators.filter(x->isid(ham,x),2:ham.odim-1));

    #validity checks
    @assert length(Bs) == len
    @assert length(mpsright) == len
    @assert ham.period == len
    @assert eltype(mpsleft)==eltype(mpsright)

    #build lBs(c)
    lBs = [ TensorMap(zeros,eltype(mpsleft),
                    space(paramsleft.lw[1,k],1)*space(paramsleft.lw[1,k],2),
                    space(Bs[1],3)'*space(paramsright.lw[1,k],3)') for k in 1:ham.odim]
    lBsc = []

    for pos = 1:len
        lBs = exci_transfer_left(lBs,ham,pos,mpsright.AR[pos],mpsleft.AL[pos])*exp(conj(1im*p))
        lBs += exci_transfer_left(leftenv(paramsleft,pos,mpsleft),ham,pos,Bs[pos],mpsleft.AL[pos])*exp(conj(1im*p))

        if trivial
            for i in ids
                @tensor lBs[i][-1,-2,-3,-4] -= lBs[i][1,-2,-3,2]*r_RL(mpsleft,pos)[2,1]*l_RL(mpsleft,pos)[-1,-4]
            end
        end

        push!(lBsc,lBs)
    end

    #build rBs(c)
    rBs = [ TensorMap(zeros,eltype(mpsleft),
                    space(paramsleft.rw[end,k],1)*space(Bs[1],3),
                    space(paramsright.rw[end,k],2)'*space(paramsright.rw[end,k],3)') for k in 1:ham.odim]
    rBsc = []

    for pos=len:-1:1
        rBs = exci_transfer_right(rBs,ham,pos,mpsleft.AL[pos],mpsright.AR[pos])*exp(1im*p)
        rBs += exci_transfer_right(rightenv(paramsright,pos,mpsright),ham,pos,Bs[pos],mpsright.AR[pos])*exp(1im*p)

        if trivial
            for i in ids
                @tensor rBs[i][-1,-2,-3,-4] -= rBs[i][1,-2,-3,2]*l_LR(mpsleft,pos)[2,1]*r_LR(mpsleft,pos)[-1,-4]
            end
        end

        push!(rBsc,rBs)
    end
    rBsc=reverse(rBsc)

    # B in same unit cell as B'
    # this is the only point where we have to take renorm into account (a constant shift in the hamiltonian will only affect the terms where both Bs are at the same position)
    for i = 1:len
        for (j,k) in keys(ham,i)
            @tensor toret[i][-1,-2,-3,-4] +=    leftenv(paramsleft,i,mpsleft)[j][-1,1,2]*
                                                Bs[i][2,3,-3,4]*
                                                ham[i,j,k][1,-2,5,3]*
                                                rightenv(paramsright,i,mpsright)[k][4,5,-4]

            if (renorm) # <B|H|B>-<H>
                en = @tensor    conj(mpsleft.AC[i][11,12,13])*
                                leftenv(paramsleft,i,mpsleft)[j][11,1,2]*
                                mpsleft.AC[i][2,3,4]*
                                ham[i,j,k][1,12,5,3]*
                                rightenv(paramsleft,i,mpsleft)[k][4,5,13]
                toret[i] -= Bs[i]*en
            end
            if i>1
                @tensor toret[i][-1,-2,-3,-4] +=    lBsc[i-1][j][-1,1,-3,2]*
                                                    mpsright.AR[i][2,3,4]*
                                                    ham[i,j,k][1,-2,5,3]*
                                                    rightenv(paramsright,i,mpsright)[k][4,5,-4]
            end
            if i<len
                @tensor toret[i][-1,-2,-3,-4] +=    leftenv(paramsleft,i,mpsleft)[j][-1,1,2]*
                                                    mpsleft.AL[i][2,3,4]*
                                                    ham[i,j,k][1,-2,5,3]*
                                                    rBsc[i+1][k][4,-3,5,-4]
            end
        end

    end

    #B left to B'; outside the unit cell
    lBsE = left_excitation_transfer_system(lBs,ham,mpsleft,mpsright,trivial,p)

    for i=1:len
        for (j,k) in keys(ham,i)
            @tensor toret[i][-1,-2,-3,-4] +=    lBsE[j][-1,1,-3,2]*
                                                mpsright.AR[i][2,3,4]*
                                                ham[i,j,k][1,-2,5,3]*
                                                rightenv(paramsright,i,mpsright)[k][4,5,-4]
        end

        lBsE = exci_transfer_left(lBsE,ham,i,mpsright.AR[i],mpsleft.AL[i])*exp(conj(1im*p))

        if trivial
            for k in ids
                @tensor lBsE[k][-1,-2,-3,-4] -= lBsE[k][1,-2,-3,2]*r_RL(mpsleft,i)[2,1]*l_RL(mpsleft,i)[-1,-4]
            end
        end
    end

    #B right to B'; outside the unit cell
    rBsE = right_excitation_transfer_system(rBs,ham,mpsleft,mpsright,trivial,p)

    for i=len:-1:1
        for (j,k) in keys(ham,i)
            @tensor toret[i][-1,-2,-3,-4] +=    leftenv(paramsleft,i,mpsleft)[j][-1,1,2]*
                                                mpsleft.AL[i][2,3,4]*
                                                ham[i,j,k][1,-2,5,3]*
                                                rBsE[k][4,-3,5,-4]
        end

        rBsE=exci_transfer_right(rBsE,ham,i,mpsleft.AL[i],mpsright.AR[i])*exp(1im*p)

        if trivial
            for k in ids
                @tensor rBsE[k][-1,-2,-3,-4]-=rBsE[k][1,-2,-3,2]*l_LR(mpsleft,i)[2,1]*r_LR(mpsleft,i)[-1,-4]
            end
        end
    end

    return toret
end


text:
using ConfParser
using Test

conf = ConfParse(joinpath(@__DIR__, "confs", "config.keyonly"))
parse_conf!(conf)

@test retrieve(conf, "key") == "e613ef71d63b84b721bdd345a5708ce5738028"
@test haskey(conf, "key") == true
@test haskey(conf, "nothere") == false
@test erase!(conf, "key")
@test commit!(conf,"key","e613ef71d63b84b721bdd345a5708ce5738028")
save!(conf, joinpath(@__DIR__, "confs", "out.conf"))

conf = ConfParse(joinpath(@__DIR__, "confs", "out.conf"))
parse_conf!(conf)
@test retrieve(conf, "key")     == "e613ef71d63b84b721bdd345a5708ce5738028"


text:
using Spring

function go()
    sp=SpringParam(ls=0.1,ms=1,ks=2.5,ns=400)
    ip=InitParam(λ=2.0,shift=0.1,pos=sp.ls*sp.ns/2)
    ap=AnimParam(δt=0.1,nδt=5000,nδtperframe=100)

    V=Array{Float64,2}
    @time animate_spring2D(sp,ip,ap,V)
end
go()


text:
using LazySets
using LazySets.Arrays
using Plots, LaTeXStrings, Plots.PlotMeasures
using LinearAlgebra

(@isdefined TARGET_FOLDER) ? nothing : TARGET_FOLDER = ""
include("Clamped_Model.jl")

NSTEPS = 12000
nMasses = 1000

# computed parameters
finalTime = NSTEPS * dtn
tdom = range(0, finalTime, length = NSTEPS + 1)

PosTwoThirds = 700
PosMiddle    = round(Int, nMasses / 2)
VelTwoThirds = PosTwoThirds + nMasses
VelMiddle    = PosMiddle + nMasses

U₀ = zeros(nMasses)
V₀ = zeros(nMasses)
p = InitialValueProblem(s1000, (U₀, V₀))

println("Solving with Bathe method")
@time solBathe = solve(p, Bathe(dtn), NSTEPS=NSTEPS)

println("Solving with Newmark method")
@time solNewmark = solve(p, Trapezoidal(dtn), NSTEPS=NSTEPS)

uBathe_PosTwoThirds = displacements(solBathe, PosTwoThirds)
uBathe_VelTwoThirds = velocities(solBathe, PosTwoThirds)

uNewmark_PosTwoThirds = displacements(solNewmark, PosTwoThirds)
uNewmark_VelTwoThirds = velocities(solNewmark, PosTwoThirds)

# =====================
# Reachability results
# =====================

p = InitialValueProblem(s1000, Singleton(zeros(2*nMasses)))
println("Solving with ORBIT method")
@time sol_nobloating = solve(p, NSTEPS=NSTEPS, alg=ORBIT(δ=dtn))

println("Solving with LGG09 method")
@time Mout = _sol_reachability_clamped(s1000);
sol_PosTwoThirds = [ReachSet(Interval(-Mout[2, k], Mout[1, k]), (0.0 .. dtn) + (k-1) * dtn) for k in 1:size(Mout, 2)] |> Flowpipe
sol_VelTwoThirds = [ReachSet(Interval(-Mout[4, k], Mout[3, k]), (0.0 .. dtn) + (k-1) * dtn) for k in 1:size(Mout, 2)] |> Flowpipe

# =========================
# graficos desplazamientos
# =========================

# time window where the analytic solution is computed
tdomasol = range(0.0, finalTime, length=10_000)

asol_PosTwoThirds, asol_VelTwoThirds = _sol_analitica_clamped(PosTwoThirds, tdomasol, 500)

function plot_displacements()

    fig = plot(xlab=L"\textrm{Time}", ylab=L"\textrm{Displacement of node } %$PosTwoThirds",
               legend=:bottomright,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(20, "Times"),
               yguidefont=font(20, "Times"),
               xtick=([0, 0.0025, 0.005, 0.0075, 0.01],[L"0.000", L"0.0025", L"0.0050", L"0.0075", L"0.01"]),
               ytick=([0, 0.02, 0.04, 0.06, 0.08, 0.1], [L"0.00", L"0.02", L"0.04", L"0.06", L"0.08", L"0.1"]),
               xlims=(0.0, 0.012), ylims=(-0.002, 0.1),
               bottom_margin=6mm, left_margin=20mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    # flowpipe
    plot!(fig, sol_PosTwoThirds, vars=(0, 1), c=:lightblue, lc=:lightblue, lw=4.0, alpha=1., lab="") # lab=L"\textrm{Flowpipe}")

    # plot!(fig, tdom, uNewmark, seriestype=:path, marker=:none, color=:red, lab=L"\textrm{Newmark}", lw=1.0)
    # plot!(fig, tdom, uBathe, seriestype=:path, marker=:none, color=:green, lab=L"\textrm{Bathe}", lw=1.0)
    # plot!(fig, sol_nobloating, vars=(0, PosTwoThirds), seriestype=:path, marker=:none, lw = 2.0, c=:blue, lab=L"\textrm{Analytic (ODE)}")

    tZoomDispLeft  = 0.0055
    tZoomDispRight = 0.0058

    yZoomDispBott  = 0.08
    yZoomDispTop = 0.098

    # == Analytic ODE ==
    plot!(fig, sol_nobloating, vars=(0, PosTwoThirds), seriestype=:path, marker=:none, lw = 1.5, alpha=1, c=:blue, lab="")
    #plot!(fig, sol_nobloating[ai:paso:bi], vars=(0, VelTwoThirds), linetype=:scatter, marker=:square, markersize=5,  c=:blue, lab=L"\textrm{Analytic (ODE)}")

    plot!(fig, tdomasol, asol_PosTwoThirds, c=:magenta, lw=1., lab="") # , lab=L"\textrm{Analytic (PDE)}")

    PBotLeft  = [tZoomDispLeft,yZoomDispBott]
    PBotRight = [tZoomDispRight,yZoomDispBott]
    PTopLeft  = [tZoomDispLeft,yZoomDispTop]
    PTopRight = [tZoomDispRight,yZoomDispTop]

    L1 = LineSegment( PBotLeft, PBotRight )
    L2 = LineSegment( PBotRight, PTopRight )
    L3 = LineSegment( PTopRight, PTopLeft )
    L4 = LineSegment( PTopLeft, PBotLeft )

    plot!(fig, L1, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L2, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L3, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L4, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
end

fig = plot_displacements()
savefig(fig, joinpath(TARGET_FOLDER, "u_node700.pdf"))

# ================================
# graficos desplazamientos (ZOOM)
# ================================

function plot_displacements_zoom()

    tZoomDispLeft  = 0.00559
    tZoomDispRight = 0.00572

    idxtleft = findfirst(t -> t > tZoomDispLeft, tdom)
    idxtright = findfirst(t -> t > tZoomDispRight, tdom)

    yZoomDispBott  = 0.0925
    yZoomDispTop = 0.0937

    fig = plot(xlab=L"\textrm{Time } (\times 10^{-3})", ylab=L"\textrm{Displacement of node } %$PosTwoThirds",
               legend=:bottomright,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(20, "Times"),
               yguidefont=font(20, "Times"),
               xtick=([0.00560, 0.00562, 0.00564, 0.00566, 0.00568, 0.00570, 0.00572], [L"5.60", L"5.62", L"5.64", L"5.66", L"5.68", L"5.70", L"5.72"]),
               ytick=([0.09250, 0.09275, 0.09300, 0.09325, 0.09350], [L"0.09250", L"0.09275", L"0.09300", L"0.09325", L"0.09350"]),
               ylims=(yZoomDispBott-0.00005, 0.09355), xlims=(5.6e-3, tZoomDispRight),
               bottom_margin=6mm, left_margin=30mm, right_margin=10mm, top_margin=3mm, size=(900, 600))

    # flowpipe
    plot!(fig, sol_PosTwoThirds(tZoomDispLeft .. tZoomDispRight), vars=(0, 1), c=:lightblue, lc=:black, lw=0.5, alpha=1., lab=L"\textrm{Flowpipe}")

    # solucion analitica
    tdomasol_zoom = range(tZoomDispLeft, tZoomDispRight, length=15_000)

    asol_PosTwoThirds_zoom, asol_VelTwoThirds_zoom = _sol_analitica_clamped(PosTwoThirds, tdomasol_zoom, 2000);

    plot!(fig, tdomasol_zoom, asol_PosTwoThirds_zoom, c=:magenta, lw=1.5, lab=L"\textrm{Analytic (PDE)}")

    paso = 8

    ai = idxtleft
    bi = idxtright

    # == Newmark ==
    # continua sin markers
    plot!(fig, tdom[ai:bi], uNewmark_PosTwoThirds[ai:bi], seriestype=:path, marker=:none, lab="", color=:red, lw=2.0)

    # con markers pero salteados
    plot!(fig, tdom[ai:paso:bi], uNewmark_PosTwoThirds[ai:paso:bi], linetype=:scatter, marker=:circle, markersize=5, lab=L"\textrm{Newmark}", color=:red)

    # == Bathe ==
    # continua sin markers
    plot!(fig, tdom[ai:bi], uBathe_PosTwoThirds[ai:bi], seriestype=:path, marker=:none, lab="", color=:green, lw=2.0)

    # con markers pero salteados
    plot!(fig, tdom[ai:paso:bi], uBathe_PosTwoThirds[ai:paso:bi], linetype=:scatter, marker=:utriangle, markersize=5, lab=L"\textrm{Bathe}", color=:green)

    # == Analytic ODE ==
    plot!(fig, sol_nobloating, vars=(0, PosTwoThirds), seriestype=:path, marker=:none, lw = 1.5, alpha=1, c=:blue, lab=L"\textrm{Analytic (ODE)}")
    plot!(fig, sol_nobloating[ai:paso:bi], vars=(0, PosTwoThirds), linetype=:scatter, marker=:square, markersize=5,  c=:blue, lab="")

end

fig = plot_displacements_zoom()
savefig(fig, joinpath(TARGET_FOLDER, "u_node700_zoom.pdf"))

# =========================
# velocidad
# =========================

function plot_velocity()

    fig = plot(xlab=L"\textrm{Time}", ylab=L"\textrm{Velocity of node } %$PosTwoThirds",
               legend=:bottomright,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(20, "Times"),
               yguidefont=font(20, "Times"),
               xtick=([0, 0.0025, 0.005, 0.0075, 0.01],[L"0.000", L"0.0025", L"0.0050", L"0.0075", L"0.01"]),
               ytick=([-75, -50, -25, 0, 25, 50, 75], [L"-75", L"-50", L"0.0", L"25", L"50", L"75"]),
               xlims=(0.0, 0.012), # , ylims=(-90.0, 90.0),
               bottom_margin=6mm, left_margin=20mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    # flowpipe
    plot!(fig, sol_VelTwoThirds, vars=(0, 1), c=:lightblue, lc=:lightblue, lw=4.0, alpha=1., lab="") # lab=L"\textrm{Flowpipe}")

    # plot!(fig, tdom, uNewmark, seriestype=:path, marker=:none, color=:red, lab=L"\textrm{Newmark}", lw=1.0)

    # plot!(fig, tdom, uBathe, seriestype=:path, marker=:none, color=:green, lab=L"\textrm{Bathe}", lw=1.0)
    # plot!(fig, sol_nobloating, vars=(0, PosTwoThirds), seriestype=:path, marker=:none, lw = 2.0, c=:blue, lab=L"\textrm{Analytic (ODE)}")

    tZoomDispLeft  = 0.008175
    tZoomDispRight = 0.008375

    yZoomDispBott  = 40
    yZoomDispTop = 92

    # == Analytic ODE ==
    plot!(fig, sol_nobloating, vars=(0, VelTwoThirds), seriestype=:path, marker=:none, lw = 1.5, alpha=1, c=:blue, lab="")
    #plot!(fig, sol_nobloating[ai:paso:bi], vars=(0, VelTwoThirds), linetype=:scatter, marker=:square, markersize=5,  c=:blue, lab=L"\textrm{Analytic (ODE)}")

    # == Analytic PDE ==
    tdomasol = range(0, finalTime, length=1000)

    _, asol_VelTwoThirds = _sol_analitica_clamped(PosTwoThirds, tdomasol, 1000);

    plot!(fig, tdomasol, asol_VelTwoThirds, c=:magenta, lw=1.5, lab="") # , lab=L"\textrm{Analytic (PDE)}")

    # ZOOM BOX
    PBotLeft  = [tZoomDispLeft,yZoomDispBott]
    PBotRight = [tZoomDispRight,yZoomDispBott]
    PTopLeft  = [tZoomDispLeft,yZoomDispTop]
    PTopRight = [tZoomDispRight,yZoomDispTop]

    L1 = LineSegment( PBotLeft, PBotRight )
    L2 = LineSegment( PBotRight, PTopRight )
    L3 = LineSegment( PTopRight, PTopLeft )
    L4 = LineSegment( PTopLeft, PBotLeft )

    plot!(fig, L1, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L2, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L3, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L4, color=:black, style=:dash, lw=3.0, seriestype=:solid, marker=:none, alpha=1.)

end

fig = plot_velocity()
savefig(fig, joinpath(TARGET_FOLDER, "v_node700.pdf"))

# ==============================
# graficos velocidad (ZOOM)
# ==============================

function plot_velocity_zoom()

    tZoomDispLeft  = 0.008175
    tZoomDispRight = 0.008375

    yZoomDispBott  = 30
    yZoomDispTop = 92

    idxtleft = findfirst(t -> t > tZoomDispLeft, tdom)
    idxtright = findfirst(t -> t > tZoomDispRight, tdom)

    fig = plot(xlab=L"\textrm{Time } (\times 10^{-3})", ylab=L"\textrm{Velocity of node } %$PosTwoThirds",
               legend=:bottomright,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(20, "Times"),
               yguidefont=font(20, "Times"),
               xtick=([0.00820, 0.00825, 0.00830, 0.00835], [L"8.20", L"8.25", L"8.30", L"8.35"]),
               ytick=([30, 40, 50, 60, 70, 80, 90], [L"30", L"40", L"50", L"60", L"70", L"80", L"90"]),
               ylims=(yZoomDispBott, yZoomDispTop),
               xlims=(tZoomDispLeft, tZoomDispRight),
               bottom_margin=6mm, left_margin=30mm, right_margin=10mm, top_margin=3mm, size=(900, 600))

    # flowpipe
    plot!(fig, sol_VelTwoThirds(tZoomDispLeft .. tZoomDispRight), vars=(0, 1), c=:lightblue, lc=:black, lw=0.5, alpha=1., lab=L"\textrm{Flowpipe}")

    # solucion analitica
    tdomasol_zoom = range(tZoomDispLeft, tZoomDispRight, length=103)

    asol_VelTwoThirds_zoom, asol_VelTwoThirds_zoom = _sol_analitica_clamped(PosTwoThirds, tdomasol_zoom, 1000);

    plot!(fig, tdomasol_zoom, asol_VelTwoThirds_zoom, c=:magenta, lw=1.5, lab=L"\textrm{Analytic (PDE)}")

    paso = 8

    ai = idxtleft
    bi = idxtright

    # == Newmark ==
    # continua sin markers
    plot!(fig, tdom[ai:bi], uNewmark_VelTwoThirds[ai:bi], seriestype=:path, marker=:none, lab="", color=:red, lw=2.0)

    # con markers pero salteados
    plot!(fig, tdom[ai:paso:bi], uNewmark_VelTwoThirds[ai:paso:bi], linetype=:scatter, marker=:circle, markersize=5, lab=L"\textrm{Newmark}", color=:red)

    # == Bathe ==
    # continua sin markers
    plot!(fig, tdom[ai:bi], uBathe_VelTwoThirds[ai:bi], seriestype=:path, marker=:none, lab="", color=:green, lw=2.0)

    # con markers pero salteados
    plot!(fig, tdom[ai:paso:bi], uBathe_VelTwoThirds[ai:paso:bi], linetype=:scatter, marker=:utriangle, markersize=5, lab=L"\textrm{Bathe}", color=:green)

    # == Analytic ODE ==
    plot!(fig, sol_nobloating, vars=(0, VelTwoThirds), seriestype=:path, marker=:none, lw = 1.5, alpha=1, c=:blue, lab=L"\textrm{Analytic (ODE)}")
    plot!(fig, sol_nobloating[ai:paso:bi], vars=(0, VelTwoThirds), linetype=:scatter, marker=:square, markersize=5,  c=:blue, lab="")
end

fig = plot_velocity_zoom()
savefig(fig, joinpath(TARGET_FOLDER, "v_node700_zoom.pdf"))

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
using Roots

struct fcOutput{T}
    Lbound::T
    Rbound::T
    Rmulti::T
    Req::Vector{T}
    Rbound_n::Vector{T}
    Rmulti_n::Vector{T}
end

function polyfc(L0::Real, Kₓ::Real, f::Number, Rtot::Vector, IgGC::Vector, Kav::AbstractMatrix)
    # Data consistency check
    ansType = promote_type(typeof(L0), typeof(Kₓ), typeof(f), eltype(Rtot), eltype(IgGC), eltype(Kav))
    @assert size(Kav) == (length(IgGC), length(Rtot))
    @assert ndims(Kav) == 2
    IgGC /= sum(IgGC)

    # Setup constant terms
    if length(IgGC) > 1
        A = vec(Kav' * IgGC)
    else
        A = vec(Kav)
    end
    L0fK = L0 * f / Kₓ
    L0fA = L0 * f * A

    # Solve for Phisum
    function phi_func(ϕs::Real)
        return ϕs - Kₓ * dot(A, Rtot ./ (1.0 .+ L0fA .* (1 .+ ϕs) .^ (f .- 1)))
    end

    bounds = (convert(ansType, 0.0), convert(ansType, Kₓ * dot(A, Rtot)))
    Phisum = find_zero(phi_func, bounds, Bisection())

    Req = Rtot ./ (1.0 .+ L0fA * (1 + Phisum) ^ (f - 1))
    Phisum_n = sum(ones(ansType, size(Kav, 1), size(Kav, 2)) .* IgGC .* Kav .* Req' .* Kₓ, dims = 1)

    w = fcOutput{ansType}(
        L0 / Kₓ * ((1 + Phisum)^f - 1),
        sum(Rtot - Req),
        L0fK * Phisum * ((1 + Phisum)^(f - 1) - 1),
        Req,
        vec(L0fK .* Phisum_n * (1 + Phisum)^(f - 1)),
        vec(L0fK .* Phisum_n * ((1 + Phisum)^(f - 1) - 1)),
    )
    return w
end

polyfcm = (Kₓ, f, Rtot, IgG, Kav) -> polyfc(sum(IgG) / f, Kₓ, f, Rtot, IgG ./ sum(IgG), Kav)


text:
module Rasters

# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end Rasters

using Dates

import Adapt,
       ArchGDAL,
       ColorTypes,
       ConstructionBase,
       DiskArrays,
       FillArrays,
       Flatten,
       GeoInterface,
       HDF5,
       PolygonInbounds,
       ProgressMeter,
       Missings,
       Mmap,
       NCDatasets,
       RecipesBase,
       Reexport,
       Setfield

Reexport.@reexport using DimensionalData, GeoFormatTypes, RasterDataSources

using DimensionalData.Tables,
      DimensionalData.LookupArrays,
      DimensionalData.Dimensions
      DimensionalData.LookupArrays.IntervalSets

using DimensionalData: Name, NoName
using .Dimensions: StandardIndices, DimTuple
using .LookupArrays: LookupArrayTuple

using RecipesBase: @recipe, @series
using Base: tail, @propagate_inbounds

using Setfield: @set, @set!
using ColorTypes: RGB

export AbstractRaster, Raster
export AbstractRasterStack, RasterStack
export AbstractRasterSeries, RasterSeries
export Projected, Mapped
export Band
export missingval, boolmask, missingmask, replace_missing, replace_missing!,
       aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!,
       resample, warp, crop, extend, trim, slice, points, subset, inpolygon,
       classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!,
       setcrs, setmappedcrs
export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
export reproject, convertlookup
export geoarray, stack, series

const DD = DimensionalData
const DA = DiskArrays
const GI = GeoInterface
const LA = LookupArrays

# DimensionalData documentation urls
const DDdocs = "https://rafaqz.github.io/DimensionalData.jl/stable/api"
const DDdimdocs = joinpath(DDdocs, "#DimensionalData.Dimension")
const DDarraydocs = joinpath(DDdocs, "#DimensionalData.AbstractDimensionalArray")
const DDabssampleddocs = joinpath(DDdocs, "#DimensionalData.AbstractSampled")
const DDsampleddocs = joinpath(DDdocs, "#DimensionalData.Sampled")
const DDlocusdocs = joinpath(DDdocs, "#DimensionalData.Locus")
const DDselectordocs = joinpath(DDdocs, "#DimensionalData.Selector")
const DDtidocs = joinpath(DDdocs, "#DimensionalData.Ti")

const EXPERIMENTAL = """
    WARNING: This feature is experimental. It may change in future versions, and may
    not be 100% reliable in all cases. Please file github issues if problems occur.
    """

# Source dispatch singletons
struct NCDfile end
struct GRDfile end
struct GDALfile end
struct SMAPfile end

include("lookup.jl")
include("dimensions.jl")
include("filearray.jl")
include("filestack.jl")
include("openstack.jl")
include("array.jl")
include("stack.jl")
include("series.jl")

const RasterStackOrArray = Union{AbstractRasterStack,AbstractRaster}
const RasterSeriesOrStack = Union{AbstractRasterSeries,AbstractRasterStack}

include("utils.jl")
include("skipmissing.jl")
include("polygon_ops.jl")
include("table_ops.jl")
include("create.jl")
include("read.jl")
include("write.jl")
include("convenience.jl")
include("show.jl")
include("plotrecipes.jl")

include("methods/aggregate.jl")
include("methods/classify.jl")
include("methods/crop_extend.jl")
include("methods/extract.jl")
include("methods/inpolygon.jl")
include("methods/mask.jl")
include("methods/mosaic.jl")
include("methods/points.jl")
include("methods/rasterize.jl")
include("methods/replace_missing.jl")
include("methods/reproject.jl")
include("methods/resample.jl")
include("methods/slice_combine.jl")
include("methods/trim.jl")
include("methods/warp.jl")

include("sources/grd.jl")
include("sources/smap.jl")
include("sources/ncdatasets.jl")
include("sources/gdal.jl")
include("sources/rasterdatasources.jl")

end


text:
using JSON
using AssetRegistry
using Sockets
using Base64: stringmime
export webio_serve

"""
    webio_serve(app, port=8000)

Serve a Mux app which might return a WebIO node.
"""
function webio_serve(app, args...)
    http = Mux.App(Mux.mux(
        Mux.defaults,
        app,
        Mux.notfound()
    ))

    websock = Mux.App(Mux.mux(
        Mux.wdefaults,
        Mux.route("/webio-socket", create_socket),
        Mux.wclose,
        Mux.notfound(),
    ))

    Mux.serve(http, websock, args...)
end

struct WebSockConnection <: WebIO.AbstractConnection
    sock
end

function create_socket(req)
    sock = req[:socket]
    conn = WebSockConnection(sock)

    t = @async while isopen(sock)
        data = read(sock)

        msg = JSON.parse(String(data))
        WebIO.dispatch(conn, msg)
    end

    wait(t)
end

function Sockets.send(p::WebSockConnection, data)
    write(p.sock, sprint(io->JSON.print(io,data)))
end

Base.isopen(p::WebSockConnection) = isopen(p.sock)

Mux.Response(o::AbstractWidget) = Mux.Response(Widgets.render(o))
function Mux.Response(o::Union{Node, Scope})
    key = AssetRegistry.register(joinpath(@__DIR__, "..", "..", "packages/mux-provider/dist"))
    Mux.Response(
        """
        <!doctype html>
        <html>
          <head>
            <meta charset="UTF-8">
            <script src="$(WebIO.baseurl[])$key/mux.js"></script>
          </head>
          <body>
            $(stringmime(MIME"text/html"(), o))
          </body>
        </html>
        """
    )
end

function WebIO.register_renderable(::Type{T}, ::Val{:mux}) where {T}
    eval(:(Mux.Response(x::$T) = Mux.Response(WebIO.render(x))))
end

WebIO.setup_provider(::Val{:mux}) = nothing # Mux setup has no side-effects
WebIO.setup(:mux)


text:
export stopAtIteration

function stopAtIteration(n)
    function(m)
        if m.iterationCount > n
            m.done = true
        end
    end
end

export stopAfter

function stopAfter(Δt)
    stopTime = time() + Δt
    function(m)
        if time() > stopTime
            m.done = true
        end
    end
end

export fixedEvidence

function fixedEvidence()
    logEv0 = -Inf
    function(m)
        logEv = logEvidence(m)

        if logEv == logEv0
            m.done = true
        end
        logEv0 = logEv
        return logEv
    end
end


text:
unpack(v::StaticVector{2}) = @inbounds return v[1], v[2]

@inline function cross2(v1::StaticVector{2}, v2::StaticVector{2})
    x1, y1 = unpack(v1)
    x2, y2 = unpack(v2)
    x1 * y2 - x2 * y1
end

normsquared(x) = sum(y -> y^2, x)


text:
using MPI
using PencilArrays
using Test

MPI.Initialized() || MPI.Init()

comm = MPI.COMM_WORLD
nprocs = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
myid = rank + 1

# This can be simplified to `redirect_stdout(devnull)` on Julia ≥ 1.6.
let dev_null = @static Sys.iswindows() ? "nul" : "/dev/null"
    MPI.Comm_rank(comm) == 0 || redirect_stdout(open(dev_null, "w"))
end

pen = Pencil((16, 32, 14), comm)
u = PencilArray{Int32}(undef, pen)
fill!(u, 2myid)

@testset "Reductions" begin
    @test minimum(u) == 2
    @test maximum(u) == 2nprocs
    @test minimum(abs2, u) == 2^2
    @test maximum(abs2, u) == (2nprocs)^2
    @test sum(u) === MPI.Allreduce(sum(parent(u)), +, comm)
    @test sum(abs2, u) === MPI.Allreduce(sum(abs2, parent(u)), +, comm)

    # These exact equalities work because we're using integers.
    # They are not guaranteed to work with floats.
    @test foldl(min, u) === minimum(u)
    @test mapfoldl(abs2, min, u) === minimum(abs2, u)

    @test foldr(min, u) === minimum(u)
    @test mapfoldr(abs2, min, u) === minimum(abs2, u)
end


text:
function random_сsrand(line::Ref{UInt}, ns::Namespace, seed::AetherlangObject)::AetherlangObject{AetherlangVoid}
    ccall(:srand, Cvoid, (Int32,), seed.dataref)
    AetherlangObject()
end

function AE_choose(line::Ref{UInt}, ns::Namespace, arr::AetherlangObject)::AetherlangObject
    AetherlangObject(arr.dataref[rand(1:length(arr.dataref))])
end

function AE_shuffled(line::Ref{UInt}, ns::Namespace, arr::AetherlangObject)::AetherlangObject
    a = [arr.dataref[i] for i in 1:length(arr.dataref)]
    for i in 1:length(a)
        q::Int = rand(1:length(a))
        a[i], a[q] = a[q], a[i]
    end
    AetherlangObject(a)
end

random_namespace_modify = Namespace(
    "c.srand" => AetherlangObject(random_сsrand),
    "c.rand" => AetherlangObject((line::Ref{UInt}, ns::Namespace, dummy::AetherlangObject)->AetherlangObject(ccall(:rand, Int32, ()))),
    "choose" => AetherlangObject(AE_choose),
    "shuffled" => AetherlangObject(AE_shuffled)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
# we are making a Julia model if this gets called, so lets load the Julia strategy -
include("./strategy/General.jl")
include("./strategy/JuliaStrategy.jl")

function make_julia_model(path_to_model_file::String, path_to_output_dir::String;
    defaults_file_name::String="Defaults.toml", model_type::Symbol=:static)

    # start -
    time_start = time();

    # Fix: Need to generate a Default.toml project first #2
    # check, do we have the Defaults file?
    path_to_defaults_file = "$(dirname(path_to_model_file))/$(defaults_file_name)"
    if (isfile(path_to_defaults_file) == false)
        throw(ArgumentError("Ooops! We are missing the Defaults.toml file. Please generate a default project using the generate_default_project function. For details see the generate_default_project help in the REPL"))
    end

    @info "Loading the $(defaults_file_name)"

    # Load the defaults file -
    tmp_dict =  TOML.parsefile(path_to_defaults_file)
    default_parameter_dictionary = tmp_dict["default_parameter_values"]
    author_metadata_dictionary = tmp_dict["author_metadata"]

    # legit path to model file?
    is_file_path_ok(path_to_model_file)

    # before we get too far along, we need to check if the user already has code in the location that we want to generate at -
    # if they do, then move it -
    if (isdir(path_to_output_dir) == true)

        @info "Ooops! We found some code where you wanted to generate your model. Moving ..."

        # ok, looks like we may have a conflict - mv the offending code
        if (move_existing_project_at_path(path_to_output_dir) == false)

            # Something happend ... the world is ending ...
            throw(ArgumentError("automatic directory conflict resolution failed. Unable to move existing directory $(path_to_output_dir)"))
        end
    end

    # let the uset know what is going on -
    @info "Starting a new code genration job: reading model file at $(path_to_model_file)"

    # initialize -
    src_component_set = Set{VLProgramComponent}()
    config_component_set = Set{VLProgramComponent}()
    root_component_set = Set{VLProgramComponent}()

    # Load the JSON configuration file -
    path_to_configuration_file = "$(path_to_package)/configuration/Configuration.json"
    config_dict = JSON.parsefile(path_to_configuration_file)

    # let the uset know what is going on -
    @info "Loaded the configuration file $(path_to_configuration_file)"

    # parse the model file -
    problem_dictionary = parse_vff_model_file(path_to_model_file, config_dict, default_parameter_dictionary)
    problem_dictionary["configuration_dictionary"] = config_dict
    problem_dictionary["default_parameter_dictionary"] = default_parameter_dictionary
    problem_dictionary["author_metadata_dictionary"] = author_metadata_dictionary

    # let the uset know what is going on -
    @info "Model file parsing complete"

    # generate a parameters toml file -
    (program_component_parameter_dictionary, runtime_model_parameter_name_array) = generate_control_parameters_dictionary_component(problem_dictionary)
    push!(config_component_set, program_component_parameter_dictionary)
    problem_dictionary["runtime_model_parameter_name_array"] = runtime_model_parameter_name_array

    @info "Control function generation complete"

    # generate the data dictionary component -
    program_component_data_dictionary = build_data_dictionary_program_component(problem_dictionary)
    push!(src_component_set, program_component_data_dictionary)

    @info "Data dictionary generation complete"

    # generate the kinetics component -
    program_component_kinetics = build_kinetics_program_component(problem_dictionary)
    push!(src_component_set, program_component_kinetics)

    @info "Kinetics function generation complete"

    # Write the Control functions -
    program_component_control = build_control_program_component(problem_dictionary)
    push!(src_component_set,program_component_control)

    @info "Control function generation complete ..."

    # Write the stoichiometric_matrix --
    program_component_stoichiometric_matrix = generate_stoichiomteric_matrix_program_component(problem_dictionary)
    push!(src_component_set,program_component_stoichiometric_matrix)

    @info "The stoichiometric matrix has been generated  ... writing model files to $(path_to_output_dir)"

    # Transfer distrubtion jl files to the output -> these files are shared between model types
    transfer_distribution_files("$(path_to_package)/distribution/julia/src", "$(path_to_output_dir)/src",".jl")

    # Transfer root jl files -> these files are moved to the root dir of the file
    transfer_distribution_files("$(path_to_package)/distribution/julia/root", "$(path_to_output_dir)",".jl")

    # Transfer dynamic or static jl files -
    if (model_type == :static)

        # Transfer root jl files -> these files are moved to the root dir of the project -
        transfer_distribution_files("$(path_to_package)/distribution/julia/root/static", "$(path_to_output_dir)",".jl")
    elseif (model_type == :dynamic)
        # Transfer root jl files -> these files are moved to the root dir of the project -
        transfer_distribution_files("$(path_to_package)/distribution/julia/root/dynamic", "$(path_to_output_dir)",".jl")
    else
        # Something happend ... the world is ending ...
        throw(ArgumentError("Oooops! Unsupported model type has been requested: $(model_type). Model type must be either static or dynamic."))
    end

    # Transfer distibution TOML files to output -
    transfer_distribution_files("$(path_to_package)/distribution","$(path_to_output_dir)/src/config",".toml")

    # transfer the README files -
    path_to_readme_file = splitdir(path_to_model_file)[1]
    transfer_distribution_files("$(path_to_package)/distribution",path_to_readme_file,".md")

    # dump src and config components to disk -
    write_program_components_to_disk("$(path_to_output_dir)/src", src_component_set)
    write_program_components_to_disk("$(path_to_output_dir)/src/config", config_component_set)

    # stop timer -
    elapsed_time = (time() - time_start)    # in ns
    problem_dictionary["elapsed_time"] = elapsed_time
    problem_dictionary["path_to_model_file"] = path_to_model_file
    problem_dictionary["path_to_output_dir"] = path_to_output_dir

    # generate model report file (useful for debugging) -
    program_component_report = build_model_generation_report(problem_dictionary)
    push!(root_component_set, program_component_report)

    # Dump the component_set to disk -
    write_program_components_to_disk("$(path_to_output_dir)", root_component_set)
    @info "Done! elapsed_time: $(time() - time_start)s"
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
"""
    AVLTree

struct
"""
mutable struct AVLTree{K,D}
    root::Union{Node{K,D},Nothing}
end

AVLTree() = AVLTree{Any,Any}(nothing)
AVLTree{K,D}() where {K,D} = AVLTree{K,D}(nothing)

Base.eltype(::Type{AVLTree{K,D}}) where {K,D} = Tuple{K,D}
Base.getindex(tr::AVLTree{K,D},k::K) where {K,D} = Base.getkey(tr, k)
Base.setindex!(tr::AVLTree{K,D},k::K,d::D) where {K,D} = AVLTrees.insert!(tr, k, d)
Base.haskey(tr::AVLTree{K,D},k::K) where {K,D} = !(find_node(tr, k) === nothing)
Base.length(tr::AVLTree{K,D}) where {K,D} = AVLTrees.size(tr)
Base.isempty(tr::AVLTree{K,D}) where {K,D} = isnothing(tr.root)

function Base.getkey(tr::AVLTree{K,D}, k::K) where {K,D}
    d = findkey(tr, k)
    if isnothing(d) throw(KeyError(k)) else d end
end

function Base.size(tree::AVLTree)
    return __size(tree.root)
end # function

@inline function __size(node::Union{Nothing,Node})
    if isnothing(node)
        return 0
    end
    return __size(node.left) + __size(node.right) + 1
end

"""
    insert!(args)

documentation
"""
function Base.insert!(tree::AVLTree{K,D}, key, data) where {K,D}
    parent = nothing
    node = tree.root

    while !isnothing(node)
        parent = node
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            node.data = data
            return
        end
    end

    if isnothing(parent)
        tree.root = Node{K,D}(key, data)
    elseif key < parent.key
        parent.left = Node{K,D}(key, data, parent)
        balance_insertion(tree, parent, true)
    elseif key > parent.key
        parent.right = Node{K,D}(key, data, parent)
        balance_insertion(tree, parent, false)
    end

    return
end # function

"""
    balance_insertion(tree::AVLTree{K,D},node::Node{K,D},left_insertion::Bool) where {K,D}

documentation
"""
function balance_insertion(
    tree::AVLTree{K,D},
    node::Node{K,D},
    left_insertion::Bool,
) where {K,D}
    while !isnothing(node)
        node.bf += ifelse(left_insertion, -1, 1)
        node, height_changed = rebalance(tree, node)
        if height_changed
            break
        end
        left_insertion = !isnothing(node.parent) && node.parent.left == node
        node = node.parent
    end
end # function

"""
    rebalance(node::Node)

documentation
"""
function rebalance(tree::AVLTree{K,D}, node::Node{K,D})::Tuple{Node{K,D},Bool,} where {K,D}
    if node.bf == 2
        height_changed = node.right.bf != 0
        if node.right.bf == -1
            rotate_right(tree, node.right)
        end
        node = rotate_left(tree, node)
        return node, height_changed
    elseif node.bf == -2
        height_changed = node.left.bf != 0
        if node.left.bf == 1
            rotate_left(tree, node.left)
        end
        node = rotate_right(tree, node)
        return node, height_changed
    else
        return node, node.bf == 0
    end
end # function

@inline function rotate_left(t::AVLTree{K,D}, x::Node{K,D}) where {K,D}
    y = x.right

    x.right = y.left
    if !isnothing(y.left)
        y.left.parent = x
    end
    y.left = x

    if isnothing(x.parent)
        t.root = y
    else
        if x.parent.left == x
            x.parent.left = y
        else
            x.parent.right = y
        end
    end

    y.parent = x.parent
    x.parent = y

    x.bf -= y.bf * (y.bf >= 0) + 1
    y.bf += x.bf * (x.bf < 0) - 1

    return y
end

@inline function rotate_right(t::AVLTree{K,D}, x::Node{K,D}) where {K,D}
    y = x.left

    x.left = y.right
    if !isnothing(y.right)
        y.right.parent = x
    end
    y.right = x

    if isnothing(x.parent)
        t.root = y
    else
        if x.parent.left == x
            x.parent.left = y
        else
            x.parent.right = y
        end
    end

    y.parent = x.parent
    x.parent = y

    x.bf -= y.bf * (y.bf < 0) - 1
    y.bf += x.bf * (x.bf >= 0) + 1

    return y
end

"""
    delete!(tree::AVLTree{K,D}, node::Node{K,D}) where {K,D}

documentation
"""
function Base.delete!(tree::AVLTree{K,D}, node::Node{K,D}) where {K,D}
    if !isnothing(node.left)
        if !isnothing(node.right)
            # left != nothing && right != nothing
            temp = node.right
            while !isnothing(temp.left)
                temp = temp.left
            end
            # switch spots completely
            node.key = temp.key
            node.data = temp.data
            delete!(tree, temp)
        else
            # left != nothing && right == nothing
            dir = __parent_replace(tree, node, node.left)
            balance_deletion(tree, node.parent, dir)
        end
    else
        if !isnothing(node.right)
            # left == nothing && right != nothing
            dir = __parent_replace(tree, node, node.right)
            balance_deletion(tree, node.parent, dir)
        else
            # left == nothing && right == nothing
            dir = __parent_replace(tree, node, nothing)
            balance_deletion(tree, node.parent, dir)
        end
    end
    return
end # function

"""
    delete!(tree::AVLTree{K,D}, key::K) where {K,D}
"""
function Base.delete!(tree::AVLTree{K,D}, key::K) where {K,D}
    node = find_node(tree, key)
    if !isnothing(node)
        delete!(tree, node)
    end
end # function

"""
    balance_deletion(args)
"""
function balance_deletion(
    tree::AVLTree{K,D},
    node::Union{Node{K,D},Nothing},
    left_delete::Union{Nothing,Bool},
) where {K,D}
    while !isnothing(node)
        node.bf += ifelse(left_delete, 1, -1)
        node, height_changed = rebalance(tree, node)
        if !height_changed
            break
        end
        left_delete = !isnothing(node.parent) && node.parent.left == node
        node = node.parent
    end
end # function

#    __parent_replace(tree::AVLTree{K,D}, node::Node{K,D}, replacement::Node{K,D})
#
# Replaces node with its only child. Used on nodes with a single child when erasing a node.
@inline function __parent_replace(
    tree::AVLTree{K,D},
    node::Node{K,D},
    replacement::Node{K,D},
) where {K,D}
    if !isnothing(node.parent)
        replacement.parent = node.parent
        if node.parent.right == node
            node.parent.right = replacement
            return false
        else
            node.parent.left = replacement
            return true
        end
    else
        replacement.parent = nothing
        tree.root = replacement
        return nothing
    end
end # function

#    __parent_replace(tree::AVLTree{K,D}, node::Node{K,D}, replacement::Nothing)
# Replaces node with nothing. Used on leaf nodes when erasing a node.
@inline function __parent_replace(
    tree::AVLTree{K,D},
    node::Node{K,D},
    replacement::Nothing,
) where {K,D}
    if !isnothing(node.parent)
        if node.parent.right == node
            node.parent.right = replacement
            return false
        else
            node.parent.left = replacement
            return true
        end
    else
        tree.root = replacement
        return nothing
    end
end # function

"""
    find(tree::AVLTree{K,D}, key::K) where {K,D}

    Warning: do not use it to check whether `key` is in the `tree`.
    It returns the node.data if found which can be `nothing`.
"""
function findkey(tree::AVLTree{K,D}, key::K) where {K,D}
    node = tree.root
    while !isnothing(node)
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            return node.data
        end
    end
    return nothing
end # function

"""
    find_node(args)
"""
function find_node(tree::AVLTree{K,D}, key::K) where {K,D}
    node = tree.root
    while !isnothing(node)
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            return node
        end
    end
    return nothing
end # function

# Iteration interface

function Base.iterate(tree::AVLTree)
    if isnothing(tree.root)
        return nothing
    end
    node = tree.root
    while !isnothing(node.left)
        node = node.left
    end
    return (node.key, node.data), node
end

function Base.iterate(tree::AVLTree, node::Node)
    if !isnothing(node.right)
        node = node.right
        while !isnothing(node.left)
            node = node.left
        end
    else
        prev = node
        while !isnothing(node) && node.left != prev
            prev = node
            node = node.parent
        end
    end

    if isnothing(node)
        return nothing
    end

    return (node.key, node.data), node
end

# Pop and get methods

function Base.popfirst!(tree::AVLTree)
    # traverse to left-most node
    if isnothing(tree.root)
        return
    end
    node = tree.root
    while !isnothing(node.left)
        node = node.left
    end
    # delete node and return data
    node_data = node.data
    delete!(tree, node)
    return node_data
end

function Base.pop!(tree::AVLTree{K,D}, key::K) where {K,D}
    node = AVLTrees.find_node(tree, key)
    if !isnothing(node)
        node_dat = node.data
        delete!(tree, node)
        return node_dat
    else
        return
    end
end

function Base.firstindex(tree::AVLTree)
    # traverse to left-most node
    if isnothing(tree.root)
        return
    end
    node = tree.root
    while !isnothing(node.left)
        node = node.left
    end
    # return node key
    return node.key
end

## Print and Show methods

function Base.print(io::IO, tree::AVLTree{K,D}) where {K,D}
    str_lst = Vector{String}()
    for (k, v) in Base.Iterators.take(tree, 10)
        push!(str_lst, "$k => $v")
    end
    print(io, "AVLTree{$K,$D}(")
    print(io, join(str_lst, ", "))
    length(str_lst) == 10 && print(io, ", ⋯ ")
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", tree::AVLTree{K,D}) where {K,D}
    str_lst = Vector{String}()
    indent_str = "  "
    for (k, v) in Base.Iterators.take(tree, 10)
        push!(str_lst, indent_str * "$k => $v")
    end
    if length(str_lst) > 0
        print(io, "AVLTree{$K,$D} with $(length(str_lst)) entries:\n")
        print(io, join(str_lst, "\n"))
    else
        print(io, "AVLTree{$K,$D}()")
    end
    length(str_lst) == 10 && print(io, "\n", indent_str * "⋮ => ⋮ \n")
end


text:
"""
    isleaf(node::Node)
Test whether a node is a leaf or not.
"""
isleaf(node::Node) = node.children === nothing

"""
    isroot(node::Node)
Return `true` if `node` is the root node (meaning, it has no parent).
"""
isroot(node::Node) = node.parent === nothing

"""
    children(node::Node)

Return the immediate children of `node`.
"""
function children(node::Node)
    isleaf(node) ? Vector{Node}() : collect(values(node.children))
end

"""
    ordered_children(node)

Return the children as an array, ordered first by "+"
"""
function ordered_children(node)
    links_chnodes = Array{Node,1}()
    for (name, chnode) in node.children
        if chnode.MTG.link == "+"
            pushfirst!(links_chnodes, chnode)
        else
            push!(links_chnodes, chnode)
        end
    end
    return links_chnodes
end

"""
    lastchild(node::Node)

Get the last child of `node`, or `nothing` if the node is a leaf.

"""
function lastchild(node::Node)
    if isleaf(node)
        return nothing
    else
        allchildren = node.children
        return allchildren[maximum(keys(allchildren))]
    end
end

"""
Add a new child to a parent node, and add the parent node as the parent.

See also [`insert_child!`](@ref) for a more user-friendly approach.
"""
function addchild!(parent::Node, id::Int, MTG::M, attributes) where {M<:AbstractNodeMTG}
    child = Node(id, parent, MTG, attributes)
    addchild!(parent, child)
end

function addchild!(parent::Node, id::Int, MTG::M) where {M<:AbstractNodeMTG}
    child = Node(id, parent, MTG)
    addchild!(parent, child)
end

function addchild!(parent::Node, child::Node; force = false)

    if child.parent === missing || force == true
        child.parent = parent
    elseif child.parent != parent && force == false
        error("The node already has a parent. Hint: use `force=true` if needed.")
    end

    if parent.children === nothing
        parent.children = Dict{Int,Node}(child.id => child)
    else
        push!(parent.children, child.id => child)
    end
end

"""
Find the root node of a tree, given any node in the tree.
"""
function get_root(node::Node)
    if isroot(node)
        return (node)
    else
        get_root(node.parent)
    end
end

"""
    siblings(node::Node)

Return the siblings of `node` as a vector of nodes (or `nothing` if non-existant).
"""
function siblings(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent === nothing && return nothing

    all_siblings = children(node.parent)

    all_siblings[findall(x -> x != node, all_siblings)]
end

"""
    nextsibling(node::Node)

Return the next sibling of `node` (or `nothing` if non-existant).
"""
function nextsibling(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent === nothing && return nothing

    all_siblings = children(node.parent)
    # Get the index of the current node in the siblings:
    node_index = findfirst(x -> x == node, all_siblings)
    if node_index < length(all_siblings)
        all_siblings[node_index+1]
    else
        nothing
    end
end

"""
    prevsibling(node::Node)

Return the previous sibling of `node` (or `nothing` if non-existant).
"""
function prevsibling(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent === nothing && return nothing

    all_siblings = children(node.parent)
    # Get the index of the current node in the siblings:
    node_index = findfirst(x -> x == node, all_siblings)
    if node_index > 0
        all_siblings[node_index-1]
    else
        nothing
    end
end

"""
    lastsibling(node::Node)

Return the last sibling of `node` (or `nothing` if non-existant).
"""
function lastsibling(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent === nothing && return nothing

    all_siblings = children(node.parent)
    # Get the index of the current node in the siblings:

    all_siblings[maximum(keys(all_siblings))]
end

"""
    new_id(mtg)
    new_id(mtg, max_id)

Make a new unique identifier by incrementing on the maximum node id.
Hint: prefer using `max_id = max_id(mtg)` and then `new_id(mtg, max_is)` for performance
if you do it repeatidely.
"""
function new_id(max_id::Int)
    max_id + 1
end

function new_id(mtg::Node)
    new_id(max_id(mtg))
end


text:
tests = [
    "hdiFromIcdf",
]

print_with_color(:blue, "Running tests:\n")

srand(345678)

for t in tests
    test_fn = "$t.jl"
    print_with_color(:green, "* $test_fn\n")
    include(test_fn)
end


text:
include("../examples/rijke.jl")
include("../src/clvs.jl")
using PyPlot
using OrdinaryDiffEq
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function compute_acoustic_energy()
    n_samples = 112
    d_u = 3
    Eac = SharedArray{Float64}(n_samples)
    les = SharedArray{Float64}(d_u, n_samples)
    les .= 0.
    Eac .= 0.
    beta = LinRange(2.0,9.0,n_samples)
    nSteps = 5000
    nRunup = 400000
    answer = @distributed for i = 1:n_samples
    println("beta = ", beta[i])
    s = [beta[i], 0.2]
    u = Rijke_ODE(rand(N), s, nRunup)
    du_trj = zeros(N, N, nSteps)
    for j = 1:nSteps
        u = Rijke_ODE(u, s, 1)
        du_trj[:,:,j] = dRijke(u, s, 1.e-6)
        Eac[i] += norm(u[1:2*Ng])^2.0/4/nSteps
    end
    les[:,i], clv_trj = clvs(du_trj, d_u)
    @show Eac[i]
    end
    wait(answer)
    save("../data/attractor/Eac1.jld", "beta", beta,
     "Eac", Eac, "les", les)
    return beta, Eac, les

end
function plot_acoustic_energy_colored()
	X = load("../data/attractor/more_les/Eac.jld")
	Y = load("../data/attractor/more_les/Eac1.jld")
	Z = load("../data/attractor/Eac.jld")
	W = load("../data/attractor/Eac_and_beta.jld")
	Eac1 = X["Eac"]
	Eac2 = Y["Eac"]
	beta1 = X["beta"]
	beta2 = Y["beta"]
	les1 = X["les"]/dt
	les2 = Y["les"]/dt
	Eac3 = Z["Eac"]
	Eac4 = W["Eac"]
	beta3 = Z["beta"]
	beta4 = W["beta"]

    fig1, ax1 = subplots(1,1)
	ax1.plot(beta1, les1[1,:], ".-")
	ax1.plot(beta1, les1[2,:], ".-")
	ax1.plot(beta1, les1[3,:], ".-")
    ax1.xaxis.set_tick_params(labelsize=28)
    ax1.yaxis.set_tick_params(labelsize=28)
    ax1.set_xlabel(L"$\beta$",fontsize=28)
    ax1.set_ylabel(L"$\lambda$",fontsize=28)
	grid(true)

    fig2, ax2 = subplots(1,1)
	ax2.plot(beta2, les2[1,:], ".-")
	ax2.plot(beta2, les2[2,:], ".-")
	ax2.plot(beta2, les2[3,:], ".-")
    ax2.xaxis.set_tick_params(labelsize=28)
    ax2.yaxis.set_tick_params(labelsize=28)
    ax2.set_xlabel(L"$\beta$",fontsize=28)
    ax2.set_ylabel(L"$\lambda$",fontsize=28)
	grid(true)

	b_c_ll = 6.2
	b_c_ul = 7.3
	b_qp_ll_1 = 5.0
	b_qp_ul_1 = 6.2
	b_qp_ll_2 = 5.0
	b_qp_ul_2 = 6.2

	beta1_c = beta1[b_c_ll .< beta1 .< b_c_ul]
	chaotic_Eac1 = Eac1[b_c_ll .< beta1 .< b_c_ul]
	beta2_c = beta2[b_c_ll .< beta2 .< b_c_ul]
	chaotic_Eac2 = Eac2[b_c_ll .< beta2 .< b_c_ul]

	beta3_c = beta3[b_c_ll .< beta3 .< b_c_ul]
	chaotic_Eac3 = Eac3[b_c_ll .< beta3 .< b_c_ul]
	beta4_c = beta4[b_c_ll .< beta4 .< b_c_ul]
	chaotic_Eac4 = Eac4[b_c_ll .< beta4 .< b_c_ul]

	beta1_qp_1 = beta1[b_qp_ll_1 .< beta1 .< b_qp_ul_1]
	qp1_Eac1 = Eac1[b_qp_ll_1 .< beta1 .< b_qp_ul_1]
	beta2_qp_1 = beta2[b_qp_ll_1 .< beta2 .< b_qp_ul_1]
	qp1_Eac2 = Eac2[b_qp_ll_1 .< beta2 .< b_qp_ul_1]

	beta1_qp_2 = beta1[b_qp_ll_2 .< beta1 .< b_qp_ul_2]
	qp2_Eac1 = Eac1[b_qp_ll_2 .< beta1 .< b_qp_ul_2]
	beta2_qp_2 = beta2[b_qp_ll_2 .< beta2 .< b_qp_ul_2]
	qp2_Eac2 = Eac2[b_qp_ll_2 .< beta2 .< b_qp_ul_2]

	beta3_qp_1 = beta3[b_qp_ll_1 .< beta3 .< b_qp_ul_1]
	qp1_Eac3 = Eac3[b_qp_ll_1 .< beta3 .< b_qp_ul_1]
	beta4_qp_1 = beta4[b_qp_ll_1 .< beta4 .< b_qp_ul_1]
	qp1_Eac4 = Eac4[b_qp_ll_1 .< beta4 .< b_qp_ul_1]

	beta3_qp_2 = beta3[b_qp_ll_2 .< beta3 .< b_qp_ul_2]
	qp2_Eac3 = Eac3[b_qp_ll_2 .< beta3 .< b_qp_ul_2]
	beta4_qp_2 = beta4[b_qp_ll_2 .< beta4 .< b_qp_ul_2]
	qp2_Eac4 = Eac4[b_qp_ll_2 .< beta4 .< b_qp_ul_2]

	fig, ax = subplots(1,1)
    ax.plot(beta1, Eac1, "b.")
    ax.plot(beta2, Eac2, "b.")
    ax.plot(beta1_c, chaotic_Eac1, "r.")
    ax.plot(beta2_c, chaotic_Eac2, "r.")
	ax.plot(beta1_qp_1, qp1_Eac1, "g.")
	ax.plot(beta1_qp_2, qp2_Eac1, "g.")
	ax.plot(beta2_qp_1, qp1_Eac2, "g.")
	ax.plot(beta2_qp_2, qp2_Eac2, "g.")
	ax.plot(beta3, Eac3, "b.")
    ax.plot(beta4, Eac4, "b.")
    ax.plot(beta3_c, chaotic_Eac3, "r.")
    ax.plot(beta4_c, chaotic_Eac4, "r.")
	ax.plot(beta3_qp_1, qp1_Eac3, "g.")
	ax.plot(beta3_qp_2, qp2_Eac3, "g.")
	ax.plot(beta4_qp_1, qp1_Eac4, "g.")
	ax.plot(beta4_qp_2, qp2_Eac4, "g.")

	ax_in = fig.add_axes([0.2,0.75,0.15,0.1])
	ax_in.plot(beta1, Eac1, "b.")
    ax_in.plot(beta2, Eac2, "b.")
    ax_in.plot(beta1_c, chaotic_Eac1, "r.")
    ax_in.plot(beta2_c, chaotic_Eac2, "r.")
	ax_in.plot(beta1_qp_1, qp1_Eac1, "g.")
	ax_in.plot(beta1_qp_2, qp2_Eac1, "g.")
	ax_in.plot(beta2_qp_1, qp1_Eac2, "g.")
	ax_in.plot(beta2_qp_2, qp2_Eac2, "g.")
	ax_in.plot(beta3, Eac3, "b.")
    ax_in.plot(beta4, Eac4, "b.")
    ax_in.plot(beta3_c, chaotic_Eac3, "r.")
    ax_in.plot(beta4_c, chaotic_Eac4, "r.")
	ax_in.plot(beta3_qp_1, qp1_Eac3, "g.")
	ax_in.plot(beta3_qp_2, qp2_Eac3, "g.")
	ax_in.plot(beta4_qp_1, qp1_Eac4, "g.")
	ax_in.plot(beta4_qp_2, qp2_Eac4, "g.")
	ax_in.set_xlim([6.6,7.35])
	ax_in.set_ylim([17,28])

	ax_in.xaxis.set_tick_params(labelsize=18)
    ax_in.yaxis.set_tick_params(labelsize=18)
	ax_in.grid(true)
    ax.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
    ax.set_xlabel(L"$\beta$",fontsize=28)
    ax.set_ylabel(L"$\langle J_{\rm ac}\rangle$",rotation=0,
				  fontsize=28)
	ax.grid(true)
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies