text:
module LatexPrint
import DataFrames.DataFrame
export latex_form, laprint, laprintln, lap, tabular

export set_nan, set_inf, set_emptyset, set_delims, set_align
export set_im, set_bool

const EOL = " \\\\\n"  # end rows in matrices
const TAB = " & "      # tab for matrices

global NAN = "\\text{NaN}"     # default for NaN
global INF = "\\infty"         # default for infinity
global TRUE = "\\mathrm{T}"    # default for boolean true
global FALSE = "\\mathrm{F}"   # default for boolean false
global IM = "i"                # default for sqrt(-1)

global EMPTYSET = "\\emptyset" # default for empty set

global LEFT = "["    # default left delimiter for matrices
global RIGHT = "]"   # default right delimiter
global ALIGN = "c"   # default alignment character

# The set_xxx functions may take a string or a single character
Ctype = Union{String,Char}

"""
`set_nan(msg)` sets the `latex_print` output for a `NaN` value.
"""
function set_nan(msg::Ctype)
    global NAN = string(msg)
end

"""
`set_inf(msg)` sets the `latex_print` output for infinity.
"""
function set_inf(msg::Ctype)
    global INF = string(msg)
end

"""
`set_emptyset(msg)` sets the `latex_print` output for the empty set.
"""
function set_emptyset(msg::Ctype)
    global EMPTYSET = string(msg)
end

"""
`set_delims(lt,rt)` sets the `latex_print` output for left and right
matrix delimeters.
"""
function set_delims(lt::Ctype, rt::Ctype)
    global LEFT = string(lt)
    global RIGHT = string(rt)
    return (lt,rt)
end

"""
`set_align(ch)` sets the alignment character for matrices. Should
be one of `l`, `r`, or `c`.
"""
function set_align(ch::Ctype)
    global ALIGN = string(ch)
    if length(ALIGN) != 1
        warn("Alignment character should be a single character")
    end
end

"""
`set_im(ch)` sets the output for \$i\$ (square root of -1).
"""
function set_im(i::Ctype)
    global IM = string(i)
end

"""
`set_bool(tmsg, fmsg)` sets the output for the Boolean values
`true` and `false`.
"""
function set_bool(t::Ctype, f::Ctype)
    global TRUE = string(t)
    global FALSE = string(f)
    return (TRUE,FALSE)
end

# The latex_form function is the central workhorse for converting a
# Julia object to a character string (of type String) that can be
# printed (and then pasted into LaTeX mathmode). We create a version
# for each Julia datatype that we might want to render in LaTeX. Users
# can define more (as described in the README document).

# character strings we wrap in "\text"
"""
`latex_form(arg)` writes `arg` to the screen in a form
suitable for pasting into a LaTeX document.
"""
latex_form(words::AbstractString) = "\\text{" * words * "}"

# floating point numbers are passed through, unless NaN or Inf
function latex_form(x::AbstractFloat)
    if isnan(x)
        return NAN
    end
    if isinf(x)
        if x>0
            return INF
        end
        return "-" * INF
    end
    return string(x)
end

# special case for MathConst's
function latex_form(x::Irrational)
    if x==pi
        return "\\pi"
    end

    if x==e
        return "e"
    end

    if x==golden
        return "\\phi"
    end

    if x==eulergamma
        return "\\gamma"
    end

    return latex_form(float(x))  # some math constant we don't know
end

# integers are easy: just convert to string
latex_form(x::Integer) = string(x)

# booleans
latex_form(x::Bool) = x ? TRUE : FALSE

# rational numbers are presented as \frac fractions
function latex_form{T}(q::Rational{T})
    a::T = num(q)
    b::T = den(q)

    if b==0  # This is some form of infinity
        if a>0
            return INF
        else
            return "-" * INF
        end
    end

    if b==1  # Just print the numerator; it's an integer value
        return latex_form(a)
    end

    return "\\frac{" * latex_form(a) * "}{" * latex_form(b) * "}"
end

# Complex numbers
function latex_form{T}(z::Complex{T})
    if isnan(z)
        return NAN
    end
    if isinf(z)
        return INF
    end

    a,b = reim(z)

    if T == Bool  # Complex{Bool} treated like Complex{Int}
        a = Int(a)
        b = Int(b)
    end

    op = b >= 0 ? "+" : "-"

    return latex_form(a) * op * latex_form(abs(b)) * IM
end

# Sets
function latex_form(A::Union{Set,IntSet})
    if isempty(A)
        return EMPTYSET
    end
    elements = collect(A)
    try
        sort!(elements)
    end
    n = length(elements)

    result = "\\left\\{"
    for k=1:n
        result *= latex_form(elements[k])
        if k<n
            result *= ","
        else
            result *= "\\right\\}"
        end
    end
    return result
end

# Vectors (1-dimensional arrays)
function latex_form{T}(A::Array{T,1})
    result = "\\left" * LEFT * "\n\\begin{array}{" * ALIGN * "}\n"
    for x in A
        result *= latex_form(x)*EOL
    end
    result *= "\\end{array}\n\\right" * RIGHT
    return result
end

# Matrices (2-dimensional arrays)
function latex_form{T}(A::Array{T,2})
    (r,c) = size(A)

    # Header
    result  = "\\left" * LEFT * "\n\\begin{array}{"
    result *= (ALIGN ^ c) * "}\n"

    # Row-by-row
    for a=1:r
        for b=1:c
            result *= latex_form(A[a,b])
            if b<c
                result *= TAB
            else
                result *= EOL
            end
        end
    end

    # close up
    result *= "\\end{array}\n\\right" * RIGHT
    return result
end

# catch all for any types we've not implemented
latex_form(x::Any) = string(x)

# The print functions in all their forms. They simply apply latex_form
# to their argument(s) and pass the resulting string(s) to the
# corresponding version of print.
"""
`laprint(x...)` is like Julia's `print` but converts each
argument using `latex_print`.
"""
 function laprint(x...)
    xs = map(latex_form,x)
    print(xs...)
end

"""
`laprintln(x...)` is like Julia's `print` but converts each
argument using `latex_print`.
"""
function laprintln(x...)
    xs = map(latex_form,x)
    println(xs...)
end

"""
`lap(x...)` is a short version of `laprintln(x...)`.
"""
lap(x...) = laprintln(x...)  # abbreviation

# Same functions, but now with IO option.

function laprint(io::IO, x...)
    xs = map(latex_form,x)
    print(io, xs...)
end

function laprintln(io::IO, x...)
    xs = map(latex_form,x)
    println(io, xs...)
end

lap(io::IO, x...) = laprintln(io, x...)

function tabular{T}(A::Array{T,2}, alignment::String)
    (r,c) = size(A)
    println("\\begin{tabular}{", alignment, "}")
    for a=1:r
        for b=1:c
            print("\$",latex_form(A[a,b]),"\$")
            if b<c
                print(" & ")
            else
                if a<r
                    println("\\\\")
                end
            end
        end
    end
    println("\n\\end{tabular}")
end

function tabular{T}(A::Array{T,2})
    (r,c) = size(A)
    alignment = ALIGN^c
    tabular(A,alignment)
end

function tabular(A::DataFrame, alignment::String; rounding::Int = 0)
    (r,c) = size(A)
    println("\\begin{tabular}{", alignment, "}")
    println("\\hline")
    for (i,name) in enumerate(names(A))
        print(name)
        if i < length(names(A))
            print(" & ")
        end
    end
    println("\\\\")
    println("\\hline")
    for a=1:r
        for b=1:c
            print("\$",latex_form((rounding > 0) && typeof(A[a,b])<:Number ? round(A[a,b],rounding) : A[a,b]),"\$")
            if b<c
                print(" & ")
            else
                if a<r
                    println("\\\\")
                end
            end
        end
    end
    println("\n\\end{tabular}")
end

function tabular(A::DataFrame; rounding::Int = 0)
    (r,c) = size(A)
    alignment = ALIGN^c
    tabular(A,alignment; rounding = rounding)
end

end # end of module LatexPrint


text:
using Test
using GeoParams
@testset "MeltingParam.jl" begin

#Make sure structure is isbits
x = MeltingParam_Caricchi()
@test isbits(x)

# This tests the various melting parameterizations
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=10km);

T       =   Vector(250:100:1250)*K .+ 273.15K
T_nd    =   Float64.(T/CharUnits_GEO.Temperature)

# Caricchi parameterization [in ND numbers, which is anyways the typical use case]
p        =  MeltingParam_Caricchi()
phi_dim  =  zeros(size(T))
compute_meltfraction!(phi_dim, p, zeros(size(T)), ustrip.(T))

phi_dim1  = zeros(size(phi_dim))
compute_meltfraction!(phi_dim1, p, zeros(size(T)), ustrip.(T)) # in-place routine

p_nd     =  p
p_nd     =  nondimensionalize(p_nd, CharUnits_GEO)
phi_nd   = zeros(size(T))
compute_meltfraction!(phi_nd, p_nd, zeros(size(T_nd)),T_nd)

# Do this computation manually, using the actual expression of Caricchi
T_C         =   Vector(250:100:1250)    # in celcius
Phi_solid   =   1.0 .- 1.0./(1.0 .+ exp.((800.0 .- T_C)./23.0));
Phi_anal    =   1.0 .- Phi_solid

@test sum(phi_dim  - Phi_anal)   < 1e-12
@test sum(phi_dim1 - Phi_anal)   < 1e-12
@test sum(phi_nd - Phi_anal)    < 1e-12

# Test computation of density for the whole computational domain, using arrays
MatParam    =   Array{MaterialParams, 1}(undef, 3);
MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                        Melting  = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"));

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                        Melting   = MeltingParam_Caricchi());

# No melting parameterization for this phase
MatParam[3] =   SetMaterialParams(Name="UpperCrust", Phase=3,
                        Density   = PT_Density());

Mat_tup = Tuple(MatParam)

# test computing material properties
n = 100;
Phases              = ones(Int64,n,n,n);
Phases[:,:,20:end] .= 2
Phases[:,:,80:end] .= 3

ϕ = zeros(size(Phases))
T =  ones(size(Phases))*1500
P =  ones(size(Phases))*10

compute_meltfraction!(ϕ, Mat_tup, Phases, P,T) #allocations coming from computing meltfraction using PhaseDiagram_LookupTable

@test sum(ϕ)/n^3 ≈ 0.6463001302812086

# test computing material properties when we have PhaseRatios, instead of Phase numbers
PhaseRatio  = zeros(n,n,n,3);
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz)
    PhaseRatio[I] = 1.0
end

compute_meltfraction!(ϕ, Mat_tup, PhaseRatio, P,T)
@test sum(ϕ)/n^3 ≈ 0.6463001302812086

end

ext: jl


text:
@testset "successor method for discrete systems" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    c = [0.1; 0.3]
    D = [1. 1 1; 1 -1 1]
    X = Hyperrectangle([0.0, 0.3], [2.0, 2.0])
    U = Interval(0.0, 1.0)
    W = Hyperrectangle([0.0, 0.3, 0.5], [2.0, 2.0, 2.0])
    f(x) = x'*x
    f(x,u) = x'*x + u'*u
    f(x,u,w) = x'*x + u'*u + w'*w
    x = [0.0; 0.3]
    u = [0.5]
    w = [0.5, 1.0, 1.3]

    sys = DiscreteIdentitySystem(size(A, 2))
    @test successor(sys, x) == x
    sys = ConstrainedDiscreteIdentitySystem(size(A, 2), X)
    @test successor(sys, x) == x
    sys = LinearDiscreteSystem(A)
    @test successor(sys, x) == A*x
    sys = AffineDiscreteSystem(A, c)
    @test successor(sys, x) == A*x + c
    sys = LinearControlDiscreteSystem(A, B)
    @test successor(sys, x, u) == A*x + B*u
    sys = ConstrainedLinearDiscreteSystem(A, X)
    @test successor(sys, x) == A*x
    sys = ConstrainedAffineDiscreteSystem(A, c, X)
    @test successor(sys, x) == A*x + c
    sys = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    @test successor(sys, x, u) == A*x + B*u
    sys = ConstrainedAffineControlDiscreteSystem(A, B, c, X, U)
    @test successor(sys, x, u) == A*x + B*u + c
    sys = BlackBoxDiscreteSystem(f, size(A, 2))
    @test successor(sys, x) == f(x)
    sys = ConstrainedBlackBoxDiscreteSystem(f, size(A, 2), X)
    @test successor(sys, x) == f(x)
    sys = ConstrainedBlackBoxControlDiscreteSystem(f, size(A, 2), size(B, 2), X, U)
    @test successor(sys, x, u) == f(x, u)
    sys = NoisyConstrainedLinearDiscreteSystem(A, D, X, W)
    @test successor(sys, x, w) == A*x + D*w
    sys = NoisyConstrainedLinearControlDiscreteSystem(A, B, D, X, U, W)
    @test successor(sys, x, u, w) == A*x + B*u + D*w
    sys = NoisyConstrainedAffineControlDiscreteSystem(A, B, c, D, X, U, W)
    @test successor(sys, x, u, w) == A*x + B*u + c + D*w
    sys = NoisyConstrainedBlackBoxControlDiscreteSystem(f, size(A, 2), size(B, 2), size(D, 2), X, U, W)
    @test successor(sys, x, u, w) == f(x, u, w)
end

@testset "successor exception handling" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    c = [0.1; 0.3]
    D = [1. 1 1; 1 -1 1]
    X = Hyperrectangle([0.0, 0.3], [2.0, 2.0])
    U = Interval(0.0, 1.0)
    W = Hyperrectangle([0.0, 0.3, 0.5], [2.0, 2.0, 2.0])
    f(x) = x'*x
    f(x,u) = x'*x + u'*u
    f(x,u,w) = x'*x + u'*u + w'*w
    x = [0.0; 0.3]
    x_wrong_dim = [0.0]
    x_not_in_set = [-2.1; 0.3]
    u = [0.5]
    u_wrong_dim = [0.5, 0.3]
    u_not_in_set = [-1.0]
    w = [0.5, 1.0, 1.3]
    w_wrong_dim = [0.5, 1.0, 1.3, 1.0]
    w_not_in_set = [-2.1, 1.0, 1.3]

    sys = DiscreteIdentitySystem(size(A,2))
    @test_throws ArgumentError("the state vector has the wrong dimensions") successor(sys, x_wrong_dim)
    sys = ConstrainedDiscreteIdentitySystem(size(A, 2), X)
    @test_throws ArgumentError("the state vector is not contained in the state set")  successor(sys, x_not_in_set)
    sys = LinearDiscreteSystem(A)
    @test_throws ArgumentError("the state vector has the wrong dimensions") successor(sys, x_wrong_dim)
    sys = AffineDiscreteSystem(A, c)
    @test_throws ArgumentError("the state vector has the wrong dimensions") successor(sys, x_wrong_dim)
    sys = LinearControlDiscreteSystem(A, B)
    @test_throws ArgumentError("the input vector has the wrong dimensions")  successor(sys, x, u_wrong_dim)
    sys = ConstrainedLinearDiscreteSystem(A, X)
    @test_throws ArgumentError("the state vector is not contained in the state set")  successor(sys, x_not_in_set)
    sys = ConstrainedAffineDiscreteSystem(A, c, X)
    @test_throws ArgumentError("the state vector is not contained in the state set")  successor(sys, x_not_in_set)
    sys = ConstrainedLinearControlDiscreteSystem(A, B, X, U)
    @test_throws ArgumentError("the input vector is not contained in the input set")  successor(sys, x, u_not_in_set)
    sys = ConstrainedAffineControlDiscreteSystem(A, B, c, X, U)
    @test_throws ArgumentError("the input vector has the wrong dimensions")  successor(sys, x, u_wrong_dim)
    sys = BlackBoxDiscreteSystem(f, size(A, 2))
    @test_throws ArgumentError("the state vector has the wrong dimensions")  successor(sys, x_wrong_dim)
    sys = ConstrainedBlackBoxDiscreteSystem(f, size(A, 2), X)
    @test_throws ArgumentError("the state vector is not contained in the state set")  successor(sys, x_not_in_set)
    sys = ConstrainedBlackBoxControlDiscreteSystem(f, size(A, 2), size(B, 2), X, U)
    @test_throws ArgumentError("the input vector is not contained in the input set")  successor(sys, x, u_not_in_set)
    sys = NoisyConstrainedLinearDiscreteSystem(A, D, X, W)
    @test_throws ArgumentError("the noise vector is not contained in the noise set")  successor(sys, x, w_not_in_set)
    sys = NoisyConstrainedLinearControlDiscreteSystem(A, B, D, X, U, W)
    @test_throws ArgumentError("the noise vector has the wrong dimensions")  successor(sys, x, u, w_wrong_dim)
    sys = NoisyConstrainedAffineControlDiscreteSystem(A, B, c, D, X, U, W)
    @test_throws ArgumentError("the noise vector is not contained in the noise set")  successor(sys, x, u, w_not_in_set)
    sys = NoisyConstrainedBlackBoxControlDiscreteSystem(f, size(A, 2), size(B, 2), size(D, 2), X, U, W)
    @test_throws ArgumentError("the noise vector has the wrong dimensions")  successor(sys, x, u, w_wrong_dim)
end


text:
# Some passes for assigning node affinity.
#
# This applies heuristics to nodes like `Broadcast` and `Result` to schedule them in more
# sensible locations

# A follower node is one that should be scheduled as soon as possible
const FOLLOWERS = [
    "Result",
    "Sum",
    "Add",
    "Subtract",
    "MaxPool",
    "AvgPool",
    #"Concat",
    "Max",
    "Min",
    "Multiply",
    "Relu",
    "Reshape",
    "Dot",
]
function is_follower(node)::Bool
    # If this is a move node into remote memory
    if (ismovesync(node) && nGraph.is_persistent(first(outputs(node)))) || ismoveasync(node)
        return true
    end

    # Default set of follower nodes
    description = nGraph.description(node)
    return any(x -> startswith(description, x), FOLLOWERS)
end

# A leader node should be scheduled as late as possible, but clumped before its trailing
# output
const LEADERS = [
    "Broadcast",
    "ConvertLayout",
    "Slice",
]
function is_leader(node)::Bool
    # If this is a move node prefetching memory, do this as soon as possible.
    if (ismovesync(node) && !nGraph.is_persistent(first(outputs(node))))
        return true
    end

    # Default leader nodes
    description = nGraph.description(node)
    return any(x -> startswith(description, x), LEADERS)
end

# Heuristic to assign priorities to nodes in the graph to yield better schedules
function priority_pass!(f::nGraph.NFunction)
    # First, construct a map from nodes to their index in the function
    node_to_index = Dict(NodeDescriptor(n) => i for (i, n) in enumerate(f))

    # The final node priorities to assign. Save this until the end do avoid accidentally
    # mutating the underlying function too early.
    #
    # We avoid assigning priorities to normal nodes, which will just default internally
    # to a "zero" priority
    priorities = Dict{NodeDescriptor, Int64}()

    # Next, assign priorities
    for node in map(NodeDescriptor, f)
        # Async Moves take priority because if they aren't placed right after their async
        # node - codegen will clump them with the wrong computation kernel.
        if ismoveasync(node)
            priorities[node] = -2

        elseif is_follower(node)
            priorities[node] = -1

        elseif is_leader(node)
            # Find the minimum index of all outputs and control dependencies of the node
            deps = nGraph.get_outputs(node)
            p = minimum(get(node_to_index, n, typemax(Int64)) for n in deps)
            priorities[node] = p
        end
    end

    for (n,p) in priorities
        nGraph.set_priority(n, p)
    end
    nGraph.get_ordered_ops!(f)
    return nothing
end

#####
##### Scratchpad
#####

# Optimization for 2LM: assign short lived tensors to a different memory pool so they all
# end up at similar memory addresses
function setup_scratchpad!(data::Profiler.FunctionData; threshold = 100)
    for tensor in tensors(data)
        start = producer(tensor).index
        stop = consumer(tensor).index

        if stop - start >= threshold
            nGraph.set_pool_number(tensor.tensor, 2)
        end
    end
    return nothing
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Test
using Mimi

@testset "Standard API" begin

# Test that the function does not error and returns a valid value
scc1 = compute_scc(year=2020)
@test scc1 isa Float64

# Test that a higher discount rate makes a lower scc value
scc2 = compute_scc(year=2020, eta=0., prtp=0.03)
@test scc2 < scc1

# Test with a modified model
m = get_model()
update_param!(m, :tcr_transientresponse, 3)
scc3 = compute_scc(m, year=2020)
@test scc3 > scc1

# Test get_marginal_model
mm = get_marginal_model(year = 2040)
mm[:ClimateTemperature, :rt_g_globaltemperature]

# Test compute_scc_mm
result = compute_scc_mm(year=2050)
@test result.scc > scc1
@test result.mm isa Mimi.MarginalModel

end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using Documenter, OAuth

makedocs(
    sitename = "OAuth.jl"
)

deploydocs(
    deps = Deps.pip("mkdocs==0.16.3", "python-markdown-math", "mkdocs-bootswatch"),
    repo = "github.com/randyzwitch/OAuth.jl.git",
    julia  = "1.3",
    osname = "linux"
)


text:
#=
acceleration.jl

    Acceleration routines for proximal optimization algorithms, such as simple
    extrapolation steps, Nesterov momentum steps, etc.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/01/18
=#

"""
    update_acc!(acc, state)

Update the extrapolation parameter `ω` using a specific acceleration scheme,
such as simple extrapolation and Nesterov momentum.

#### Arguments
  - `acc::AbstractAccelScheme`	: acceleration scheme variables
  - `state::ProxGradState`		: state variables
"""
function update_acc! end

function update_acc!(acc::NoAccel, state::ProxGradState)
	# Do nothing...

	return nothing
end

function update_acc!(acc::Simple, state::ProxGradState)
	# Gradient based adaptive restart
    acc.k= dot(state.∇f, state.Δ) > zero(eltype(state.Δ)) ? 1 : acc.k

    # Simple extrapolation
    acc.ω= (acc.k - 1)*inv(acc.k + 2)

	return nothing
end

function update_acc!(acc::Nesterov, state::ProxGradState)
	# Gradient based adaptive restart
    acc.θ= dot(state.∇f, state.Δ) > zero(eltype(state.Δ)) ? one(acc.θ) : acc.θ

	# Clean up notation
	λ_prev= acc.ls.λ_prev
	λ= acc.ls.λ
	θ_prev= acc.θ

    # Quadratic formula
    b = θ_prev^2 * inv(λ_prev) - acc.m
    D= b^2 + 4 * inv(λ_prev * λ) * θ_prev^2

    # Update Nesterov momentum
    acc.θ= .5 * λ * (-b + sqrt(D))

    # Nesterov extrapolation
    acc.ω= λ * θ_prev * (one(θ_prev) - θ_prev) * inv(λ_prev * acc.θ + λ * θ_prev^2)

	return nothing
end


text:
@testset "ProcessEstimation.jl" begin
    df = CSV.read("Daily__Oct-12-2019_03_30_12PM.csv", DataFrame)
    dropmissing!(df)
    df.Value = - df.Value
    df.Date = (s -> Date(s, "yyyy/mm/dd")).(df.Date)
    df.Exceedances = (v -> max(0, v + 13.66)).(df.Value)

    times = df[df.Exceedances .> 0, :Date]
    marks = df[df.Exceedances .> 0, :Exceedances]

    ts = TimeSeries(times)
    mts = MarkedTimeSeries(times, marks)

    μ = 2.0
    ϕ = 1
    γ = 1
    δ = 1.0

    ξ = 1.5
    β = 2
    α = 1.1
    κ = 1.5

    @testset "negloglik(ts ; μ, ϕ, γ)" begin
        @test_logs (:warn, "ϕ must be positive or zero, taking absolute value") SEMPP.negloglik(ts, μ = μ, ϕ = -ϕ, γ = γ)
        @test SEMPP.negloglik(ts, μ = μ, ϕ = ϕ, γ = γ) isa Real
    end

    sepp = SEMPP.SEPPExpKern(ts, μ = μ, ϕ = ϕ, γ = γ)

    @testset "negloglik(sepp)" begin
        @test SEMPP.negloglik(sepp) == SEMPP.negloglik(ts, μ = μ, ϕ = ϕ, γ = γ)
    end

    @testset "negloglik(mts, mardens ; μ, ϕ, γ, δ, ξ, α, β, κ)" begin
        @test_logs (:warn, "μ γ must be positive or zero, taking absolute value") SEMPP.negloglik(mts, markdens = Distributions.GeneralizedPareto, μ = -μ, ϕ = ϕ, γ = -γ, ξ = ξ, α = α, β = β)
        @test SEMPP.negloglik(mts, markdens = Distributions.GeneralizedPareto, μ = -μ, ϕ = ϕ, γ = -γ, ξ = ξ, α = α, β = β) isa Real
    end

    sempp_egpd = SEMPPExpKern(mts, μ = μ, ϕ = ϕ, γ = γ, δ = δ, markdens = EGPD.EGPpower, ξ = ξ, α = α, β = β, κ = κ)
    sempp_gpd = SEMPPExpKern(mts, μ = μ, ϕ = ϕ, γ = γ, δ = δ, markdens = Distributions.GeneralizedPareto, ξ = ξ, α = α, β = β)

    @testset "negloglik(sempp)" begin
        @test SEMPP.negloglik(sempp_egpd) == SEMPP.negloglik(mts, markdens = EGPD.EGPpower, μ = μ, ϕ = ϕ, γ = γ, δ = δ, ξ = ξ, α = α, β = β, κ = κ)
    end

    @testset "fit!(sepp)" begin
        @test fit!(sepp) isa Real
        @test μ != sepp.μ
        @test ϕ != sepp.ϕ
        @test γ != sepp.γ
        @test sepp.μ >= 0
        @test sepp.ϕ >= 0
        @test sepp.γ >= 0
    end

    @testset "fit!(sempp)" begin
        @test fit!(sempp_egpd) isa Real
        @test fit!(sempp_gpd) isa Real
        @test μ != sempp_egpd.μ
        @test ϕ != sempp_egpd.ϕ
        @test γ != sempp_egpd.γ
        @test δ != sempp_egpd.δ
        @test ξ != sempp_egpd.ξ
        @test α != sempp_egpd.α
        @test β != sempp_egpd.β
        @test κ != sempp_egpd.κ

        @test sempp_egpd.μ >=0
        @test sempp_egpd.ϕ >=0
        @test sempp_egpd.γ >=0
        @test sempp_egpd.δ >=0
        @test sempp_egpd.ξ >=0
        @test sempp_egpd.α >=0
        @test sempp_egpd.β >=0
        @test sempp_egpd.κ >=0

        @test μ != sempp_gpd.μ
        @test ϕ != sempp_gpd.ϕ
        @test γ != sempp_gpd.γ
        @test δ != sempp_gpd.δ
        @test ξ != sempp_gpd.ξ
        @test α != sempp_gpd.α
        @test β != sempp_gpd.β

        @test sempp_gpd.μ >=0
        @test sempp_gpd.ϕ >=0
        @test sempp_gpd.γ >=0
        @test sempp_gpd.δ >=0
        @test sempp_gpd.ξ >=0
        @test sempp_gpd.α >=0
        @test sempp_gpd.β >=0
    end

    @testset "Simulation.jl" begin

        @testset "PointProcess" begin
            simulated_ts = simulation(sepp, 100*365)

            @test simulated_ts isa TimeSeries

            sepp_simul = SEPPExpKern(simulated_ts)
            fit!(sepp_simul)

            @test abs(sepp_simul.μ - sepp.μ) < 1.96 * sepp_simul.cov_mat[1,1]
            @test abs(sepp_simul.ϕ - sepp.ϕ) < 1.96 * sepp_simul.cov_mat[2,2]
            @test abs(sepp_simul.γ - sepp.γ) < 1.96 * sepp_simul.cov_mat[3,3]
        end

        @testset "MarkedPointProcess GPD" begin
            simulated_mts = simulation(sempp_gpd, 100*365)

            @test simulated_mts isa MarkedTimeSeries

            sempp_gpd_simul = SEMPPExpKern(simulated_mts)
            fit!(sempp_gpd_simul)

            @test abs(sempp_gpd_simul.μ - sempp_gpd.μ) < 1.96 * sempp_gpd_simul.cov_mat[1,1]
            @test abs(sempp_gpd_simul.ϕ - smpp_gpd.ϕ) < 1.96 * sempp_gpd_simul.cov_mat[2,2]
            @test abs(sempp_gpd_simul.γ - sempp_gpd.γ) < 1.96 * sempp_gpd_simul.cov_mat[3,3]
        end

        @testset "MarkedPointProcess EGPD" begin
            simulated_mts = simulation(sempp_egpd, 100*365)

            @test simulated_mts isa MarkedTimeSeries

            sempp_egpd_simul = SEMPPExpKern(simulated_mts)
            fit!(sempp_egpd_simul)

            @test abs(sempp_egpd_simul.μ - sempp_egpd.μ) < 1.96 * sempp_egpd_simul.cov_mat[1,1]
            @test abs(sempp_egpd_simul.ϕ - smpp_egpd.ϕ) < 1.96 * sempp_egpd_simul.cov_mat[2,2]
            @test abs(sempp_egpd_simul.γ - sempp_egpd.γ) < 1.96 * sempp_egpd_simul.cov_mat[3,3]
        end
    end
end


text:
using Coverage

# coverage = process_folder(ARGS[1])
coverage = process_folder()

print("Overall coverage: ")
println(get_summary(coverage))

covered_lines, total_lines = get_summary(coverage)

for c in coverage
  print(c.filename * ": ")
  println(get_summary(c))
end


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.MVLinearGaussian
using StaticArrays

include("test.jl")

setRNGs(0)

model, theta, ys, ko = MVLinearGaussian.defaultMVLGModel(2, 10)
println(ko.logZhats)

## just run the algorithm a few times

testSMC(model, 1024*1024, 2, false)
testSMC(model, 1024*1024, 2, false, 0.5)
testSMCParallel(model, 1024*1024, 2, false)
testSMCParallel(model, 1024*1024, 2, false, 0.5)


text:
const Cftok = UInt32

const IPC_RMID = 0 # remove identifier
const IPC_SET  = 1 # set options
const IPC_STAT = 2 # get options

const IPC_CREAT  = 0o1000 # create entry if key does not exist
const IPC_EXCL   = 0o2000 # fail if key exists
const IPC_NOWAIT = 0o4000 # error if request must wait

const SEM_UNDO = 0x1000 # undo the operation on exit

@static if Sys.isapple()
    const GETPID   = 4 # get sempid
    const GETVAL   = 5 # get semval
    const GETALL   = 6 # get all semval's
    const GETNCNT  = 3 # get semncnt
    const GETZCNT  = 7 # get semzcnt
    const SETVAL   = 8 # set semval
    const SETALL   = 9 # set all semval's
else
    const GETPID   = 11 # get sempid  # 4
    const GETVAL   = 12 # get semval  # 5
    const GETALL   = 13 # get all semval's  # 6
    const GETNCNT  = 14 # get semncnt # 3
    const GETZCNT  = 15 # get semzcnt # 7
    const SETVAL   = 16 # set semval # 8
    const SETALL   = 17 # set all semval's # 9
end

const SEM_STAT = 18 # linux specific
const SEM_INFO = 19 # linux specific

struct SemBuf
    sem_num::Cushort  # semaphore index in array
    sem_op::Cshort    # semaphore operation
    sem_flg::Cshort   # operation flags
end

function SemBuf(sem_num, sem_op; wait::Bool=true, undo::Bool=false)
    flg = 0
    wait || (flg |= IPC_NOWAIT)
    undo && (flg |= SEM_UNDO)
    SemBuf(sem_num, sem_op, flg)
end

ftok(path::String, id::Int) = ccall(:ftok, Cftok, (Cstring,Cint), path, id)

function semcreate(tok::Cftok, nsems::Int; create::Bool=true, create_exclusive::Bool=false, permission::Integer=0o660)
    flags = Cint(permission)
    if create
        flags |= IPC_CREAT
        if create_exclusive
            flags |= IPC_EXCL
        end
    end
    id = ccall(:semget, Cint, (Cftok, Cint, Cint), tok, nsems, flags)
    systemerror("error creating semaphore", id < 0)
    id
end

function semrm(id::Cint)
    ret = ccall(:semctl, Cint, (Cint,Cint,Cint), id, Cint(0), IPC_RMID)
    systemerror("error deleting semaphore", ret < 0)
end

function semset(id::Cint, val::Cint, which=0)
    ret = ccall(:semctl, Cint, (Cint,Cint,Cint,Cint), id, which, SETVAL, val);
    systemerror("error setting semaphore", ret < 0)
end

function semset(id::Cint, vals::Vector{Cushort})
    ret = ccall(:semctl, Cint, (Cint,Cint,Cint,Ptr{Cushort}), id, 0, SETALL, vals);
    systemerror("error setting semaphore", ret < 0)
end

function semget(id::Cint, which=0)
    ccall(:semctl, Cint, (Cint,Cint,Cint), id, which, GETVAL);
end

function semget(id::Cint, vals::Vector{Cushort})
    ret = ccall(:semctl, Cint, (Cint,Cint,Cint,Ptr{Cushort}), id, 0, GETALL, vals);
    systemerror("error getting semaphore", ret < 0)
end

function semop(id::Cint, ops::Vector{SemBuf})
    ret = ccall(:semop, Cint, (Cint,Ptr{Nothing},Csize_t), id, ops, length(ops))
    systemerror("error in semaphore operation", ret < 0)
end


text:
struct ICKey{T <: InitialConditionType, U <: PSY.Component} <: OptimizationContainerKey
    meta::String
end

function ICKey(
    ::Type{T},
    ::Type{U},
    meta=CONTAINER_KEY_EMPTY_META,
) where {T <: InitialConditionType, U <: PSY.Component}
    if isabstracttype(U)
        error("Type $U can't be abstract")
    end
    return ICKey{T, U}(meta)
end

get_entry_type(::ICKey{T, U}) where {T <: InitialConditionType, U <: PSY.Component} = T
get_component_type(::ICKey{T, U}) where {T <: InitialConditionType, U <: PSY.Component} = U

"""
Container for the initial condition data
"""
struct InitialCondition{T <: InitialConditionType, U <: Union{PJ.ParameterRef, Float64}}
    component::PSY.Component
    value::U
end

function InitialCondition(
    ::Type{T},
    component::PSY.Component,
    value::U,
) where {T <: InitialConditionType, U <: Union{PJ.ParameterRef, Float64}}
    return InitialCondition{T, U}(component, value)
end

function InitialCondition(
    ::ICKey{T, U},
    component::U,
    value::V,
) where {
    T <: InitialConditionType,
    U <: PSY.Component,
    V <: Union{PJ.ParameterRef, Float64},
}
    return InitialCondition{T, U}(component, value)
end

function get_condition(p::InitialCondition{T, Float64}) where {T <: InitialConditionType}
    return p.value
end

function get_condition(
    p::InitialCondition{T, PJ.ParameterRef},
) where {T <: InitialConditionType}
    return jump_value(p.value)
end

get_component(ic::InitialCondition) = ic.component
get_value(ic::InitialCondition) = ic.value
get_component_name(ic::InitialCondition) = PSY.get_name(ic.component)
get_component_type(ic::InitialCondition) = typeof(ic.component)
get_ic_type(
    ::Type{InitialCondition{T, U}},
) where {T <: InitialConditionType, U <: Union{PJ.ParameterRef, Float64}} = T
get_ic_type(
    ::InitialCondition{T, U},
) where {T <: InitialConditionType, U <: Union{PJ.ParameterRef, Float64}} = T

"""
Stores data to populate initial conditions before the build call
"""
mutable struct InitialConditionsData
    duals::Dict{ConstraintKey, DataFrames.DataFrame}
    parameters::Dict{ParameterKey, DataFrames.DataFrame}
    variables::Dict{VariableKey, DataFrames.DataFrame}
    aux_variables::Dict{AuxVarKey, DataFrames.DataFrame}
end

function InitialConditionsData()
    return InitialConditionsData(
        Dict{ConstraintKey, DataFrames.DataFrame}(),
        Dict{ParameterKey, DataFrames.DataFrame}(),
        Dict{VariableKey, DataFrames.DataFrame}(),
        Dict{AuxVarKey, DataFrames.DataFrame}(),
    )
end

function get_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: VariableType, U <: Union{PSY.Component, PSY.System}}
    return ic_data.variables[VariableKey(T, U)]
end

function get_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: AuxVariableType, U <: Union{PSY.Component, PSY.System}}
    return ic_data.aux_variables[AuxVarKey(T, U)]
end

function get_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    return ic_data.duals[ConstraintKey(T, U)]
end

function get_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: ParameterType, U <: Union{PSY.Component, PSY.System}}
    return ic_data.parameters[ParameterKey(T, U)]
end

function has_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: VariableType, U <: Union{PSY.Component, PSY.System}}
    return haskey(ic_data.variables, VariableKey(T, U))
end

function has_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: AuxVariableType, U <: Union{PSY.Component, PSY.System}}
    return haskey(ic_data.aux_variables, AuxVarKey(T, U))
end

function has_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    return haskey(ic_data.duals, ConstraintKey(T, U))
end

function has_initial_condition_value(
    ic_data::InitialConditionsData,
    ::T,
    ::Type{U},
) where {T <: ParameterType, U <: Union{PSY.Component, PSY.System}}
    return haskey(ic_data.parameters, ParameterKey(T, U))
end

######################### Initial Conditions Definitions#####################################
struct DevicePower <: InitialConditionType end
struct DeviceAboveMinPower <: InitialConditionType end
struct DeviceStatus <: InitialConditionType end
struct InitialTimeDurationOn <: InitialConditionType end
struct InitialTimeDurationOff <: InitialConditionType end
struct InitialEnergyLevel <: InitialConditionType end
struct InitialEnergyLevelUp <: InitialConditionType end
struct InitialEnergyLevelDown <: InitialConditionType end
struct AreaControlError <: InitialConditionType end


text:
using IterableTables
using StatsModels
using GLM
using DataValues
using Test

@testset "StatsModels" begin

source_array = [(a=DataValue(1),b=DataValue(1.),c=DataValue("A")), (a=DataValue(2),b=DataValue(2.),c=DataValue("B")), (a=DataValue(3),b=DataValue(3.),c=DataValue("C"))]

# TODO add some test beyond just creating a ModelFrame
mf_array = StatsModels.ModelFrame(StatsModels.@formula(a~b), source_array)

# lm(StatsModels.@formula(a~b), source_df)

end


text:
import StaticArrays as SA
import MiniFB as MFB

const T = Float32

# img

const height_img = 512
const width_img = 256

const img = zeros(UInt32, height_img, width_img)
const fb = zeros(UInt32, width_img, height_img)

# colors

const black = MFB.mfb_rgb(0, 0, 0)
const white = MFB.mfb_rgb(255, 255, 255)
const gray = MFB.mfb_rgb(127, 127, 127)
const green = MFB.mfb_rgb(0, 255, 0)
const red = MFB.mfb_rgb(255, 0, 0)

# main

count = 128

function test_mfb()
    window = MFB.mfb_open("Test", width_img, height_img)
    MFB.mfb_set_keyboard_callback(window, keyboard_callback)

    while MFB.mfb_wait_sync(window)
        img[1:128, 1:64] .= green

        state = MFB.mfb_update(window, permutedims!(fb, img, (2, 1)))

        if state != MFB.STATE_OK
            break;
        end
    end

    MFB.mfb_close(window)
end

function keyboard_callback(window, key, mod, isPressed)::Cvoid
    if isPressed
        if key == MFB.KB_KEY_UP
            display(key)
            global count = count + 1
            img[count, 1:64] .= red
        elseif key == MFB.KB_KEY_DOWN
            display(key)
        elseif key == MFB.KB_KEY_LEFT
            display(key)
        elseif key == MFB.KB_KEY_RIGHT
            display(key)
        elseif key == MFB.KB_KEY_ESCAPE
            display(key)
            MFB.mfb_close(window)
        end
    end

    return nothing
end

test_mfb()


text:
using GeometricIntegrators
using Weave

ENV["GKSwstype"] = "100"

tutorial_path = joinpath("tutorial", "tutorial.jmd")
build_path = joinpath("src", "tutorial")

weave(tutorial_path,
  out_path=build_path,
  doctype = "github")


text:
include("colorscheme.jl")
using PyPlot
using3D()

DrWatson._wsave(s, fig::Figure) = fig.savefig(s, dpi = 600, transparent = false)

PyPlot.rc("lines", lw = 2.6)
PyPlot.rc("errorbar", capsize = 6)
PyPlot.rc("axes", grid = true)
PyPlot.rc("grid", color = "0.75", alpha = 0.75)

PyPlot.rc("font", size = 26) # set default fontsize
PyPlot.rc("xtick", labelsize = 24)
PyPlot.rc("ytick", labelsize = 24)
PyPlot.rc("axes", labelsize = 32)
PyPlot.rc("legend", fontsize = 26)
# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

for z in ("x", "y")
    PyPlot.rc("$(z)tick.major", size = 7, width = 1.2)
    PyPlot.rc("$(z)tick.minor", size = 3, visible = false)
end

figx = 16 # default width correspoding to full text width
figy = 5  # default height corresponding to 1 row of plots
PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("savefig", dpi = 600, transparent = true, format = "png")

# set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORS.c))

if false # test font
    figure(figsize=(10,8)); plot(rand(10), label = "\$a=\\int_0^\\infty xdx\$")
    xlabel("x label"); ylabel("\$H_2\$"); legend(); tight_layout()
end

if false # test color scheme
    figure(figsize = (20, 10)) # show colors
    ax1 = subplot(131)
    ax2 = subplot(132)
    ax3 = subplot(133)
    lw = 120
    for (i, c) in enumerate(COLORS)
        chsv = matplotlib.colors.rgb_to_hsv(matplotlib.colors.to_rgb(c))
        ax1.plot([0, 1], [0, 0] .+ i, color = c, lw = lw)
        ax1.set_title("color")
        ax2.plot([0, 1], [0, 0] .+ i, color = string(chsv[3]), lw = lw)
        ax2.set_title("brightness")
        ax3.plot([0, 1], [0, 0] .+ i, color = string(chsv[2]), lw = lw)
        ax3.set_title("saturation")
    end
end

bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)

function add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc = 0.985, yloc = 0.975)
    bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)
    for (i, ax) in enumerate(axs)
        l = collect('a':'z')[i]
        ax.text(xloc, yloc, "$(l)"; transform = ax.transAxes, bbox = bbox, zorder = 99)
    end
end

function coolfill!(ax, x, y, dy, c, label = "")
    α = 0.25
    ax.plot(x, y, label = label, color = c, lw = 2.0)
    ax.fill_between(x, y .- dy, y .+ dy, alpha = α, color = c)
    lw = 0.5
    α2 = 0.5
    ax.plot(x, y .+ dy,  color = c, lw = lw, alpha = α2)
    ax.plot(x, y .- dy,  color = c, lw = lw, alpha = α2)
end

function nice_arrow!(ax, xc, yc, xspan, yspan;
    style = "<->", tex = "", xo = 0.2, yo = -0.2)
    ax.annotate("",  xy=(xc-xspan/2, yc - yspan/2), xycoords="data",
                xytext=(xc+xspan/2, yc + yspan/2), textcoords="data",
                arrowprops = (Dict(:arrowstyle=>style,
                                :connectionstyle=>"arc3",
                                :lw=>1.5, :facecolor => "black")), zorder = 99)
    if tex != ""
        ax.text(xc + xo, yc + yo, tex, size = 24)
    end
end

function coolhist!(ax, data, bins, color, label = "", alpha = 0.25)
    h, b, = ax.hist(data, bins, density = true, color = color,
    alpha = alpha)

    b = 0.5(b[1:end-1] .+ b[2:end])
    ax.plot(b, h, color = color, lw = 1.0, label = label)
end

function add_grid!(ax, nx::Int, ny = nx; kwargs...)
    @assert n ≥ 3
    x = ax.get_xlim()
    y = ax.get_ylim()
    dx = (x[2]-x[1])/nx; dy = (y[2]-y[1])/ny
    for i in 1:(n-1)
        ax.axhline(y[1]+n*dy, color = "gray", alpha = 0.5, kwargs...)
        ax.axvline(x[1]+n*dx, color = "gray", alpha = 0.5, kwargs...)
    end
end

"""
    axis_zoomin!(zoomin, origin, zbox, rbox, co = "C0"; kwargs...)
Create a zoomin box connecting two axes, the `zoomin` and `origin`.
The zoom-in box is in the `origin` axis, while the `zoomin` axis is the
box. `rbox` is the enclosing box of the `zoomin` axes, while `zbox`
is the small "zoom-in" box of the `origin` axis. They must be in the form
((x1, y1), (x2, y2)). `co` is color `α` the alpha setting.
"""
function axis_zoomin!(zoomin, origin, zbox, rbox, co = "C0";
    connect_lines = true, lw = 2.0, α = 1.0)
    # plot box in zoomin axis
    line, = zoomin.plot(
    [rbox[1][1], rbox[2][1], rbox[2][1], rbox[1][1], rbox[1][1]],
    [rbox[1][2], rbox[1][2], rbox[2][2], rbox[2][2], rbox[1][2]],
    color=co,  lw = lw, alpha = α)
    line.set_clip_on(false)

    line, = origin.plot(
    [zbox[1][1], zbox[2][1], zbox[2][1], zbox[1][1], zbox[1][1]],
    [zbox[1][2], zbox[1][2], zbox[2][2], zbox[2][2], zbox[1][2]],
    color=co,  lw = lw, alpha = α)
    line.set_clip_on(false)

    if connect_lines
    for e in 1:2
        con = matplotlib.patches.ConnectionPatch(
        xyA = (rbox[1][1], rbox[e][2]), xyB=(zbox[2][1], zbox[e][2]),
        coordsA="data", coordsB="data",
        axesA = zoomin, axesB=origin, color=co, lw = lw, alpha = α)
        con.set_clip_on(false)
        zoomin.add_artist(con)
    end
    end
end

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
using ConfidenceWeighted
using Base.Test

using .SoftConfidenceWeighted

scw = SCWParameter(1.0, 1.0)

@test classify(scw, Dict()) == -1

@test_throws ErrorException update(scw, Dict(), -2)
@test_throws ErrorException update(scw, Dict(), 0)
@test_throws ErrorException update(scw, Dict(), 2)

@test scw == update(scw, Dict(), -1)
@test scw == update(scw, Dict(), 1)


text:
"""
Polynomial mutation as presented in the paper:
    Deb and Deb (2012), "Analyzing Mutation Schemes for Real-Parameter Genetic Algorithms"
"""
struct PolynomialMutation{SS<:SearchSpace} <: GibbsMutationOperator
    ss::SS
    η::Float64

    PolynomialMutation(ss::SS, η = 50.0) where {SS<:SearchSpace} = new{SS}(ss, η)
    PolynomialMutation(ss::SS, options::Parameters) where {SS<:SearchSpace} =
        new{SS}(ss, options[:PM_η])
end

"""
Default parameters for `PolynomialMutation`.
"""
const PM_DefaultOptions = ParamsDict(
    :PM_η => 50.0,
)

function apply(m::PolynomialMutation, v::Number, dim::Int, target_index::Int)
    u = 2.0*rand()
    if u <= 1.0
        ΔL = u^(1.0/(1.0+m.η)) - 1.0
        return (v + ΔL * (v - mins(m.ss)[dim]))
    else
        ΔR = 1.0 - (2.0-u)^(1.0/(1.0+m.η))
        return (v + ΔR * (maxs(m.ss)[dim] - v))
    end
end


text:
__precompile__()
module ReinforcementLearning

using POMDPs
using POMDPToolbox
using Discretizers

import POMDPs: isterminal, discount
import POMDPs: action, actions, action_index, n_actions
import POMDPs: states, state_index, n_states

export

    isterminal,
    actions,
    step!,
    reset!,

    DiscretizedEnvironment

include("abstract.jl")
include("environments.jl")

end # module


text:
@testset "Set parameters" begin
    re = new(:Re)
    @parameters r, R, T
    set_param!(re,r,1.0)
    set_param!(re,R,8.314)
    set_param!(re,T,310.0)

    @test default_value(re,r) == 1.0
    @test default_value(re,R) == 8.314
    @test default_value(re,T) == 310.0
end

@testset "Incompatible parameter fail" begin
    re = new(:Re)
    @parameters s
    @test_throws ErrorException set_param!(re,s,1.0)
end

@testset "Set initial conditions" begin
    c = new(:C)
    @variables q(t)
    set_initial_value!(c,q,2.0)
    @test default_value(c,q) == 2.0
end

@testset "Missing state variable fail" begin
    c = new(:C)
    @parameters q
    @variables p(t)
    @test_throws ErrorException set_initial_value!(c,q,2.0)
    @test_throws ErrorException set_initial_value!(c,p,2.0)
end

@testset "Simulate RC circuit" begin
    r = new(:R)
    c = new(:C)
    bg = BondGraph(:RC)

    add_node!(bg, [c, r])
    connect!(bg, r, c)

    @parameters R, C
    @variables q(t)
    set_param!(r,R,2.0)
    set_param!(c,C,1.0)
    set_initial_value!(c,q,10.0)

    f(x,a,τ) = a*exp(-x/τ)

    tspan = (0.0,10.0)
    sol = simulate(bg,tspan)
    for t in [0.5,1.0,5.0,10.0]
        @test isapprox(sol(t)[1], f(t,10,2), atol=1e-5)
    end

    sol = simulate(bg,tspan; u0=[5.0])
    for t in [0.5,1.0,5.0,10.0]
        @test isapprox(sol(t)[1], f(t,5,2), atol=1e-5)
    end

    sol = simulate(bg,tspan; pmap=[1.0,3.0])
    for t in [0.5,1.0,5.0,10.0]
        @test isapprox(sol(t)[1], f(t,10,3), atol=1e-5)
    end
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
export argmin_x

"""
solve the x-minimization step in PARSDMM, i.e., x-minimization w.r.t. the augmented Lagrangian
"""
function argmin_x(
                  Q               ::Union{Array{TF,2},SparseMatrixCSC{TF,Int64},joAbstractLinearOperator{TF, TF}},
                  rhs             ::Vector{TF},
                  x               ::Vector{TF},
                  x_solve_tol_ref ::TF,
                  i               ::Integer,
                  log_PARSDMM,
                  Q_offsets=[],
                  Ax_out=zeros(TF,length(x)) ::Vector{TF},
                  comp_grid=[]
                  ) where {TF<:Real}

  #Initialize
  flag    = 0
  reslres = 0.0
  iter    = 0

  if typeof(Q)==Array{TF,2} #set up multi-threaded matrix-vector product in compressed diagonal storage format

    # w = CDS2stencil(Q, Q_offsets, comp_grid.n)
    # w = w |> gpu
    # cdims = DenseConvDims(reshape(x, comp_grid.n ...,1,1), w, padding=Int(size(w,1)-1)/2)

    Af1(in) =  Ax_CDS_MT(in, Ax_out, Q, Q_offsets)
    #Af1(in) =  Ax_stencil_cuda(in, Ax_out, comp_grid, w, cdims)

    #determine what relative residual CG needs to reach
    if i<3 #i is the PARSDMM iteration counter
      x_solve_tol_ref = TF(max(0.1*norm(Af1(x)-rhs)/norm(rhs),10*eps(TF))) #10% of current relative residual
    else
      x_solve_tol_ref = TF(min(max(0.1*norm(Af1(x)-rhs)/norm(rhs),10*eps(TF)),x_solve_tol_ref)) #10% of current relative residual, but not larger than previous relative residual
    end

    (x,flag,relres,iter) = cg(Af1,rhs,tol=x_solve_tol_ref,maxIter=1000,x=x,out=0);
    #(x,flag,relres,iter) = cg(Q,Q_offsets,rhs, x_solve_tol_ref, 1000, x, 0)

  elseif typeof(Q)==SparseMatrixCSC{TF,Int64} #CG with native Julia sparse CSC format MPVs

    #determine what relative residual CG needs to reach
    if i<3 #i is the PARSDMM iteration counter
      x_solve_tol_ref = TF(max(0.1*norm(Q*x-rhs)/norm(rhs),10*eps(TF))) #10% of current relative residual
    else
      x_solve_tol_ref = TF(min(max(0.1*norm(Q*x-rhs)/norm(rhs),10*eps(TF)),x_solve_tol_ref)) #10% of current relative residual, but not larger than previous relative residual
    end

    (x,flag,relres,iter) = cg(Q,rhs,tol=x_solve_tol_ref,maxIter=1000,x=x,out=0);

  else #typeof(Q) <: joAbstractLinearOperator{TF, TF}#CG with JOLI operators
    Af2(in) =  return Q*in
    #determine what relative residual CG needs to reach
    if i<3 #i is the PARSDMM iteration counter
      x_solve_tol_ref = TF(max(0.1*norm(Q*x-rhs)/norm(rhs),10*eps(TF))) #10% of current relative residual
    else
      x_solve_tol_ref = TF(min(max(0.1*norm(Q*x-rhs)/norm(rhs),10*eps(TF)),x_solve_tol_ref)) #10% of current relative residual, but not larger than previous relative residual
    end

    (x,flag,relres,iter) = cg(Af2,rhs,tol=x_solve_tol_ref,maxIter=1000,x=x,out=0);

  end

  x_solve_tol_ref = TF(x_solve_tol_ref)

  return x,iter,relres,x_solve_tol_ref
end

function Ax_CDS_MT(in::Vector{TF}, Ax_out::Vector{TF}, Q::Array{TF,2}, Q_offsets::Vector{Int}) where TF <: Real

  fill!(Ax_out,TF(0))
  CDS_MVp_MT(size(Q,1),size(Q,2),Q,Q_offsets,in,Ax_out);

  return Ax_out
end

# function Ax_stencil_cuda(in::Vector{TF}, Ax_out::Vector{TF}, comp_grid, w::CUDA.CuArray{TF, 5}, cdims) where TF <: Real
#   in     = in |> gpu
#   Ax_out = Ax_out |>gpu
#   fill!(Ax_out,TF(0))
#   Ax_out = reshape(Ax_out,comp_grid.n ...,1,1)
#   in     = reshape(in,comp_grid.n ...,1,1)
#   conv!(Ax_out, in, w, cdims)

#   # #for now, boundary conditions are an issue, manually correct it
#   # kernel_size = size(w,1)

#   Ax_out = Ax_out |> cpu

#   return vec(Ax_out)

# end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# function gen_ip_chol_quote!(qa,N,L,U = :U, S = :S, extract = extract_symbol)

#     for i ∈ 1:N
#         lti = small_triangle(i)
#         bti = lti+i
#         Ui_i = extract(U, bti)
#         U_i = extract(S, bti)
#         if i > 1
#             Uj_i = extract(U, lti+1)
#             Ui_i_expr = :($Ui_i = sqrt($U_i - $Uj_i*$Uj_i))
#             for j ∈ 2:i - 1
#                 Uj_i = extract(U, lti+j)
#                 push!(Ui_i_expr.args[2].args[2].args, :($Uj_i*$Uj_i))
#             end
#             push!(qa, Ui_i_expr)
#         else
#             push!(qa, :($Ui_i = sqrt($U_i)))
#         end

#         for j ∈ i+1:N
#             ltj = small_triangle(j)
#             Uj_i = extract(U, ltj+i)
#             U_i = extract(S, ltj+i)
#             if i > 1
#                 Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
#                 Uj_i_expr = :($Uj_i = ($U_i - $Ujk * $Uik) / $Ui_i)
#                 for k ∈ 2:i - 1
#                     Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
#                     push!(Uj_i_expr.args[2].args[2].args, :($Ujk * $Uik))
#                 end
#                 push!(qa, Uj_i_expr )
#             else
#                 push!(qa, :($U_i /= $Ui_i) )
#             end
#         end
#     end

#     # push!(q.args, U)
#     qa
# end

# function gen_ip_chol_quote!(qa,N,L,U = :U, S = :S, extract = extract_symbol)

#     for i ∈ 1:N
#         lti = small_triangle(i)
#         bti = lti+i
#         Ui_i = extract(U, bti)
#         U_i = extract(S, bti)
#         if i == 1
#             push!(qa, :($Ui_i = sqrt($U_i)))
#         elseif i == 2
#             Uj_i = extract(U, lti+1)
#             push!(qa, :($Ui_i = sqrt($U_i - $Uj_i*$Uj_i)))
#         else
#             Uj_i = extract(U, lti+1)
#             push!(qa, :($Ui_i = $U_i - $Uj_i*$Uj_i))
#             for j ∈ 2:i - 1
#                 Uj_i = extract(U, lti+j)
#                 push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
#             end
#             push!(qa, :($Ui_i = sqrt($Ui_i)))
#         end

#         for j ∈ i+1:N
#             ltj = small_triangle(j)
#             Uj_i = extract(U, ltj+i)
#             U_i = extract(S, ltj+i)
#             push!(qa, :($Uj_i = $U_i))
#             for k ∈ 1:i - 1
#                 Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
#                 push!(qa, :($Uj_i -= $Ujk * $Uik))
#             end
#             push!(qa, :($Uj_i /= $Ui_i) )
#         end
#     end

#     # push!(q.args, U)
#     qa
# end
function gen_ip_chol_quote!(qa,::Type{T},N::Int,L,U = :U, S = :S, extract = extract_symbol) where T

    chunksize = 32 ÷ sizeof(T)
    # @show qa

    for i ∈ 1:N
        lti = small_triangle(i)
        bti = lti+i
        Ui_i = extract(U, bti)
        U_i = extract(S, bti)
        if i == 1
            push!(qa, :($Ui_i = sqrt($U_i)))
        else
            nchunks, r = divrem(i-1, chunksize)
            if r == 0
                push!(qa, :($Ui_i = $U_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1:chunksize ]...)  ))
                for chunk ∈ 2:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                end

            else
                push!(qa, :($Ui_i = $U_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1:r ]...)  ))
                for chunk ∈ 1:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                end
            end
            # for j ∈ 1:i - 1
            #     Uj_i = extract(U, lti+j)
            #     push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
            # end
            push!(qa, :($Ui_i = sqrt($Ui_i)))
        end

        for j ∈ i+1:N
            ltj = small_triangle(j)
            Uj_i = extract(U, ltj+i)
            U_i = extract(S, ltj+i)
            if i == 1
                push!(qa, :($Uj_i = $U_i / $Ui_i) )
            else

                nchunks, r = divrem(i-1, chunksize)
                if r == 0
                    push!(qa, :($Uj_i = $U_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1:chunksize ]...)  ))
                    for chunk ∈ 2:nchunks
                        push!(qa, :($Uj_i= $Uj_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                    end

                else
                    push!(qa, :($Uj_i = $U_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1:r ]...)  ))
                    for chunk ∈ 1:nchunks
                        push!(qa, :($Uj_i = $Uj_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                    end
                end

                # push!(qa, :($Uj_i = $U_i))
                # for k ∈ 1:i - 1
                #     Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
                #     push!(qa, :($Uj_i -= $Ujk * $Uik))
                # end
                push!(qa, :($Uj_i /= $Ui_i) )
            end
        end
    end

    # push!(q.args, U)
    qa
end
function gen_ip_choldet_quote!(qa,::Type{T},N::Int,L,U = :U, S = :S, extract = extract_symbol) where T

    chunksize = 32 ÷ sizeof(T)
    # @show qa
    push!(qa, :(root_det = $(one(T))))
    for i ∈ 1:N
        lti = small_triangle(i)
        bti = lti+i
        Ui_i = extract(U, bti)
        U_i = extract(S, bti)
        if i == 1
            push!(qa, :($Ui_i = sqrt($U_i)))
        else
            nchunks, r = divrem(i-1, chunksize)
            if r == 0
                push!(qa, :($Ui_i = $U_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1:chunksize ]...)  ))
                for chunk ∈ 2:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                end

            else
                push!(qa, :($Ui_i = $U_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1:r ]...)  ))
                for chunk ∈ 1:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(extract(U, lti+j)) * $(extract(U, lti+j))) for j ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                end
            end
            # for j ∈ 1:i - 1
            #     Uj_i = extract(U, lti+j)
            #     push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
            # end
            push!(qa, :($Ui_i = sqrt($Ui_i)))
            push!(qa, :(root_det *= $Ui_i))
        end

        for j ∈ i+1:N
            ltj = small_triangle(j)
            Uj_i = extract(U, ltj+i)
            U_i = extract(S, ltj+i)
            if i == 1
                push!(qa, :($Uj_i = $U_i / $Ui_i) )
            else

                nchunks, r = divrem(i-1, chunksize)
                if r == 0
                    push!(qa, :($Uj_i = $U_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1:chunksize ]...)  ))
                    for chunk ∈ 2:nchunks
                        push!(qa, :($Uj_i= $Uj_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                    end

                else
                    push!(qa, :($Uj_i = $U_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1:r ]...)  ))
                    for chunk ∈ 1:nchunks
                        push!(qa, :($Uj_i = $Uj_i - +$([:( $(extract(U, ltj+k)) * $(extract(U, lti+k))) for k ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                    end
                end

                # push!(qa, :($Uj_i = $U_i))
                # for k ∈ 1:i - 1
                #     Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
                #     push!(qa, :($Uj_i -= $Ujk * $Uik))
                # end
                push!(qa, :($Uj_i /= $Ui_i) )
            end
        end
    end

    # push!(q.args, U)
    qa
end
# function gen_ip_chol_quote!(qa,N,L,U = :U, S = :S, extract = extract_symbol)

#     for i ∈ 1:N
#         lti = small_triangle(i)
#         bti = lti+i
#         Ui_i = extract(U, bti)
#         U_i = extract(S, bti)
#         push!(qa, :($Ui_i = $U_i))
#         for j ∈ 1:i - 1
#             Uj_i = extract(U, lti+j)
#             push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
#         end
#         push!(qa, :($Ui_i = sqrt($Ui_i)))

#         for j ∈ i+1:N
#             ltj = small_triangle(j)
#             Uj_i = extract(U, ltj+i)
#             U_i = extract(S, ltj+i)
#             push!(qa, :($Uj_i = $U_i))
#             for k ∈ 1:i - 1
#                 Ujk, Uik = extract(U, ltj+k), extract(U, lti+k)
#                 push!(qa, :($Uj_i -= $Ujk * $Uik))
#             end
#             push!(qa, :($Uj_i /= $Ui_i) )
#         end
#     end

#     # push!(q.args, U)
#     qa
# end
function gen_extracted_ip_chol_quote!(qa, ::Type{T}, N, L, U = :U, S = :S) where T
    for i ∈ 1:L
        push!(qa, :($(id_symbol(S, i)) = $(S)[$i]))
    end
    gen_ip_chol_quote!(qa, T, N, L, U, S, id_symbol)
    for i ∈ 1:L
        push!(qa, :($(U)[$i] = $(id_symbol(U, i))))
    end
end
# function chol_quote(S = :S, ::Type{T} = Float64) where T
#     q, qa = create_quote()
#     for i ∈ 1:L
#         push!(qa, :($(id_symbol(S, i)) = $(S)[$i]))
#     end
#     gen_ip_chol_quote!(qa, T, N, L, :U, :S, id_symbol)
#     q
# end
@generated function LinearAlgebra.chol(S::AbstractSymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    gen_extracted_ip_chol_quote!(qa, T, N, L, :S, :S)
    push!(q.args, :(StaticUpperTriangularMatrix{$T,$N,$L}( ( @ntuple $L S )  )))
    q
end
@generated function cholesky!(S::SymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    gen_extracted_ip_chol_quote!(qa, T, N, L, :S, :S)
    push!(q.args, :(UpperTriangularMatrix{$T,$N,$L}(S.data)))
    q
end
@generated function cholesky!(U::UpperTriangularMatrix{T,N,L}, S::AbstractSymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    gen_extracted_ip_chol_quote!(qa, T, N, L, :U, :S)
    push!(q.args, :U)
    q
end

@generated function choldet!(S::SymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    extract_linear!(qa, L, :S)
    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    insert_linear!(qa, L, :S)
    push!(q.args, :(root_det))
    q
end
@generated function choldet!(S::LinearAlgebra.Symmetric{T,RecursiveMatrix{T,N,N,N2}}) where {T,N,N2}
    q, qa = create_quote()
    L = big_triangle(N)
    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc] ))
    end

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl

    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( S.data[$nr, $nc] = $(Symbol(:S_, ind)) ))
    end
    push!(q.args, :(root_det))
    q
end
@generated function choldet!(S::SymmetricMatrix{T,M,L2}, ::Val{N}) where {T,M,N,L2}
    @assert N <= M
    L = big_triangle(N)
    q, qa = create_quote()
    extract_linear!(qa, L, :S)
    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    insert_linear!(qa, L, :S, :S)
    push!(q.args, :(root_det))
    q
end
@generated function choldet!(S::LinearAlgebra.Symmetric{T,RecursiveMatrix{T,M,M,N2}}, ::Val{N}) where {T,M,N,N2}
    @assert N <= M
    q, qa = create_quote()
    L = big_triangle(N)
    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc]) )
    end

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl

    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( S.data[$nr, $nc] = $(Symbol(:S_, ind)) ))
    end
    push!(q.args, :(root_det))
    q
end
@generated function choldet!(Si::Union{SymmetricMatrix{T,N,L1},Symmetric{T,RecursiveMatrix{T,N,N,L1}}},
                            S::Union{SymmetricMatrix{T,N,L2},Symmetric{T,RecursiveMatrix{T,N,N,L2}}}) where {T,N,L1,L2}
    q, qa = create_quote()
    if S <: Symmetric
        L_triangle = big_triangle(N)
        ind = 0
        for nc ∈ 1:N, nr ∈ 1:nc
            ind += 1
            push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc]) )
        end
    else
        L_triangle = L1
        extract_linear!(qa, L1, :S)
    end

    gen_ip_choldet_quote!(qa, T, N, L_triangle, :S, :S, id_symbol) #in this file, cholesky.jl #creates root_det variable, root of the determinant of S.

    if Si <: Symmetric
        ind = 0
        for nc ∈ 1:N, nr ∈ 1:nc
            ind += 1
            push!(qa, :( Si.data[$nr, $nc] = $(Symbol(:S_, ind)) ))
        end
    else
        insert_linear!(qa, L_triangle, :S)
    end
    push!(q.args, :(root_det))
    q
end
@generated function inv!(S::SymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    extract_linear!(qa, L, :S)

    gen_ip_chol_quote!(qa, T, N, L, :S, :S, id_symbol)

    insert_linear!(qa, L, :S)
    push!(q.args, :(S))
    q
end
@generated function invdet!(S::SymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    extract_linear!(qa, L, :S)

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    gen_ip_inv_quote!(qa, N, L, :S, :S, id_symbol) #in inverse.jl # triangular matrix inversion
    gen_upper_xxt_quote!(qa,T,N,L, :S, :Si, id_symbol) #in arithmetic.jl

    insert_linear!(qa, L, :S, :Si)
    push!(q.args, :(root_det))
    q
end
@generated function invdet!(S::SymmetricMatrix{T,M,L2}, ::Val{N}) where {T,N,L2,M}
    @assert N <= M
    L = big_triangle(N)
    q, qa = create_quote()
    extract_linear!(qa, L, :S)

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    gen_ip_inv_quote!(qa, N, L, :S, :S, id_symbol) #in inverse.jl # triangular matrix inversion
    gen_upper_xxt_quote!(qa,T,N,L, :S, :Si, id_symbol) #in arithmetic.jl

    insert_linear!(qa, L, :S, :Si)
    push!(q.args, :(root_det))
    q
end
@generated function invdet!(S::LinearAlgebra.Symmetric{T,RecursiveMatrix{T,N,N,N2}}) where {T,N,N2}
    q, qa = create_quote()
    L = big_triangle(N)
    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc]) )
    end
    # extract_linear!(qa, L, :S)

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    gen_ip_inv_quote!(qa, N, L, :S, :S, id_symbol) #in inverse.jl # triangular matrix inversion
    gen_upper_xxt_quote!(qa,T,N,L, :S, :Si, id_symbol) #in arithmetic.jl

    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( S.data[$nr, $nc] = $(Symbol(:Si_, ind)) ))
    end
    push!(q.args, :(root_det))
    q
end
@generated function invdet!(S::LinearAlgebra.Symmetric{T,RecursiveMatrix{T,M,M,N2}}, ::Val{N}) where {T,N,N2,M}
    @assert N <= M
    q, qa = create_quote()
    L = big_triangle(N)
    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc]) )
    end
    # extract_linear!(qa, L, :S)

    gen_ip_choldet_quote!(qa, T, N, L, :S, :S, id_symbol) #in this file, cholesky.jl
    gen_ip_inv_quote!(qa, N, L, :S, :S, id_symbol) #in inverse.jl # triangular matrix inversion
    gen_upper_xxt_quote!(qa,T,N,L, :S, :Si, id_symbol) #in arithmetic.jl

    ind = 0
    for nc ∈ 1:N, nr ∈ 1:nc
        ind += 1
        push!(qa, :( S.data[$nr, $nc] = $(Symbol(:Si_, ind)) ))
    end
    push!(q.args, :(root_det))
    q
end
@generated function invdet!(Si::Union{SymmetricMatrix{T,N,L1},Symmetric{T,RecursiveMatrix{T,N,N,L1}}},
                            S::Union{SymmetricMatrix{T,N,L2},Symmetric{T,RecursiveMatrix{T,N,N,L2}}}) where {T,N,L1,L2}
    q, qa = create_quote()
    if S <: Symmetric
        L_triangle = big_triangle(N)
        ind = 0
        for nc ∈ 1:N, nr ∈ 1:nc
            ind += 1
            push!(qa, :( $(Symbol(:S_, ind)) = S.data[$nr, $nc]) )
        end
    else
        L_triangle = L1
        extract_linear!(qa, L1, :S)
    end

    gen_ip_choldet_quote!(qa, T, N, L_triangle, :S, :S, id_symbol) #in this file, cholesky.jl #creates root_det variable, root of the determinant of S.
    gen_ip_inv_quote!(qa, N, L_triangle, :S, :S, id_symbol) #in inverse.jl # triangular matrix inversion
    gen_upper_xxt_quote!(qa,T,N, L_triangle, :S, :Si, id_symbol) #in arithmetic.jl

    if Si <: Symmetric
        ind = 0
        for nc ∈ 1:N, nr ∈ 1:nc
            ind += 1
            push!(qa, :( Si.data[$nr, $nc] = $(Symbol(:Si_, ind)) ))
        end
    else
        insert_linear!(qa, L_triangle, :Si)
    end
    push!(q.args, :(root_det))
    q
end
# @generated function cholesky!(U::UpperTriangularMatrix{T,N,L}, S::AbstractSymmetricMatrix{T,N,L}) where {T,N,L}
#     q, qa = create_quote()
#     gen_extracted_ip_chol_quote!(qa, T, N, L, :U, :S)
#     push!(q.args, :U)
#     q
# end

function gen_ip_revchol_quote!(qa, ::Type{T}, N,L, U = :U, S = :S, extract = extract_symbol) where T

    for i ∈ N:-1:1
        lti = small_triangle(i)
        bti = lti+i
        Ui_i = extract(U, bti)
        U_i = extract(S, bti)
        if i == N
            push!(qa, :($Ui_i = sqrt($U_i)))
        else
            push!(qa, :($Ui_i = $U_i))
            for j ∈ i+1:N
                Uj_i = extract(U, small_triangle(j) + i)
                push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
            end
            push!(qa, :($Ui_i = sqrt($Ui_i)))
        end

        for j ∈ 1:i-1
            ltj = small_triangle(j)
            Uj_i = extract(U, lti+j)
            U_i  = extract(S, lti+j)
            if i == N
                push!(qa, :($Uj_i = $U_i / $Ui_i) )
            else
                push!(qa, :($Uj_i = $U_i))
                for k ∈ i+1:N
                    ltk = small_triangle(k)
                    Ujk, Uik = extract(U, ltk+j), extract(U, ltk+i)
                    push!(qa, :($Uj_i -= $Ujk * $Uik))
                end
                push!(qa, :($Uj_i /= $Ui_i) )
            end
        end
    end
    # push!(q.args, :($T( ( @ntuple $L Ui ) )))
    qa

end

"""
This directly calculates the inverse of the Cholesky decomposition of the inverse.
Or, for positive definite input matrix S, it calculates the upper triangular matrix U such that
S = U * U'
"""
@generated function reverse_cholesky(S::SymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    for i ∈ 1:L
        push!(qa, :($(id_symbol(:S, i)) = S[$i]))
    end
    gen_ip_revchol_quote!(qa, T, N, L, :U, :S, id_symbol)
    push!(q.args, :(UpperTriangularMatrix{$T,$N,$L}( ( @ntuple $L U ) )))
    q
end

@generated function reverse_cholesky!(U::UpperTriangularMatrix{T,N,L}, S::AbstractSymmetricMatrix{T,N,L}) where {T,N,L}
    q, qa = create_quote()
    for i ∈ 1:L
        push!(qa, :($(id_symbol(:S, i)) = S[$i]))
    end
    gen_ip_revchol_quote!(qa, T, N, L, :U, :S, id_symbol)
    for i ∈ 1:L
        push!(qa, :(U[$i] = $(id_symbol(:U, i))))
    end
    push!(q.args, :U)
    q
end


text:
function balances(dx, x, m, t)

    # get data from the m (model) dictionary -
    S = m["S"]
    number_of_dynamic_states = m["number_of_dynamic_states"]

    # compute the kinetics - powerlaw
    rV = compute_kinetics_vector(x,m,t)

    # compute the rhs -> store in a temp vector
    tmp = S*rV

    # populate the dx vector with the tmp vector -
    for i ∈ 1:number_of_dynamic_states
        dx[i] = tmp[i]
    end
    abcd = dx[13]
end


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
pyplot()
Random.seed!(11)
using BSON: @save
include("../scripts/paper_settings.jl");
# include("../../utils/sindy_turing_utils.jl")
# include("../../utils/main_algorithm_utils.jl")
# include("../../utils/sindy_statistics_utils.jl")

################################################################################
#####
##### SIR example: get data
#####
include(EXAMPLE_FOLDERNAME*"example_fitzHughNagumo.jl")

# Define the multivariate forcing function
num_ctr = 3;
    U_starts = rand(1, num_ctr) .* tspan[2]/2
    U_widths = 0.6;#0;
    amplitude = 1.0
my_U_func_time2(t) = U_func_time(t, u0,
                        U_widths, U_starts,
                        F_dim=1,
                        amplitude=amplitude)

# Get data
sol = solve_fhn_system(U_func_time=my_U_func_time2)
dat = Array(sol)
numerical_grad = numerical_derivative(dat, ts)
true_grad = core_dyn_true(dat)

U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time2(t)
end

## Also get baseline true/ideal cases
# Uncontrolled
dat_raw = Array(solve_fhn_system())
numerical_grad_raw = numerical_derivative(dat_raw, ts)
true_grad_raw = zeros(size(dat))

# Intialize truth object
this_truth = sra_truth_object(dat_raw, true_grad, U_true, core_dyn_true)

#####
##### Build SRA object and analyze
#####
# Initialize
this_model = sra_stateful_object(ts, tspan, dat, u0, numerical_grad)
prams = this_model.parameters
prams.sindyc_ensemble_parameters[:selection_criterion] =
    sindy_cross_validate;
# Several changes are required because there are only 2 variables here
prams.variable_names = ["v", "ω"];
prams.sindy_terms_list = Iterators.product(1:4, 1:4)
prams.sindy_library["cross_terms"] = [2, 3] # Also include cubic terms

fit_first_model(this_model, 1);
print_true_equations(this_truth)
print_current_equations(this_model)

plot_sindy_model(this_model, 1)

#################################################################
### Iterate
calculate_subsampled_ind(this_model);
all_models = fit_model(this_model);

print_true_equations(this_truth)
print_current_equations(this_model)

# plot_subsampled_points(this_model)
plot_subsampled_simulation(this_model, 1)
plot_subsampled_derivatives(this_model, 1)
# plot_residual(this_model)

# this_model.U_true

# plot_subsampled_points_and_control(this_model, this_truth)

#####
##### Save FitzHugh data
#####
this_dat_name = "dat_library_of_examples_fhn_"

save_for_plotting(this_model, this_truth, this_dat_name)


text:
"""
Modia module for elaboration of models.

* Developer: Hilding Elmqvist, Mogram AB
* First version: July-August 2016
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)
"""
module ModelElaboration

const PrintOriginalModel = false # Needed for the moment to set partial attribute.
const PrintInstantiated = false
const PrintFlattened = true
const PrintElaborated = true
const PrintJSON = false
const fileStdOut = false

# ----------------------------------------------------------------------------

#using Debug
using Zygote
using ..Instantiation
using ..Execution
using Base.Meta: quot, isexpr
using Base: isapprox
using DataStructures
import ..Instantiation: GetField, This, Der, Symbolic, vars_of, Instance
import ..Execution: logTiming
using ..Synchronous
using ..BLTandPantelides
using ..BLTandPantelidesUtilities
using ..BasicStructuralTransform
using ..StructuralTransform
using ..SymbolicTransform
using ..Utilities
using ..ModiaLogging
#= Temporarily removed due to problem with PyPlot
using PyPlot
=#
import ModiaMath

using JSON

@static if VERSION < v"0.7.0-DEV.2005"
else
    using Printf
end

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

export elaborate, prettyPrint, simulateModel, simulate, simulate_der, skewCoords, skew, residue, residue_der, modiaCross, showExpr, transformModel, simulateMultiModeModel, allInstances
# export modiaSwitches, defineSwitch, setSwitch, showSwitches, getSwitch,
export checkSimulation

const useIncidenceMatrix = BasicStructuralTransform.useIncidenceMatrix
const elaborate = true
#@show elaborate

const noResult = Dict{Symbol,AbstractArray{T,1} where T}()

allInstances(s) = return s

function dummyDer(e)
    f = e
    if typeof(f) == Instantiation.GetField
        f = f.name
    end

    if typeof(f) == Instantiation.Der
        f = dummyDer(f.base)
    end
    s = string("der_", f)
    return Symbol(s)
end

function dummyState(e)
    f = e
    if typeof(f) == Instantiation.GetField
        f = f.name
    end

    if typeof(f) == Instantiation.Der
        f = dummyState(f.base)
    end
    string(f)
end

function substituteDer(ex, dummyDerivatives, unassigned)
    if typeof(ex) == Instantiation.Der
        newname = dummyDer(ex.base)
        dumSt = dummyState(ex)
        unass = map(string, unassigned)

        if dumSt in unass
            # println("substituteDer for unassigned: ", newname)
            ex = GetField(This(), newname)
        elseif newname in dummyDerivatives
            #  println("substituteDer: ", newname)
            ex = GetField(This(), newname)
        else
            ex
        end

    elseif typeof(ex) == Expr
        Expr(ex.head, [substituteDer(arg, dummyDerivatives, unassigned) for arg in ex.args]...)
    else
        ex
    end
end

function elaborateModel(flat_model)
    loglnModia("elaborate:")
    solved_model = flat_model

    #  print("Elaborate:          ")
    #  @time
    (equations, variables, assign, components, Avar, Bequ, states, deriv, unassignedNames, incidenceMatrix, varSizes, varTypes, equSizes, equTypes) = StructuralTransform.transformStructurally(flat_model)

    if equations == nothing # Connector
        return
    end

    flat_model.equations = equations

    loglnModia("\nSYMBOLIC PROCESSING")

    if logTiming
        print("Symbolic processing:   ")
        @time solvedComponents, algebraic = differentiateSortedEquations(equations, variables, components, assign, Avar, Bequ)
    else
        solvedComponents, algebraic = differentiateSortedEquations(equations, variables, components, assign, Avar, Bequ)
    end

    if length(SymbolicTransform.dummyDerivatives.keys) > 0
        printSymbolList("\nDummy derivatives", SymbolicTransform.dummyDerivatives.keys)
    end

    unassignedVar = []
    solved_model.equations = []
    initSynchronousCounters()
    previousVariables = []

    for c in solvedComponents
        for (e, solved) in c
            e = substituteClocked(previousVariables, e)
            e = substituteDer(e, SymbolicTransform.dummyDerivatives.keys, unassignedVar)

            if !solved # typeof(e.args[1]) == Instantiation.Der || length(c) > 1   #
                # Use equation
                push!(solved_model.equations, Expr(:(=), e.args[1], e.args[2]))
            else
                # Use assignment
                push!(solved_model.equations, Expr(:(:=), e.args[1], e.args[2]))
            end
        end
    end

    # solved_model.variables = [solved_model.variables; SymbolicTransform.dummyDerivatives]

    for d in SymbolicTransform.dummyDerivatives.keys
    # Remove der_
        v = string(d)
        if length(v) > 4 && v[1:4] == "der_"
            v = v[5:end]
        end

        if length(v) > 4 && v[1:4] == "der_"
            v = v[5:end]
        end

        vsize = varSizes[Symbol(v)]
        st = if vsize == (); 0.0 else zeros(vsize) end
        solved_model.variables[d] = Instantiation.Variable(start=st)
    end

    #= For test:
    for d in map(string, unassignedVar)
        @show Symbol(d)
        solved_model.variables[Symbol("der_"*d)] = Instantiation.Variable()  # add typ
        solved_model.variables[Symbol("der_der_"*d)] = Instantiation.Variable()
    end
    =#

    solved_model.initial_pre = []
    solved_model.F_pre = []
    solved_model.F_post = []

    for p in previousVariables
        push!(solved_model.initial_pre, :($(p.name) = 0))
        push!(solved_model.F_pre, :($(p.name) = 0))
        push!(solved_model.F_post, :(dummy = $(Synchronous.updatePrevious)($(p.name), $(SymbolicTransform.simulationModel_symbol), 1)))
    end

    if PrintElaborated
        loglnModia("\nSTRUCTURALLY AND SYMBOLICALLY PROCESSED MODEL")
        loglnModia("  Sorted equations, := used for solved algebraic variables")
        showInstance(solved_model)
        loglnModia()
    end

    return solved_model
end

# Handle allInstances operator
function allInstances!(all, instName, inst, var, class)
    for key in keys(inst.variables)
        if isa(inst.variables[key], Instance)
            allInstances!(all, key, inst.variables[key], var, class)
        else
            if inst.model_name == class && key == var
                push!(all, instName)
            end
        end
    end
end

function substituteAllInstances(ex, modified_model, class, flat)
    if typeof(ex) == Expr && ex.head == :call && ex.args[1] == allInstances
        all = []
        name = string(ex.args[2].name)

@static if VERSION < v"0.7.0-DEV.2005"
        name = name[rsearchindex(name, ".") + 1:end]
else
        name = name[first(something(findlast(".", name), 0:-1)) + 1:end]
end

        name = Symbol(name)

        if !flat
            loglnModia("Found name: $name in $class")
            Expr(ex.head, ex.args..., class)
        else
            allInstances!(all, :(), modified_model, name, ex.args[3])
            loglnModia("all = $all")
            ref = [GetField(This(), Symbol(string(b) * "." * string(name))) for b in all]
            loglnModia("ref = $ref")
            Expr(:vect, ref...)
        end

    elseif typeof(ex) == Expr
        Expr(ex.head, [substituteAllInstances(arg, modified_model, class, flat) for arg in ex.args]...)
    else
        ex
    end
end

function traverseAndSubstituteAllInstances(instName, modified_model, inst, flat)
    for key in keys(inst.variables)
        if isa(inst.variables[key], Instance)
            traverseAndSubstituteAllInstances(key, modified_model, inst.variables[key], flat)
        end
    end

    newEqu = []
    for e in inst.equations
        newE = substituteAllInstances(e, modified_model, inst.model_name, flat)
        if newE != e
            # @show e
            # @show newE
    end
        push!(newEqu, newE)
    end
    inst.equations = newEqu
end

function simulateModelWithOptions(model, t; options=Dict())
    opt = Dict(options)
    ModiaLogging.setDefaultLogName(string(model.name))
    ModiaLogging.setOptions(opt)
    StructuralTransform.setOptions(opt)
    BasicStructuralTransform.setOptions(opt)
    Instantiation.setOptions(opt)
    Execution.setOptions(opt)

    logSimulation = false
    if haskey(opt, :logSimulation)
        logSimulation = opt[:logSimulation]
        @show logSimulation
        delete!(opt, :logSimulation)
    end

    relTol = 1e-4
    if haskey(opt, :relTol)
        relTol = opt[:relTol]
        @show relTol
        delete!(opt, :relTol)
    end

    hev = 1e-8
    if haskey(opt, :hev)
        hev = opt[:hev]
        @show hev
        delete!(opt, :hev)
    end

    if length(keys(opt)) > 0
        println("Option(s) not found: ", keys(opt))
    end

    if fileStdOut
        @static if VERSION < v"0.7.0-DEV.2005"
            originalSTDOUT = STDOUT
        else
            originalSTDOUT = stdout
        end
        (outRead, outWrite) = redirect_stdout()
    end

    openLogModia()

    start = time_ns()
    println("check")
    setDerAsFunction(false)
    if BasicStructuralTransform.logStatistics
        println("\nSimulating model: ", model.name)
    end
    loglnModia("\nSimulating model: ", model.name)

    if PrintOriginalModel
        loglnModia("ORIGINAL MODEL:")
        showModel(model)
    end

    if logTiming
        print("Instantiate:           ")
        @time modified_model = instantiate(model, first(t))
    else
        modified_model = instantiate(model, first(t))
    end
    traverseAndSubstituteAllInstances(:(), modified_model, modified_model, false)

    if PrintInstantiated
        loglnModia("INSTANTIATED MODEL:")
        showInstance(modified_model)
        loglnModia()
    end

    if PrintJSON
        name = string(modified_model.model_name)
        file = open(name * ".JSON", "w")
        printJSON(file, modified_model, name, name)
        close(file)
    end

    if logTiming
        print("Flatten:               ")
        @time flat_model = flatten(modified_model)
    else
        flat_model = flatten(modified_model)
    end
    traverseAndSubstituteAllInstances(:(), modified_model, flat_model, true)

    if PrintFlattened
        loglnModia("\nFLATTENED MODEL")
        showInstance(flat_model)
        loglnModia()
    end

    if elaborate
        solved_model = elaborateModel(flat_model)
        if solved_model == nothing
            return nothing
        end
    else
        solved_model = flat_model
    end

    #=
    #  JSON.json(t::This) = "this."

    if PrintJSONsolved
       jsonSolved = JSON.json(solved_model.equations[1])
        @show jsonSolved
    end
    =#
    @static if VERSION < v"0.7.0-DEV.2005"
        disableSimulation = false
    else
        disableSimulation = false
    end

    if !BasicStructuralTransform.newStateSelection && !disableSimulation ## Disable simulation for the moment
        loglnModia("\nSIMULATION")
        # @show incidenceMatrix
        t = vec(collect(Float64, t))
        if logTiming
            print("Code generation and simulation:         ")
            # @show solved_model t useIncidenceMatrix logSimulation
            @time res = simulate_ida(solved_model, t, if useIncidenceMatrix; incidenceMatrix else nothing end, log=logSimulation, relTol=relTol, hev=hev)
        else
            res = simulate_ida(solved_model, t, if useIncidenceMatrix; incidenceMatrix else nothing end, log=logSimulation, relTol=relTol, hev=hev)
        end
    else
        res = Dict{Symbol,AbstractArray{T,1} where T}()
    end

    setDerAsFunction(true)

    if fileStdOut
        #  close(outWrite)
        output = readavailable(outRead)
        close(outRead)
        redirect_stdout(originalSTDOUT)
        print(String(output))
    end

    if logTiming
        println(@sprintf("Total time: %0.3f", (time_ns() - start) * 1E-9), " seconds")
    end
    closeLogModia()
    println("something")
    println("noRes = $noResult")
    println("res = $res")
    #test_expr!("@test", ex, kws...)
    #orig_ex = Expr(:inert, ex)
    #result = get_test_result(ex, __source__)
    #:(do_test($result, $orig_ex))
    if res != noResult
        return res
    else
        throw(DomainError("res = $res"))
    end
    return res
end

function simulateModelWithOptionsDer(model, t, var, var_val, key; options=Dict())
    opt = Dict(options)
    ModiaLogging.setDefaultLogName(string(model.name))
    ModiaLogging.setOptions(opt)
    StructuralTransform.setOptions(opt)
    BasicStructuralTransform.setOptions(opt)
    Instantiation.setOptions(opt)
    Execution.setOptions(opt)

    logSimulation = false
    if haskey(opt, :logSimulation)
        logSimulation = opt[:logSimulation]
        @show logSimulation
        delete!(opt, :logSimulation)
    end

    relTol = 1e-4
    if haskey(opt, :relTol)
        relTol = opt[:relTol]
        @show relTol
        delete!(opt, :relTol)
    end

    hev = 1e-8
    if haskey(opt, :hev)
        hev = opt[:hev]
        @show hev
        delete!(opt, :hev)
    end

    if length(keys(opt)) > 0
        println("Option(s) not found: ", keys(opt))
    end

    if fileStdOut
        @static if VERSION < v"0.7.0-DEV.2005"
            originalSTDOUT = STDOUT
        else
            originalSTDOUT = stdout
        end
        (outRead, outWrite) = redirect_stdout()
    end

    openLogModia()

    start = time_ns()
    println("check")
    setDerAsFunction(false)
    if BasicStructuralTransform.logStatistics
        println("\nSimulating model: ", model.name)
    end
    loglnModia("\nSimulating model: ", model.name)

    if PrintOriginalModel
        loglnModia("ORIGINAL MODEL:")
        showModel(model)
    end

    if logTiming
        print("Instantiate:           ")
        @time modified_model = instantiate(model, first(t))
    else
        modified_model = instantiate(model, first(t))
    end
    traverseAndSubstituteAllInstances(:(), modified_model, modified_model, false)

    if PrintInstantiated
        loglnModia("INSTANTIATED MODEL:")
        showInstance(modified_model)
        loglnModia()
    end

    if PrintJSON
        name = string(modified_model.model_name)
        file = open(name * ".JSON", "w")
        printJSON(file, modified_model, name, name)
        close(file)
    end

    if logTiming
        print("Flatten:               ")
        @time flat_model = flatten(modified_model)
    else
        flat_model = flatten(modified_model)
    end
    traverseAndSubstituteAllInstances(:(), modified_model, flat_model, true)

    if PrintFlattened
        loglnModia("\nFLATTENED MODEL")
        showInstance(flat_model)
        loglnModia()
    end

    if elaborate
        solved_model = elaborateModel(flat_model)
        if solved_model == nothing
            return nothing
        end
    else
        solved_model = flat_model
    end

    #=
    #  JSON.json(t::This) = "this."

    if PrintJSONsolved
       jsonSolved = JSON.json(solved_model.equations[1])
        @show jsonSolved
    end
    =#
    @static if VERSION < v"0.7.0-DEV.2005"
        disableSimulation = false
    else
        disableSimulation = false
    end

    if !BasicStructuralTransform.newStateSelection && !disableSimulation ## Disable simulation for the moment
        loglnModia("\nSIMULATION")
        # @show incidenceMatrix
        println("model = ", typeof(solved_model.variables))
        res = Dict{AbstractString,Any}
        println("key = $key")
        k = String(key)
        t = vec(collect(Float64, t))

        sim_der(a) = simulate_ida_der(solved_model, var, a, k, t, if useIncidenceMatrix; incidenceMatrix else nothing end, log=logSimulation, relTol=relTol, hev=hev)

        g = gradient(Params(t)) do
         sim_der(var_val)
       end

          if logTiming
                print("Code generation and simulation:         ")
                # @show solved_model t useIncidenceMatrix logSimulation
                @time res = g(t)#simulate_ida(solved_model, t, if useIncidenceMatrix; incidenceMatrix else nothing end, log=logSimulation, relTol=relTol, hev=hev)
            else
                s = g(t)#sim_der'(var_val)
                println(s)
                res = s#simulate_ida(solved_model, t, if useIncidenceMatrix; incidenceMatrix else nothing end, log=logSimulation, relTol=relTol, hev=hev)
            end

    else
        res = Dict{Symbol,AbstractArray{T,1} where T}()
    end

    setDerAsFunction(true)

    if fileStdOut
        #  close(outWrite)
        output = readavailable(outRead)
        close(outRead)
        redirect_stdout(originalSTDOUT)
        print(String(output))
    end

    if logTiming
        println(@sprintf("Total time: %0.3f", (time_ns() - start) * 1E-9), " seconds")
    end
    closeLogModia()
    println("something")
    println("noRes = $noResult")
    println("res = $res")
    #test_expr!("@test", ex, kws...)
    #orig_ex = Expr(:inert, ex)
    #result = get_test_result(ex, __source__)
    #:(do_test($result, $orig_ex))
    if res != noResult
        return res
    else
        throw(DomainError("res = $res"))
    end
    return res
end

function simulateModel(model, t; options...)
    return simulateModelWithOptions(model, t, options=options)
end

function simulate(model, stopTime; startTime=0, options...)
    nSteps = 1000
    @static if VERSION < v"0.7.0-DEV.2005"
        t = linspace(startTime, stopTime, nSteps)
    else
        t = range(startTime, stop=stopTime, length=nSteps)
    end

    return simulateModelWithOptions(model, t, options=options)
end

function simulate_der(model, stopTime, var, var_val, key; startTime=0, options...)
    nSteps = 1000
    @static if VERSION < v"0.7.0-DEV.2005"
        t = linspace(startTime, stopTime, nSteps)
    else
        t = range(startTime, stop=stopTime, length=nSteps)
    end

    return simulateModelWithOptionsDer(model, t, var, var_val, key, options=options)
end

function print_rgb(r, g, b, t)
    println("\e[1m\e[38;2;$r;$g;$b;249m", t)
end

@static if VERSION < v"0.7.0-DEV.2005"
    printstyled(s; bold=false, color=:black) = print_with_color(color, s, bold=bold)
end
"""
    function checkSimulation(mod, stopTime, observer, finalSolution; options...)
Simulates model mod until stopTime and checks that the final value of the observer variable is approximately finalSolution. The following options are available:

* `logTranslation`: determines if logging of the translation is performed
* ...
"""
function checkSimulation(mod, stopTime, observer="", finalSolution=0.0; startTime=0, options...)
    nSteps = 1000
    @static if VERSION < v"0.7.0-DEV.2005"
        t = linspace(startTime, stopTime, nSteps)
    else
        t = range(startTime, stop=stopTime, length=nSteps)
    end
    res = nothing

    try
        res = simulateModelWithOptions(mod, t, options=options)
    catch err
        st = stacktrace(catch_backtrace())
        closeLogModia()

        setTestStatus(false)
        println()
        println("\n----------------------\n")
        println()
            printstyled("Simulation FAILED:", bold=true, color=:red); println()
        if isa(err, ErrorException)
            printstyled(err.msg, bold=true, color=:red); println()
        elseif isa(err, UndefVarError)
            printstyled(err, bold=true, color=:red); println()
            ModiaLogging.increaseLogCategory(:(UndefinedSymbol))
        else
            printstyled(err, bold=true, color=:red); println()
        end
        println()
        println("\n----------------------\n")
        println()
        println("Stack trace: ------------------------------------------------")
        for l in st
            println(l)
        end
        println("End stack trace: --------------- ----------------------------")
        @test false
        return noResult
    end

    final = nothing
    if res != nothing && res != noResult && observer != "" && haskey(res, observer)
        obs = res[observer]
        if length(obs) > 0
            final = obs[end]
        end
    end

    if final != nothing
        println("final $observer = $final")
        ok = finalSolution == nothing || typeof(final) == Float64 && Base.isapprox(final, finalSolution, rtol=1.0E-3) || final == finalSolution
        @test ok

        if !ok
            setTestStatus(false)
            println("final solution $observer = $finalSolution")
            printstyled("Simulation NOT OK", bold=true, color=:red); println()
            println()
#= Replaced PyPlot solution with ModiaMath.plot(..) call
            figure()
            title("Simulation NOT OK in " * string(mod.name))
            plot(res["time"], res[observer])
            grid(true)
            xlabel("time [s]")
            legend([observer],  loc="center right")
=#
            ModiaMath.plot(res, observer, heading="Simulation NOT OK in " * string(mod.name))
            false
        else
            ModiaLogging.increaseLogCategory(:(CalculatedResult))
            setTestStatus(true)
            printstyled("Simulation OK", bold=true, color=:green); println()
            println()
            true
        end

    else
        setTestStatus(true)
        printstyled("Simulation RAN", bold=true, color=:green); println()
        println()
        true
        @test true
    end
    res
end

# ----------------------------

"""
Experimental code for multi-mode handling with impulses.
"""
function simulateMultiModeModel(model, t0, t1; n=1000, m=100, options...)
    opt = Dict(options)
    ModiaLogging.setDefaultLogName(string(model.name))
    ModiaLogging.setOptions(opt)
    StructuralTransform.setOptions(opt)
    BasicStructuralTransform.setOptions(opt)
    Instiantiation.setOptions(opt)
    Execution.setOptions(opt)
    if length(keys(opt)) >= 0
        println("Option(s) not found: ", keys(opt))
    end

    println("\nSimulating model: ", model.name)
    loglnModia("\nSimulating model: ", model.name)

    dt = (t1 - t0) / m
    @static if VERSION < v"0.7.0-DEV.2005"
        times1 = linspace(t0, t0 + dt, 100)
    else
        times1 = range(t0, stop=t0 + dt, length=100)
    end

    model1 = flatten(instantiate(model, t0, Dict()))
    #  model1 = elaborateModel(flatten(instantiate(model, t0, Dict())))
    println("\nSub-simulation 1; startTime=$t0, interval=$dt")
    res1 = simulate_ida(model1, times1, if false; nothing else nothing end, log=false)
    tn = t0 + dt

    for step in 2:m
        println("\nSub-simulation $step; startTime=$tn, interval=$dt")

        # extract final values

        # finals = [name => values[end] for (name,values) in res1]
        finals = Dict(name => if length(values) > 0; values[end] else false end for (name, values) in res1)

        @static if VERSION < v"0.7.0-DEV.2005"
            times2 = linspace(tn, tn + dt, 100)
        else
            times2 = range(tn, stop=tn + dt, length=100)
        end
        model2 = flatten(instantiate(model, tn, Dict()))

        # carry over final values
        vars = vars_of(model2)

        for (name, final) in finals
            name = Symbol(name)
            if haskey(vars, name)
                var = vars[name]
                if isa(var, Variable)
                    var.start = final
#                   @show name final
                end
            end
        end

        res2 = simulate_ida(model2, times2, if false; nothing else nothing end, log=false)
        for (name, values) in res1
            res1[name] = vcat(res1[name], res2[name])
        end
        tn += dt
    end

    return res1
end

end


text:
module OneDimensionalNodes

export spectralderivative, spectralinterpolation, legendregauss,
    legendregausslobatto

include("operators.jl")
include("quadratures.jl")

end # module


text:
# This function solves for multivariable roots of a function g
using LinearAlgebra
using ForwardDiff

function NewtonRoot_mult(g, init_guess)
"""
    #### Fields

    - 'g': a function we want to find root of
    - 'init_guess': initial guess

    #### Returns

    - 'xn1' : a root of g
"""

    xn1 = copy(init_guess);
    gprime = x -> ForwardDiff.jacobian(g, x)
    xn = zeros(length(init_guess));

    tol = 1e-10;
    maxiter = 1000;
    iter = 0;
    distance = 10;

    while distance >= tol && iter < maxiter
        iter = iter + 1;
        xn_temp = xn1 - inv(gprime(xn1))*g(xn1);
        for i in 1:length(xn)
            xn[i] = maximum([xn_temp[i], 1e-10])
        end
        distance = norm(xn - xn1);
        xn1 = copy(xn);
    end

    return xn1
end

## Checks
# f1(x) = log.(x)
# NewtonRoot(f1, [10.0,10.0])