text:
module CustomSaltTest

using Hashids
using Test

const TESTDATA = [
    [0],
    [1, 2, 3],
    [10000000000, 0, 0, 0, 999999999999999],
    [827364827638428762424, 0x410498019283098083],
    [340282366920938463463374607431768211456]
]

function test_with_custom_salt(salt::String)
    config = Hashids.configure(salt)
    for numbers in TESTDATA
        @test decode(config, encode(config, numbers)) == numbers
    end
end

@testset "Custom Salt" begin

@testset "\"\"" begin
    test_with_custom_salt("")
end  # of "\"\""

@testset "\"   \"" begin
    test_with_custom_salt("   ")
end  # of "\"   \""

@testset "long salt" begin
    test_with_custom_salt(raw"abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890`~!@#$%^&*()-_=+\\|'\";:/?.>,<{[}]")
end  # of "long salt"

@testset "weird salt" begin
    test_with_custom_salt(raw"`~!@#$%^&*()-_=+\\|'\";:/?.>,<{[}]")
end  # of "weird salt"

@testset "ultra weird salt" begin
    test_with_custom_salt("🤺👩🏿🦳🛁👩🏻🦷🤦♂🐁☝🏼✍🏾👉🏽🇸🇰❤🍭")
end  # of "ultra weird salt"

end  # of "Custom Salt"

end  # module


text:
using Test, GHOST

pats = [ GitHubPersonalAccessToken(split(pat, ':')...) for pat in filter!(!isempty, split(ENV["GH_PAT"], '\n')) ]
@testset "Basic Pipeline" begin
    setup(pats = pats)
    licenses()
    setup_parallel(2)
    spdxs = execute(GHOST.PARALLELENABLER.conn,
                    "SELECT spdx FROM $(GHOST.PARALLELENABLER.schema).licenses where spdx = ANY('{EUPL-1.2,UPL-1.0}'::text[]) ORDER BY spdx;",
                    not_null = true) |>
        (obj -> getproperty.(obj, :spdx))
    for spdx in spdxs
        queries(spdx)
    end
    data = execute(GHOST.PARALLELENABLER.conn,
                   String(read(joinpath(pkgdir(GHOST), "src",  "assets", "sql", "queries_batches.sql"))),
                   not_null = true) |>
        DataFrame |>
        (df -> groupby(df, [:queries, :query_group]));
    low = minimum(sum(subdf.queries) for subdf in data)
    data = data[findfirst(subdf -> sum(subdf.queries) == low, data)]
    find_repos(data)
    data = execute(GHOST.PARALLELENABLER.conn,
                   String(read(joinpath(pkgdir(GHOST), "src", "assets", "sql", "branches_min_max.sql"))) |>
                       (obj -> replace(obj, "schema" => GHOST.PARALLELENABLER.schema)) |>
                       (obj -> replace(obj, "min_lim" => 1)) |>
                       (obj -> replace(obj, "max_lim" => 2)),
                   not_null = true) |>
        (obj -> getproperty.(obj, :branch))
    query_commits(data, 3)
    chk = execute(GHOST.PARALLELENABLER.conn,
                  "SELECT count(*) > 0 success FROM $(GHOST.PARALLELENABLER.schema).queries WHERE done;",
                  not_null = true) |>
        (obj -> getproperty.(obj, :success)) |>
        only
    @test chk
end
# Documentation
@testset "Documentation" begin
    using Documenter, GHOST

    prefix = ".."

    DocMeta.setdocmeta!(GHOST,
                       :DocTestSetup,
                       :(using GHOST, Documenter;
                         ENV["COLUMNS"] = 120;
                         ENV["LINES"] = 30;),
                       recursive = true)
    # doctest(GHOST, fix = true)
    makedocs(sitename = "GHOST",
             format = Documenter.HTML(assets = [joinpath("assets", "custom.css")]),
             modules = [GHOST],
             pages = [
                 "Introduction" => "index.md",
                 "Manual" => "manual.md",
                 "API" => "api.md",
                 ],
             source = joinpath(prefix, "docs", "src"),
             build = joinpath(prefix, "docs", "build"),
             )
    @test true
end

ext: jl


text:
"""
Structural and symbolic processing.

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

Examples of use can be found in TestSymbolic.jl

"""
module Symbolic

export removeBlock, removeQuoteNode, makeDerVar, append, prepend, Incidence, findIncidence!, linearFactor, solveEquation,
    isLinear, getCoefficients, substitute, removeUnits, resetEventCounters, getEventCounters, substituteForEvents

using Base.Meta: isexpr
using OrderedCollections
using ModiaBase.Simplify
using Unitful
using Measurements
using MonteCarloMeasurements

"""
    e = removeBlock(ex)

Remove :block with LineNumberNode from expressions

* `ex`: Expression or array of expressions
* `return `e`: ex with :block removed
"""
removeBlock(ex) = ex
removeBlock(arr::Array{Any,1}) = [removeBlock(a) for a in arr]
removeBlock(arr::Array{Expr,1}) = [removeBlock(a) for a in arr]
removeBlock(q::QuoteNode) = q
    # if typeof(q.value) == Symbol; q.value else q end # begin print("QuoteNode: "); dump(q); q; end # q.value
removeBlock(d::OrderedDict) = OrderedDict{Symbol, Any}([(k,removeBlock(v)) for (k,v) in d])

function removeBlock(ex::Expr)
    if isexpr(ex, :block) && length(ex.args) == 2 && typeof(ex.args[1]) == LineNumberNode
        ex.args[2]
    elseif isexpr(ex, :block) && length(ex.args) == 3 && typeof(ex.args[1]) == LineNumberNode && typeof(ex.args[2]) == LineNumberNode
        ex.args[3]
    else
        Expr(ex.head, [removeBlock(arg) for arg in ex.args]...)
    end
end

removeQuoteNode(ex) = ex
removeQuoteNode(arr::Array{Any,1}) = [removeQuoteNode(a) for a in arr]
removeQuoteNode(arr::Array{Expr,1}) = [removeQuoteNode(a) for a in arr]
removeQuoteNode(q::QuoteNode) = q.value
    # if typeof(q.value) == Symbol; q.value else q end # begin print("QuoteNode: "); dump(q); q; end # q.value
removeQuoteNode(d::OrderedDict) = OrderedDict{Symbol, Any}([(k,removeQuoteNode(v)) for (k,v) in d])

function removeQuoteNode(ex::Expr)
    Expr(ex.head, [removeQuoteNode(arg) for arg in ex.args]...)
end

removeUnits(ex) = if typeof(ex) <: Unitful.Quantity; @show ex; ustrip(ex) else ex end
removeUnits(ex::Expr) = if ex.head == :macrocall && ex.args[1] == Symbol("@u_str"); 1 else Expr(ex.head, [removeUnits(arg) for arg in ex.args]...) end

append(ex, suffix) = if ex == nothing; suffix else Expr(:., ex, QuoteNode(suffix)) end

prepend(ex, prefix) = ex
#prepend(ex, prefix::Array{Any,1}) = if length(prefix) == 0; ex else prepend(prepend(ex, prefix[end]), prefix[1:end-1]) end
#prepend(ex, prefix::Array{Symbol,1}) = if length(prefix) == 0; ex else prepend(prepend(ex, prefix[end]), prefix[1:end-1]) end
#prepend(ex::Symbol, prefix::Array{Symbol,1}) = if length(prefix) == 0; ex else prepend(prepend(ex, prefix[end]), prefix[1:end-1]) end
#prepend(ex::QuoteNode, prefix::Symbol) = Expr(:., prefix, QuoteNode(ex))
prepend(ex::Symbol, prefix) = if prefix == nothing; ex elseif ex in [:time, :instantiatedModel, :_leq_mode, :_x]; ex else Expr(:., prefix, QuoteNode(ex)) end
prepend(ex::Symbol, prefix::Nothing) = ex
prepend(arr::Array{Expr,1}, prefix) = [prepend(a, prefix) for a in arr]
#prepend(dict::OrderedCollections.OrderedDict{Symbol,Expr}, prefix) = OrderedDict([prepend(k, prefix) => prepend(k, prefix) for (k,v) in dict])

function prepend(ex::Expr, prefix)
    if typeof(prefix) == Array{Any,1} && length(prefix) == 0
        ex
    elseif ex.head == :. && ex.args[1] == :up
        ex.args[2].value
    elseif ex.head in [:call, :kw]
        if false #ex.head == :call && ex.args[1] == :der
            e = Symbol(ex)
            :($prefix.$e)
        else
            Expr(ex.head, ex.args[1], [prepend(arg, prefix) for arg in ex.args[2:end]]...)
        end
    elseif ex.head == :macrocall
        if length(ex.args) >= 1 && ex.args[1] == Symbol("@u_str")
            # Do not expand  units, such as u"s", because otherwise logCode=true results in wrong code, because show(u"N") is N and not u"N".
            ex
        else
            eval(ex)
        end
    else
        Expr(ex.head, [prepend(arg, prefix) for arg in ex.args]...)
    end
end

nCrossingFunctions = 0
nAfter = 0
nClocks = 0
nSamples = 0
previousVars = []
preVars = []
holdVars = []

function resetEventCounters()
    global nCrossingFunctions
    global nAfter
    global nClocks
    global nSamples
	global previousVars
    global preVars
    global holdVars
    nCrossingFunctions = 0
    nAfter = 0
    nClocks = 0
    nSamples = 0
	previousVars = []
    preVars = []
    holdVars = []
end

function getEventCounters()
    global nCrossingFunctions
    global nAfter
    global nClocks
    global nSamples
	global previousVars
    global preVars
    global holdVars
    return (nCrossingFunctions, nAfter, nClocks, nSamples, previousVars, preVars, holdVars)
end

substituteForEvents(ex) = ex

function substituteForEvents(ex::Expr)
    global nCrossingFunctions
    global nAfter
    global nClocks
    global nSamples
	global previousVar
    global preVars
    global holdVars
    if ex.head in [:call, :kw]
        if ex.head == :call && ex.args[1] == :positive
            nCrossingFunctions += 1
            :(positive(instantiatedModel, $nCrossingFunctions, ustrip($(substituteForEvents(ex.args[2]))), $(string(substituteForEvents(ex.args[2]))), _leq_mode))
        elseif ex.head == :call && ex.args[1] == :Clock
            @assert 2<=length(ex.args)<=3 "The Clock function takes one or two arguments: $ex"
             nClocks += 1
             if length(ex.args) == 2
                :(Clock(ustrip($(substituteForEvents(ex.args[2]))), instantiatedModel, $nClocks))
             else
                :(Clock(ustrip($(substituteForEvents(ex.args[2]))), ustrip($(substituteForEvents(ex.args[3]))), instantiatedModel, $nClocks))
             end
        elseif ex.head == :call && ex.args[1] == :sample
            nSamples += 1
            :(sample($(substituteForEvents(ex.args[2])), $(substituteForEvents(ex.args[3])), instantiatedModel, $nSamples))
        elseif ex.head == :call && ex.args[1] == :pre
            if length(ex.args) == 2
                push!(preVars, ex.args[2])
                nPre = length(preVars)
                :(pre(instantiatedModel, $nPre))
            else
                error("The pre function takes one arguments: $ex")
            end
        elseif ex.head == :call && ex.args[1] == :previous
            if length(ex.args) == 3
                push!(previousVars, ex.args[2])
                nPrevious = length(previousVars)
                :(previous($(substituteForEvents(ex.args[3])), instantiatedModel, $nPrevious))
            else
                error("The previous function presently takes two arguments: $ex")
            end
        elseif ex.head == :call && ex.args[1] == :hold
            push!(holdVars, ex.args[2])
            nHold = length(holdVars)
            if length(ex.args) == 3
                :(hold($(substituteForEvents(ex.args[2])), $(substituteForEvents(ex.args[3])), instantiatedModel, $nHold))
            else
#                error("The hold function takes two or three arguments, hold(v, clock) or hold(expr, start, clock) : $ex")
                error("The hold function takes two arguments, hold(v, clock): $ex")
            end
        elseif ex.head == :call && ex.args[1] in [:initial, :terminal]
            if length(ex.args) == 1
                :($(ex.args[1])(instantiatedModel))
            else
                error("The $(ex.args[1]) function don't take any arguments: $ex")
            end
        elseif ex.head == :call && ex.args[1] == :after
            # after(instantiatedModel, nr, t, tAsString, leq_mode)
            nAfter += 1
            :(after(instantiatedModel, $nAfter, ustrip($(substituteForEvents(ex.args[2]))), $(string(substituteForEvents(ex.args[2]))), _leq_mode))
        else
            Expr(ex.head, ex.args[1], [substituteForEvents(arg) for arg in ex.args[2:end]]...)
        end
    else
        Expr(ex.head, [substituteForEvents(arg) for arg in ex.args]...)
    end
end

Incidence = Union{Symbol, Expr}

"""
    findIncidence!(ex, incidence::Array{Incidence,1})

Traverses an expression and finds incidences of Symbols and der(...)

* `ex`: Expression or array of expressions
* `incidence`: array of incidences. New incidences of `ex` are pushed.
"""
findIncidence!(ex, incidence::Array{Incidence,1}, includeX::Bool=true) = nothing
findIncidence!(s::Symbol, incidence::Array{Incidence,1}, includeX::Bool=true) = begin if ! (s in [:(:), :end]); push!(incidence, s) end end
findIncidence!(arr::Array{Any,1}, incidence::Array{Incidence,1}, includeX::Bool=true) = [findIncidence!(a, incidence, includeX) for a in arr]
findIncidence!(arr::Array{Expr,1}, incidence::Array{Incidence,1}, includeX::Bool=true) = [findIncidence!(a, incidence, includeX) for a in arr]
function findIncidence!(ex::Expr, incidence::Array{Incidence,1}, includeX::Bool=true)
    if ex.head == :macrocall && ex.args[1] == Symbol("@u_str")
        nothing
    elseif isexpr(ex, :function)
        nothing
    elseif ex.head == :call
        if ex.args[1] == :der
            push!(incidence, ex) # der(x)
            if includeX
                push!(incidence, ex.args[2]) # x
            end
        elseif ex.args[1] in [:pre, :previous]
            [findIncidence!(e, incidence, includeX) for e in ex.args[3:end]] # skip operator/function name and first argument
		else
            [findIncidence!(e, incidence, includeX) for e in ex.args[2:end]] # skip operator/function name
        end
    elseif ex.head == :.
        push!(incidence, ex)
#		if ex.args[2].value != :all
#			push!(incidence, ex.args[1])
#		end
    elseif ex.head == :generator
        vars = [v.args[1] for v in ex.args[2:end]]
        incid = Incidence[]
        [findIncidence!(e, incid, includeX) for e in ex.args]
        unique!(incid)
        setdiff!(incid, vars)
        push!(incidence, incid...)
    else
        # For example: =, vect, hcat, block, ref
        [findIncidence!(e, incidence, includeX) for e in ex.args]
    end
    nothing
end

"""
    (rest, factor, linear) = linearFactor(ex, x)

Finds the linear `factor` and `rest` if `ex` is `linear` with regards to `x` (ex == factor*x + rest)

* `ex`: Expression
* `return (rest, factor, linear)`:
"""
linearFactor(ex, x) = (ex, 0, true)
linearFactor(ex::Symbol, x::Incidence) =  if ex == x; (0, 1, true) else (ex, 0, true) end
function linearFactor(ex::Expr, x::Incidence)
    if ex.head == :block
        linearFactor(ex.args[1], x)
    elseif ex.head == :macrocall && ex.args[1] == Symbol("@u_str")
        (ex, 0, true)
    elseif isexpr(ex, :call) && ex.args[1] == :der
        if ex == x; (0, 1, true) else (ex, 0, true) end
    elseif isexpr(ex, :call) && ex.args[1] in [:positive, :previous]
        (ex, 0, true)
    elseif isexpr(ex, :call)
        func = ex.args[1]
        if func in [:_DUPLICATEEQUATION, :_DUPLICATIMPLICITDEPENDENCY, :implicitDependency]
            return (ex, 0, false)
        end
        arguments = ex.args[2:end]
        factored = [linearFactor(a, x) for a in arguments]
        rests = [f[1] for f in factored]
        factors = [f[2] for f in factored]
        linears = [f[3] for f in factored]
        if func == :+
            rest = foldl(add, rests)
            factor = foldl(add, factors)
            (rest, factor, all(linears))
        elseif func == :-
            if length(rests) == 1
                rest = sub(0, rests[1])
            else
                rest = foldl(sub, rests)
            end
            if length(factors) == 1
                factor = sub(0, factors[1])
            else
                factor = foldl(sub, factors)
            end
            (rest, factor, all(linears))
        elseif func == :*
            if length(arguments) > 2
                linearFactor(foldl(mult, arguments), x)
            else
                # (r1 + f1*x)*(r2 + f2*x) = (r1*r2 + r1*f2*x + f1*r2*x + ...)
                rest = foldl(mult, rests)
                factor = 0
                if factors[1] == 0
                    factor = mult(rests[1], factors[2])
                    (rest, factor, all(linears))
                elseif length(factors) == 2 && factors[2] == 0
                    factor = mult(rests[2], factors[1])
                    (rest, factor, all(linears))
                else
                    (0, 0, false)
                end
            end
        elseif func == :/
            # (r1 + f1*x)/r2 = (r1/r2 + f1/r2*x)
            rest = foldl(divide, rests)
            @assert length(factors) == 2 "Non-binary division is not handled."
            factor = divide(factors[1], rests[2])
            (rest, factor, all(linears) && all(factors[2:end] .== 0))
        elseif func == :\
            # r1 \ (r2 + f2*x) = (r1\r2 + r1\f2*x)
            rest = divide(rests[2], rests[1])
            @assert length(factors) == 2 "Non-binary \\ is not handled."
            factor = divide(factors[2], rests[1])
            (rest, factor, all(linears) && factors[1] == 0)
        else
            (ex, 0, all(linears) && all(factors .== 0))
        end
    elseif ex.head == :.
        if ex == x; (0, 1, true) else (ex, 0, true) end
    elseif isexpr(ex, :if) || isexpr(ex, :elseif)
        cond = linearFactor(ex.args[1], x)
        then = linearFactor(ex.args[2], x)
        els = linearFactor(ex.args[3], x)
        if then[1] == els[1]
            rest = then[1]
        else
            rest = :(if $(cond[1]); $(then[1]) else $(els[1]) end)
        end
        if then[2] == els[2]
            factor = then[2]
        else
            factor = :(if $(cond[1]); $(then[2]) else $(els[2]) end)
        end
        (rest, factor, cond[3] && then[3] && els[3])
    elseif isexpr(ex, :(=))
        LHS = linearFactor(ex.args[1], x)
        RHS = linearFactor(ex.args[2], x)
        rest = sub(LHS[1], RHS[1])
        factor = sub(LHS[2], RHS[2])
        (rest, factor, LHS[3] && RHS[3])
    elseif isexpr(ex, :vect) || isexpr(ex, :vcat) || isexpr(ex, :hcat) || isexpr(ex, :row)
        arguments = ex.args[2:end]
        factored = [linearFactor(a, x) for a in arguments]
        linears = [f[3] for f in factored]
        (ex, 0, all(linears))
    else
#        @warn "Unknown expression type" ex
#        dump(ex)
        incidence = Incidence[]
        findIncidence!(ex, incidence)
        if x in incidence
            (ex, 0, false)
        else
            (ex, 0, true)
        end
    end
end

"""
    (solved_equation, solved) = solveEquation(equ::Expr, x)

Solves `equ` for `x` if `ex` is linear with regards to `x` (equ == factor*x + rest = 0)

* `equ`: Equation (expression = expression)
* `return (solved_equation, solved)`: If the equation is `solved`, solved_equation constains the corresponding Expr
"""
function solveEquation(equ::Expr, x)
    (rest, factor, linear) = linearFactor(equ, x)
    (:($x = $(divide(sub(0, rest), factor))), linear)
end

"""
    (linear, constant) = isLinear(ex::Expr, x)

Tests if `ex` is `linear` with regards to `x` (ex == factor*x + rest) and checks if the factor is `constant`

* `ex`: Expression
* `return (linear, constant)`
"""
function isLinear(equ::Expr, x::Incidence)
    (rest, factor, linear) = linearFactor(equ, x)
    (linear, typeof(factor) != Expr)
end

"""
    (incidence, coefficients, rest, linear) = getCoefficients(ex)

If `ex` is `linear` with regards to all `incidence` (Symbols and der(...)), the `coefficients` and `rest` are returned

* `ex`: Expression
* `return (incidence, coefficients, rest, linear)`
"""

function getCoefficients(ex::Expr)
    incidence = Incidence[]
    findIncidence!(ex, incidence)
    coefficients = []
    rest = ex
    linear = true
    for x in incidence
        crossIncidence = Incidence[]
        findIncidence!(coefficients, crossIncidence)
        if x in crossIncidence
            linear = false
            break
        end
        (rest, factor, linearRest) = linearFactor(rest, x)
        if !linearRest
            linear = false
            break
        end
        push!(coefficients, factor)
    end
    incidence, coefficients, rest, linear
end

substitute(substitutions, ex) = begin nex = get(substitutions, ex, ex); if nex != ex; nex = substitute(substitutions, nex) else nex end; nex end

substitute(substitutions, ex::Vector{Symbol}) = [substitute(substitutions, e) for e in ex]

substitute(substitutions, ex::Vector{Expr}) = [substitute(substitutions, e) for e in ex]

substitute(substitutions, ex::MonteCarloMeasurements.StaticParticles) = ex

substitute(substitutions, ex::Measurements.Measurement{Float64}) = ex

function substitute(substitutions, ex::Expr)
    if isexpr(ex, :quote)
        ex
    elseif ex.head == :.
        nex = get(substitutions, ex, ex)
        if nex != ex
            nex = substitute(substitutions, nex)
        end
        nex
    elseif ex.head == :call && ex.args[1] == :+ && length(ex.args) ==3
        add(substitute(substitutions, ex.args[2]), substitute(substitutions,ex.args[3]))
    elseif ex.head == :call && ex.args[1] == :- && length(ex.args) ==3
        sub(substitute(substitutions, ex.args[2]), substitute(substitutions,ex.args[3]))
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) ==3
        mult(substitute(substitutions, ex.args[2]), substitute(substitutions,ex.args[3]))
    else
        Expr(ex.head, [substitute(substitutions, arg) for arg in ex.args]...)
    end
end

function testSubstitutions()
    ex = substitute(Dict(:a=>:b), :(a+b=0))
    @show ex
    ex = substitute(Dict(:a=>:b), :(a+b+c=0))
    @show ex
    ex = substitute(Dict(:b=>0), :(a+b=0))
    @show ex
    ex = substitute(Dict(:b=>0.0), :(a+b=0))
    @show ex
    ex = substitute(Dict(:b=>0), :(a+b+c=0))
    @show ex
    ex = substitute(Dict(:b=>0.0), :(a+b+c=0))
    @show ex
    ex = substitute(Dict(:b=>0.0), :(a+b*(c+d)=0))
    @show ex

    println("TEST getCoefficients")
    n = 10000
    @show n
    e = Expr(:call, :f, fill(:x, n)...)
#    @show e
    @time incidence, coefficients, rest, linear = getCoefficients(e)
#    @show incidence coefficients rest linear

    e = Expr(:call, :_DUPLICATEEQUATION, fill(:x, n)...)
#    @show e
    @time incidence, coefficients, rest, linear = getCoefficients(e)
#    @show incidence coefficients rest linear
end

# testSubstitutions()

end


text:
 function j2(n,k,m)
    p, i, seq=[0:n-1], 0, Int[]
    while length(p)>m
        i=(i+k-1)%length(p)
        push!(seq,splice!(p,i+1))
    end
    prt_array(x)=replace(chomp(string(x)),"\n",", ")
    @sprintf("Prisoner killing order: %s.\nSurvivors: %s",prt_array(seq),"["*prt_array(p)*"]")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("ghtype_tests.jl")
include("event_tests.jl")
include("read_only_api_tests.jl")


text:
######################################################################
# se2.jl
#
######################################################################

const enum = enumerate

export ErgodicManagerSE2

"""
`ErgodicManagerSE2(L::Float64, K::Int, bins::Int)`

`ErgodicManagerSE2(example_name::String; K::Int=5, bins::Int=100)`

Valid `example_name` entries are:

* "single gaussian"
* "double gaussian"
"""
struct ErgodicManagerSE2{PL} <: ErgodicManager
    domain::Domain

    M::UnitRange{Int}               # Fourier coefficient indices
    N::UnitRange{Int}
    P::PL

    phi::Array{Float64,3}				# spatial distribution
    phik::Array{Complex{Float64},3}		# spatial Fourier coefficients

    # user should never have to interact with these
    Lambda::Array{Float64,3}			# constants
    cache::Array{Complex{Float64},6}    # stores values of F_mnp(x,y,z)
end

function ErgodicManagerSE2(d::Domain, K::Int=5; M::Int=K, N::Int=K, P=K)
    domain = deepcopy(d)

    M = 0:K
    N = 0:K
    if isa(P, Real)
        P = 1:P
    end

    phi = zeros(x_cells(d), y_cells(d), z_cells(d))
    phik = zeros(Complex{Float64}, length(M), length(N), length(P))

    Lambda = zeros(length(M), length(N), length(P))
    cache = zeros(Complex{Float64}, z_cells(d), y_cells(d), x_cells(d), length(P), length(N), length(M))

    em = ErgodicManagerSE2(domain, M, N, P, phi, phik, Lambda, cache)

    Lambda!(em)
    tic()
    cache!(em)
    println("SE2 cache time = ", round(toq(),3), " seconds")

    return em
end

function ErgodicManagerSE2(d::Domain, phi::Array{Float64,3}, K::Int=5)
    em = ErgodicManagerSE2(d, K)
    copy!(em.phi, phi)
    decompose!(em)
    return em
end

# Calls to F_mnp are expensive as they rely on the bessel function.
# Therefore, we store calls to this function in an array beforehand.
function naive_cache!(em::ErgodicManagerSE2)
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        for xi = 1:x_cells(em), yi = 1:y_cells(em), zi=1:z_cells(em)
            x = x_min(em) + (xi-0.5)*x_size(em)
            y = y_min(em) + (yi-0.5)*y_size(em)
            z = z_min(em) + (zi-0.5)*z_size(em)
            em.cache[zi,yi,xi,pl,ni,mi] = F_mnp(m,n,p,x,y,z)
        end
    end
end

# more efficient way to cache
# reduces calls to bessel function
function cache!(em::ErgodicManagerSE2)
    i = float(im)
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        inm = i^(n-m)
        for xi = 1:x_cells(em), yi = 1:y_cells(em)
            x = x_min(em) + (xi-0.5)*x_size(em)
            y = y_min(em) + (yi-0.5)*y_size(em)
            r = sqrt(x*x + y*y)
            psi = atan2(y,x)    # atan(y/x)
            bj = besselj(m-n, p*r)
            for zi = 1:z_cells(em)
                z = z_min(em) + (zi-0.5)*z_size(em)
                rar = exp(i*(m*psi + (n-m)*z))
                em.cache[zi,yi,xi,pl,ni,mi] = inm * rar * bj
            end
        end
    end
end

# fills the matrix Lambda_{m,n,p}
function Lambda!(em::ErgodicManagerSE2)
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        den_sqrt = (1.0 + m*m + n*n + p*p)
        em.Lambda[mi, ni, pl] = 1.0 / (den_sqrt * den_sqrt)
    end
end

# Setting and computing phi
# TODO: this is a nightmare and probably doesn't work
function phi!(em::ErgodicManagerSE2, dm::VF, ds::MF)
    # first, generate d
    d = zeros(x_cells(em), y_cells(em), z_cells(em))
    d_sum = 0.0
    for xi = 1:x_cells(em)
		x = x_min(em) + (xi-0.5)*x_size(em)
		println("xi = ", xi)
		for yi = 1:y_cells(em)
			y = y_min(em) + (yi-0.5)*y_size(em)
			for zi = 1:z_cells(em)
				z = z_min(em) + (zi-0.5)*z_size(em)
				d[xi,yi,zi] = my_pdf([x,y,z], dm, ds)
				d_sum += d[xi,yi,zi]
			end
		end
	end
    normalize!(d, em.domain.cell_size)
    em.phi = d
end

######################################################################
# Computing Fourier coefficients
######################################################################
# update the Fourier coefficients based on some distribution
# Here, I assume it is discrete, but I should change this...
# TODO: maybe I should do some bounds checking?
function decompose!(em::ErgodicManagerSE2, phi::Array{Float64,3})
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        em.phik[mi,ni,pl] = phi_mnp(em, mi, ni, pl, phi)
    end
end

# computes phi_mnp by iterating over the state space
# uses the cache to evaluate F_mnp
function phi_mnp(em::ErgodicManagerSE2, mi::Int, ni::Int, pl::Int, phi::Array{Float64,3})
	val = 0.0im
	for xi = 1:x_cells(em), yi = 1:y_cells(em), zi = 1:z_cells(em)
        val += phi[xi,yi,zi] * em.cache[zi, yi, xi, pl, ni, mi]
	end
	return val * em.domain.cell_size
end

function F_mnp(m::Int, n::Int, p::Real, x::Float64, y::Float64, z::Float64)
    # compute psi, r
    r = sqrt(x*x + y*y)
    psi = atan2(y,x)	# atan(y/x)
    i = float(im)
    return i^(n-m) * exp(i*(m*psi + (n-m)z)) * besselj(m-n, p*r)
end

function decompose(em::ErgodicManagerSE2, traj::VVF)
    N = length(traj)-1

    ck = zeros(Complex{Float64}, length(em.M), length(em.N), length(em.P))

    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        fk_sum = 0.0im
        # now loop over time
        for i = 0:N-1
            xi = traj[i + 1]
            fk_sum += F_mnp(m, n, p, xi[1], xi[2], xi[3])
        end
        # TODO: check that this is right
        ck[mi, ni, pl] = fk_sum / N
    end
    return ck
end

# reconstructs from Fourier coefficients in ck
# This comes from last equation in chapter 10.3 of
#  Engineering Applications of Noncommutative Harmonic Analysis
#
# Returns an array of complex numbers.
# Even if you expect an array of reals, there will be small complex values
# call real(output) if you want just the real parts
function reconstruct2(em::ErgodicManagerSE2, ck::Array{Complex{Float64},3})

    vals = zeros(Complex{Float64}, x_cells(em), y_cells(em), z_cells(em))

    # this order is much more efficient
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        c = ck[mi, ni, pl]
        for xi = 1:x_cells(em), yi = 1:y_cells(em), zi=1:z_cells(em)
            f = conj(em.cache[zi, yi, xi, pl, ni, mi])
            vals[xi,yi,zi] += c * f * p
        end
	end
	return vals
end

# This is a test... multiplying coefficients when M\neq0,N\neq0 by 2
# Doing this because negatives end up being conjugates
function reconstruct(em::ErgodicManagerSE2, ck::Array{Complex{Float64},3})

    vals = zeros(Complex{Float64}, x_cells(em), y_cells(em), z_cells(em))

    # this order is much better
    for (mi,m) in enum(em.M), (ni,n) in enum(em.N), (pl,p) in enum(em.P)
        c = ck[mi, ni, pl]
        for xi = 1:x_cells(em), yi = 1:y_cells(em), zi=1:z_cells(em)
            f = conj(em.cache[zi, yi, xi, pl, ni, mi])
            vals[xi,yi,zi] += 2*c * f * p
        end
	end
    for (pl,p) in enum(em.P)
        c = ck[1, 1, pl]
        for xi = 1:x_cells(em), yi = 1:y_cells(em), zi=1:z_cells(em)
            f = conj(em.cache[zi, yi, xi, pl, 1, 1])
            vals[xi,yi,zi] -= c * f * p
        end
    end
	return vals
end


text:
module TestPack

greet() = println("Hello")
greetagain() = printlin("Hello for the second time.")

end

max_stars_repo_name: moesphere/TestPack

max_issues_repo_name: moesphere/TestPack

max_forks_repo_name: moesphere/TestPack


text:
using SatelliteDynamics, LinearAlgebra

function kron_delta(a::Real, b::Real)
    return a == b ? 1 : 0
end

"""
Internal helper function to aid in denormalization of gravity field coefficients.
Provides method co compute the factorial ratio: (n-m)!/(n+m)!
Arguments:
- `n::Real`: Gravity model degree, n.
- `m::Real`: Gravity model order, m.
Returns:
- `p::Real`: Factorial product
"""
function facprod(n::Int, m::Int)
    p = 1.0

    for i in (n-m+1):(n+m+1) # Is this correct? Should it be (n-m+1):(n+m)
        p = p/i
    end

    return p
end
"""
Compute the gravitational acceleration at a model given a spherical harmonic
gravity field model.
Arguments:
- `r::Array{<:Real, 1}`: Position of the point in the body (field) fixed frame. [m]
- `coef::Array{<:Real, 2}`: Gravity coefficients stored in dense matrix form. C_nm terns are stored along rows indexed from C_00 in coef[1, 1] to C_nm in coef[n+1, m+1]. S_nm terms are stored along matrix columns with S_n0 stored in coef[0, n+1]
- `n_max::Integer`: Maximum degree coefficient to use in expansion
- `m_max::Integer`: Maximum order coefficient to use in the expansion. Must be less than the degree.
- `r_ref::Real`: Reference distance of the gravity field.
- `GM::Real`: Gravitational constant of central body
- `normralized::Bool`: Whether the input gravity field coefficients are normalized coefficients (Default: true)
Returns:
- `a::Array{<:Real, 1}`: Acceleration in X, Y, and Z inertial directions [m/s^2]
"""
function spherical_harmonic_gravity2(r::Array{<:Real, 1}, coef::Array{<:Real, 2}, n_max::Integer, m_max::Integer, r_ref::Real, GM::Real; normalized::Bool=true)
    # Intermediate computations
    r_sqr = dot(r, r)
    rho   = r_ref^2/r_sqr
    x0    = r_ref * r[1] / r_sqr
    y0    = r_ref * r[2] / r_sqr
    z0    = r_ref * r[3] / r_sqr

    # Initialize Intermediary Matrices
    V = zeros(eltype(r), n_max+2, n_max+2)
    W = zeros(eltype(r), n_max+2, n_max+2)

    # Calculate zonal terms V(n, 0). Set W(n,0)=0.0
    V[0+1, 0+1] = r_ref /sqrt(r_sqr)
    W[0+1, 0+1] = 0.0

    V[1+1, 0+1] = z0 * V[0+1, 0+1]
    W[1+1, 0+1] = 0.0

    for n in 2:(n_max+1)
        V[n+1, 0+1] = ((2*n-1)*z0*V[n+1-1, 0+1] - (n-1)*rho*V[n+1-2,0+1])/n
        W[n+1, 0+1] = 0.0
    end

    # Calculate tesseral and sectoral terms
    for m in 1:m_max+1
        # Calculate V(m,m) to V(n_max+1,m)
        V[m+1, m+1] = (2*m-1)*(x0*V[m+1-1, m+1-1] - y0*W[m+1-1, m+1-1])
        W[m+1, m+1] = (2*m-1)*(x0*W[m+1-1, m+1-1] + y0*V[m+1-1, m+1-1])

        if m <= m_max
            V[m+1+1, m+1] = (2*m+1)*z0*V[m+1,m+1]
            W[m+1+1, m+1] = (2*m+1)*z0*W[m+1,m+1]
        end

        for n in (m+2):(n_max+1)
            V[n+1,m+1] = ((2*n-1)*z0*V[n+1-1,m+1]-(n+m-1)*rho*V[n+1-2,m+1])/(n-m)
            W[n+1,m+1] = ((2*n-1)*z0*W[n+1-1,m+1]-(n+m-1)*rho*W[n+1-2,m+1])/(n-m)
        end
    end

    # Calculate accelerations
    ax = 0.0
    ay = 0.0
    az = 0.0

    for m in 0:m_max
        for n in m:n_max
            C = 0.0
            S = 0.0
            if m == 0
                # Denormalize Coeeficients
                if normalized
                    N = sqrt(2*n+1)
                    C = N*coef[n+1, 0+1]
                else
                    C = coef[n+1, 0+1]
                end

                ax -= C*V[n+1+1, 1+1]
                ay -= C*W[n+1+1, 1+1]
                az -= (n+1)*C*V[n+1+1, 0+1]

            else
                if normalized
                    N = sqrt((2 - kron_delta(0,m))*(2*n+1)*facprod(n,m))
                    C = N*coef[n+1,   m+1]
                    S = N*coef[m+1-1, n+1]
                else
                    C = coef[n+1,   m+1]
                    S = coef[m+1-1, n+1]
                end

                fac =  0.5 * (n-m+1)*(n-m+2)
                ax += +0.5 * (-C * V[n+1+1, m+1+1] - S * W[n+1+1, m+1+1])
                      +fac * (+C * V[n+1+1, m+1-1] + S * W[n+1+1, m+1-1])
                ay += +0.5 * (-C * W[n+1+1, m+1+1] + S * V[n+1+1, m+1+1])
                      +fac * (-C * W[n+1+1, m+1-1] + S * V[n+1+1, m+1-1])
                az += (n-m+1)*(-C*V[n+1+1, m+1] - S*W[n+1+1, m+1])
            end
        end
    end

    a = (GM / (r_ref^2)) * [ax, ay, az]

    return a
end

# export accel_gravity
"""
Computes the accleration caused by Earth gravity as modeled by a spherical
harmonic gravity field.
Arguments:
- `r_sat::Array{<:Real, 1}`: Satellite position in the inertial frame [m]
- `R_eci_ecef::Array{<:Real, 2}`: Rotation matrix transforming a vector from the inertial to body-fixed reference frames.
- `n_max::Integer`: Maximum degree coefficient to use in expansion
- `m_max::Integer`: Maximum order coefficient to use in the expansion. Must be less than the degree.

Return:
- `a::Array{<:Real, 1}`: Gravitational acceleration in X, Y, and Z inertial directions [m/s^2]
References:
1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.56-68.
"""
function accel_gravity2(x, R_eci_ecef, n_max::Int=5, m_max::Int=5)

    # Gravitational parameters of primary body
    # GM    = GRAVITY_MODEL.GM
    # R_ref = GRAVITY_MODEL.R

    # Body-fixed position
    r_bf = R_eci_ecef * x[1:3]

    # Compute spherical harmonic acceleration
    a_ecef = spherical_harmonic_gravity2(r_bf, GRAVITY_MODEL.data, n_max, m_max,
     GRAVITY_MODEL.R, GRAVITY_MODEL.GM, normalized=GRAVITY_MODEL.normalized)

    # Inertial acceleration
    a_eci = R_eci_ecef' * a_ecef

    # Finished
    return a_eci
end

function accel_point_mass2(x, gm_body::Real=GM_EARTH)
    # Restrict inputs to position only. Considered in body frame
    r  = x[1:3]

    # Acceleration
    a = -gm_body * r/norm(r)^3

    return a
end

function myRz(θ)
    return [cos(θ) sin(θ) 0;
           -sin(θ) cos(θ) 0;
            0      0      1]
end
function spherical_harmonic_gravity_wrapper(r_eci,R_eci_ecef)

    # R_eci_ecef = rECItoECEF(epc)

    # spherical harmonic gravity (ECI) - point mass gravity (ECI)
    return accel_gravity2(r_eci,R_eci_ecef) - accel_point_mass2(r_eci)
end

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
"""

    struct RegressionSetting
        formula::FormulaTerm
        data::DataFrame
    end

    Immutable data structure for a regression setting.

# Arguments
- `formula::FormulaTerm`: A formula object describes the linear regression model.
- `data::DataFrame`: DataFrame object holds the data

# Notes
    Implemented methods in this packages accepts linear models as RegressionSetting objects.
    This objects holds the model formula and the data used in regression estimations.

# Examples
```julia-repl
julia> setting = RegressionSetting(@formula(calls ~ year), phones)
RegressionSetting(calls ~ year, 24×2 DataFrame
│ Row │ year  │ calls   │
│     │ Int64 │ Float64 │
├─────┼───────┼─────────┤
│ 1   │ 50    │ 4.4     │
│ 2   │ 51    │ 4.7     │
│ 3   │ 52    │ 4.7     │
│ 4   │ 53    │ 5.9     │
│ 5   │ 54    │ 6.6     │
│ 6   │ 55    │ 7.3     │
│ 7   │ 56    │ 8.1     │
│ 8   │ 57    │ 8.8     │
│ 9   │ 58    │ 10.6    │
│ 10  │ 59    │ 12.0    │
⋮
│ 14  │ 63    │ 21.2    │
│ 15  │ 64    │ 119.0   │
│ 16  │ 65    │ 124.0   │
│ 17  │ 66    │ 142.0   │
│ 18  │ 67    │ 159.0   │
│ 19  │ 68    │ 182.0   │
│ 20  │ 69    │ 212.0   │
│ 21  │ 70    │ 43.0    │
│ 22  │ 71    │ 24.0    │
│ 23  │ 72    │ 27.0    │
│ 24  │ 73    │ 29.0    │)
```
"""
struct RegressionSetting
    formula::FormulaTerm
    data::DataFrame
end

"""

    createRegressionSetting(formula, data)

    Create a regression setting for a given formula and data

# Arguments
- `formula::FormulaTerm`: A formula object describes the linear regression model.
- `data::DataFrame`: DataFrame object holds the data

# Notes
    Implemented methods in this packages accepts linear models as RegressionSetting objects.
    This objects holds the model formula and the data used in regression estimations.
    createRegressionSetting is a helper function for creating RegressionSetting objects.

# Examples
```julia-repl
julia> setting = createRegressionSetting(@formula(calls ~ year), phones)
RegressionSetting(calls ~ year, 24×2 DataFrame
│ Row │ year  │ calls   │
│     │ Int64 │ Float64 │
├─────┼───────┼─────────┤
│ 1   │ 50    │ 4.4     │
│ 2   │ 51    │ 4.7     │
│ 3   │ 52    │ 4.7     │
│ 4   │ 53    │ 5.9     │
│ 5   │ 54    │ 6.6     │
│ 6   │ 55    │ 7.3     │
│ 7   │ 56    │ 8.1     │
│ 8   │ 57    │ 8.8     │
│ 9   │ 58    │ 10.6    │
│ 10  │ 59    │ 12.0    │
⋮
│ 14  │ 63    │ 21.2    │
│ 15  │ 64    │ 119.0   │
│ 16  │ 65    │ 124.0   │
│ 17  │ 66    │ 142.0   │
│ 18  │ 67    │ 159.0   │
│ 19  │ 68    │ 182.0   │
│ 20  │ 69    │ 212.0   │
│ 21  │ 70    │ 43.0    │
│ 22  │ 71    │ 24.0    │
│ 23  │ 72    │ 27.0    │
│ 24  │ 73    │ 29.0    │)
```
"""
function createRegressionSetting(formula::FormulaTerm, data::DataFrame)::RegressionSetting
    return RegressionSetting(formula, data)
end

"""

    designMatrix(setting)

    Return matrix of independent variables including the variable (ones) of the constanst term for a given regression setting.

# Arguments
- `setting::RegressionSetting`: A regression setting object.

# Notes
    Design matrix is a matrix which holds values of independent variables on its columns for a given linear regression model.

# Examples
```julia-repl
julia> setting = createRegressionSetting(@formula(calls ~ year), phones);
julia> designMatrix(setting)
24×2 Array{Float64,2}:
 1.0  50.0
 1.0  51.0
 1.0  52.0
 1.0  53.0
 1.0  54.0
 1.0  55.0
 1.0  56.0
 1.0  57.0
 1.0  58.0
 1.0  59.0
 1.0  60.0
 1.0  61.0
 1.0  62.0
 1.0  63.0
 1.0  64.0
 1.0  65.0
 1.0  66.0
 1.0  67.0
 1.0  68.0
 1.0  69.0
 1.0  70.0
 1.0  71.0
 1.0  72.0
 1.0  73.0
```
"""
function designMatrix(setting::RegressionSetting)::Array{Float64,2}
    mf = ModelFrame(setting.formula, setting.data)
    mm = ModelMatrix(mf)
    return convert(Array{Float64, 2}, mm.m)
end

"""

    responseVector(setting)

    Return vector of dependent variable of a given regression setting.

# Arguments
- `setting::RegressionSetting`: A regression setting object.

# Examples
```julia-repl
julia> setting = createRegressionSetting(@formula(calls ~ year), phones);
julia> responseVector(setting)
24-element Array{Float64,1}:
   4.4
   4.7
   4.7
   5.9
   6.6
   7.3
   8.1
   8.8
  10.6
  12.0
  13.5
  14.9
  16.1
  21.2
 119.0
 124.0
 142.0
 159.0
 182.0
 212.0
  43.0
  24.0
  27.0
  29.0
```
"""
function responseVector(setting::RegressionSetting)::Array{Float64,1}
    mf = ModelFrame(setting.formula, setting.data)
return convert(Array{Float64,1}, setting.data[:,mf.f.lhs.sym])
end

"""

    applyColumns(f, data)

    Apply function f to each columns of data.

# Arguments
- `f::Function`: A function that takes a one dimensional array as argument.
- `data::DataFrame`: A DataFrame object.

"""
function applyColumns(f::Function, data::DataFrame)
    return [f(col) for col = eachcol(data)]
end

"""

    applyColumns(f, data)

    Apply function f to each columns of data.

# Arguments
- `f::Function`: A function that takes a one dimensional array as argument.
- `data::Matrix`: A Matrix object.
"""
function applyColumns(f::Function, data::Matrix)
    return [f(col) for col = eachcol(data)]
end

"""

    find_minimum_nonzero(arr)

    Return minimum of numbers greater than zero.

# Arguments
- `arr::Array{Float64, 1}`: A function that takes a one dimensional array as argument.

# Example
```julia-repl
julia> find_minimum_nonzero([0.0, 0.0, 5.0, 1.0])
1.0
```
"""
function find_minimum_nonzero(arr::Array{Float64,1})
    return minimum(filter(x -> x > 0, arr))
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
module steadystate

using ..states, ..operators, ..operators_dense, ..superoperators, ..bases
using ..timeevolution
using Arpack, LinearAlgebra

"""
    steadystate.master(H, J; <keyword arguments>)

Calculate steady state using long time master equation evolution.

# Arguments
* `H`: Arbitrary operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
        operator type.
* `rho0=dm(basisstate(b))`: Initial density operator. If not given the
        ``|0⟩⟨0|`` state in respect to the choosen basis is used.
* `tol=1e-3`: Tracedistance used as termination criterion.
* `hmin=1e-7`: Minimal time step used in the time evolution.
* `rates=ones(N)`: Vector or matrix specifying the coefficients for the
        jump operators.
* `Jdagger=dagger.(Jdagger)`: Vector containing the hermitian conjugates of the
        jump operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given this function `fout(t, rho)` is called every time an
        output should be displayed. To limit copying to a minimum the given
        density operator `rho` is further used and therefore must not be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(H::AbstractOperator{B,B}, J::Vector;
                rho0::DenseOperator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                hmin=1e-7, tol=1e-3,
                rates::Union{Vector{Float64}, Matrix{Float64}, Nothing}=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where B<:Basis
    t,u = timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger,
                        hmin=hmin, hmax=Inf,
                        display_initialvalue=false,
                        display_finalvalue=false,
                        display_intermediatesteps=true,
                        fout=fout,
                        steady_state = true,
                        tol = tol, kwargs...)
end

"""
    steadystate.liouvillianspectrum(L)
    steadystate.liouvillianspectrum(H, J)

Calculate eigenspectrum of the Liouvillian matrix `L`. The eigenvalues and -states are
sorted according to the absolute value of the eigenvalues.

# Keyword arguments:
* `nev = min(10, length(L.basis_r[1])*length(L.basis_r[2]))`: Number of eigenvalues.
* `which = :LR`: Find eigenvalues with largest real part. Keyword for `eigs`
    function (ineffective for DenseSuperOperator).
* `kwargs...`:  Keyword arguments for the Julia `eigen` or `eigens` function.
"""
function liouvillianspectrum(L::DenseSuperOperator; nev::Int = min(10, length(L.basis_r[1])*length(L.basis_r[2])), which::Symbol = :LR, kwargs...)
    d, v = eigen(L.data; kwargs...)
    indices = sortperm(abs.(d))[1:nev]
    ops = DenseOperator[]
    for i in indices
        data = reshape(v[:,i], length(L.basis_r[1]), length(L.basis_r[2]))
        op = DenseOperator(L.basis_r[1], L.basis_r[2], data)
        push!(ops, op)
    end
    return d[indices], ops
end

function liouvillianspectrum(L::SparseSuperOperator; nev::Int = min(10, length(L.basis_r[1])*length(L.basis_r[2])), which::Symbol = :LR, kwargs...)
    d, v, nconv, niter, nmult, resid = try
        eigs(L.data; nev = nev, which = which, kwargs...)
    catch err
        if isa(err, SingularException) || isa(err, ARPACKException)
            error("Arpack's eigs() algorithm failed; try using DenseOperators or change nev.")
        else
            rethrow(err)
        end
    end
    indices = sortperm(abs.(d))[1:nev]
    ops = DenseOperator[]
    for i in indices
        data = reshape(v[:,i], length(L.basis_r[1]), length(L.basis_r[2]))
        op = DenseOperator(L.basis_r[1], L.basis_r[2], data)
        push!(ops, op)
    end
    return d[indices], ops
end

liouvillianspectrum(H::AbstractOperator{B,B}, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) where B<:Basis = liouvillianspectrum(liouvillian(H, J; rates=rates); kwargs...)

"""
    steadystate.eigenvector(L)
    steadystate.eigenvector(H, J)

Find steady state by calculating the eigenstate with eigenvalue 0 of the Liouvillian matrix `L`, if it exists.

# Keyword arguments:
* `tol = 1e-9`: Check `abs(eigenvalue) < tol` to determine zero eigenvalue.
* `nev = 2`: Number of calculated eigenvalues. If `nev > 1` it is checked if there
    is only one eigenvalue with real part 0. No checks for `nev = 1`: use if
    faster or for avoiding convergence errors of `eigs`. Changing `nev` thus only
    makes sense when using SparseSuperOperator.
* `which = :LR`: Find eigenvalues with largest real part. Keyword for `eigs` function (ineffective for DenseSuperOperator).
* `kwargs...`:  Keyword arguments for the Julia `eigen` or `eigs` function.
"""
function eigenvector(L::SuperOperator; tol::Real = 1e-9, nev::Int = 2, which::Symbol = :LR, kwargs...)
    d, ops = liouvillianspectrum(L; nev = nev, which = which, kwargs...)
    if abs(d[1]) > tol
        error("Eigenvalue with smallest absolute value is not zero.")
    end
    if nev > 1
        if abs(real(d[2])) < tol
            @warn("Several eigenvalues with real part 0 detected; use steadystate.liouvillianspectrum to find out more.")
        end
    end
    return ops[1]/tr(ops[1])
end

eigenvector(H::AbstractOperator, J::Vector; rates::Union{Vector{Float64}, Matrix{Float64}}=ones(Float64, length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...)

end # module


text:
using ConstrainedDynamics
using ConstrainedControl

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1

p2 = [0.0;0.0;length1 / 2] # joint connection point

# Initial orientation
ϕ1 = 0
q1 = UnitQuaternion(RotX(ϕ1))

# Links
origin = Origin{Float64}()
link1 = Box(width, depth, length1, length1, color = RGBA(1., 1., 0.))

# Constraints
joint_between_origin_and_link1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=p2))

links = [link1]
constraints = [joint_between_origin_and_link1]

mech = Mechanism(origin, links, constraints)
setPosition!(origin,link1,p2 = p2,Δq = q1)

pid = PID(mech, joint_between_origin_and_link1.id, pi/2, P = 10., I = 10., D = 5.)
@test true


text:
"""
    function angle_between(v1, v2)

Calculates the angle between two vectors
"""
function angle_between(v1, v2)
    _v1 = normalize(v1)
    _v2 = normalize(v2)
    angle = acos(min(dot(_v1, _v2), 1))
end

"""
    function azimuth(d::Direction)

Calculate the azimuth angle for a given direction.
"""
azimuth(d) = atan(d[2], d[1])

"""
    function zenith(d::Direction)

Calculate the zenith angle for a given direction.
"""
zenith(d) = acos(-d[3]/norm(d))

"""
    function cartesian(ϕ, θ; r=1.0)

Calculate the cartesian coordinates for given `ϕ`, `θ` and radius `r`.
"""
function cartesian(ϕ, θ; r=1.0)
    return Position(r*sin(θ) * cos(ϕ), r*sin(θ) * sin(ϕ), cos(θ))
end

"""
    function pld3(p1, p2, d2)

Calculate the distance between a point (p1) and a line (given by p2 and d2).
"""
function pld3(p1, p2, d2)
    norm(cross(d2, (p2 - p1))) / norm(d2)
end

"""
    function lld3(P, u, Q, v)

Calculate the distance between two lines.
"""
function lld3(P, u, Q, v)
    R = Q - P
    n = cross(u, v)
    return norm(R⋅n) / norm(n)
end

"""
    function lld3(t₁::T, t₂::T) where T::Track

Calculate the distance between two tracks.
"""
function lld3(t₁::Track, t₂::Track)
    return lld3(t₁.pos, t₁.dir, t₂.pos, t₂.dir)
end

"""
    function project(P, t::Track)

Projects a point to a track.
"""
function project(P, t::Track)
    A = t.pos
    B = A + t.dir
    project(P, A, B)
end

"""
    function project(P, A, B)

Project P onto a line spanned by A and B.
"""
function project(P, A, B)
    Position(A + ((P-A)⋅(B-A))/((B-A)⋅(B-A)) * (B-A))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "URCU"
version = v"0.13.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://git.lttng.org/userspace-rcu.git", "c19d74445b6ed2e6840038afc2f974e3e098334b")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/userspace-rcu/
./bootstrap
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(!Sys.iswindows, supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("liburcu-memb", :liburcu_memb),
    LibraryProduct("liburcu-mb", :liburcu_mb),
    LibraryProduct("liburcu-qsbr", :liburcu_qsbr),
    LibraryProduct("liburcu-signal", :liburcu_signal),
    LibraryProduct("liburcu-cds", :liburcu_cds),
    LibraryProduct("liburcu", :liburcu),
    LibraryProduct("liburcu-common", :liburcu_common),
    LibraryProduct("liburcu-bp", :liburcu_bp)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@safetestset "Fullstack app" begin

  testdir = pwd()
  using Pkg

  @safetestset "Create and run a full stack app with resources" begin
    using Genie

    content = "Test OK!"

    workdir = Base.Filesystem.mktempdir()
    cd(workdir)

    Genie.newapp("fullstack_test", fullstack = true, testmode = true)

    Genie.Generator.newcontroller("Foo", pluralize = false)
    @test isfile(joinpath("app", "resources", "foo", "FooController.jl")) == true

    mkpath(joinpath("app", "resources", "foo", "views"))
    @test isdir(joinpath("app", "resources", "foo", "views")) == true

    open(joinpath("app", "resources", "foo", "views", "foo.jl.html"), "w") do io
      write(io, content)
    end
    @test isfile(joinpath("app", "resources", "foo", "views", "foo.jl.html")) == true

    Genie.Router.route("/test") do
      Genie.Renderer.Html.html(:foo, :foo)
    end

    r = Genie.Requests.HTTP.request("GET", "http://localhost:8000/test")

    @test occursin(content, String(r.body)) == true

    down()
  end;

  cd(testdir)
  Pkg.activate(".")

end;


text:
include("../src/BlueMesh7.jl")
using .BlueMesh7: generate_graph, initialize_mesh, start, plotgrid
using Random

minutes = 20

Random.seed!(7)
positions = generate_positions()

node_roles = zeros(Int, length(positions))
mesh = initialize_mesh(positions, node_roles)

Random.seed!(7)
produced, received = start(mesh; minutes)
println("None relays had $(round(received / produced; digits=3)) success of delivery")

node_roles = ones(Int, length(positions))
mesh = initialize_mesh(positions, node_roles)

Random.seed!(7)
produced, received = start(mesh; minutes)
println("All relays had $(round(received / produced; digits=3)) success of delivery")

node_roles = rand(0:1, length(positions))
mesh = initialize_mesh(positions, node_roles)

Random.seed!(7)
produced, received = start(mesh; minutes)
println("Half relays had $(round(received / produced; digits=3)) success of delivery")


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
using StatsPlots, StatsBase
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")

################################################################################
#####
##### Lotka-Volterra example: get data
#####
include(EXAMPLE_FOLDERNAME*"example_lotkaVolterra.jl")

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

# Get data
sol = solve_lv_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

plot(sol)
    plot!(ts, U_true[1,:])
    title!("Data")

## Also get baseline true/ideal cases
# Uncontrolled
dat_raw = Array(solve_lv_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)

#####
##### Loop over additive noise
#####
noise_vals = 0.0:0.05:0.3
num_models = 20
all_err = zeros(length(noise_vals), num_models)
all_err_deriv = zeros(length(noise_vals), num_models)
all_err_deriv_subsample = zeros(length(noise_vals), num_models)
all_naive_err = zeros(length(noise_vals), num_models)
all_naive_err_deriv = zeros(length(noise_vals), num_models)

all_i = zeros(length(noise_vals), num_models)

# Make model template
variable_names = ["x", "y"]
optimizer = slstNumber((1, 1));
library = Dict("cross_terms"=>[2], "constant"=>nothing)
model_template = sindyModel(convert_string2function(library),
    variable_names, optimizer)

for (i,σ) in enumerate(noise_vals)
    for j in 1:num_models
        noisy_grad = numerical_grad .+ σ.*randn(size(dat))
        # Initialize
        global this_model = sra_stateful_object(ts, tspan, dat, u0, noisy_grad)
        # Reset parameters because we are NOT using control
        this_model.parameters = get_sra_defaults(false)
        prams = this_model.parameters
        prams.model_template = model_template
        prams.sindy_terms_list = Iterators.product(1:3, 1:3)
        prams.initial_subsampling = true

        fit_first_model(this_model, 30);
        print_true_equations(this_truth)
        print_current_equations(this_model)
        all_naive_err[i, j] = calc_coefficient_error(this_model, this_truth)
        all_naive_err_deriv[i, j] = rss_sindy_derivs(
                this_model.sindy_model, dat, numerical_grad) ./ length(dat)

        #################################################################
        ### Iterate
        is_improved, max_iter = true, 10
        while is_improved
            calculate_subsampled_ind(this_model);
            all_models, is_improved = fit_model(this_model);
            is_improved && print_current_equations(this_model)
            this_model.i > 10 && break
        end
        println("Finished with $(this_model.i) iterations")

        all_err[i, j] = calc_coefficient_error(this_model, this_truth)
        all_err_deriv[i, j] = rss_sindy_derivs(
                this_model.sindy_model, dat, numerical_grad)
        ind = this_model.subsample_ind
        all_err_deriv_subsample[i, j] = rss_sindy_derivs(
                this_model.sindy_model, dat[:,ind], numerical_grad[:,ind]) ./
                (length(ind)*size(dat,1))
        all_i[i, j] = this_model.i
        println("Final error in coefficients = $(all_err[i, j])")
        println("")
    end
    println("Finished noise level $σ")
    println("=====================================")
end

coef_norm = sum(core_dyn_true.A.^2)
vec_naive, std_naive = mean_and_std(all_naive_err./coef_norm, 2)
vec_err, std_err = mean_and_std(all_err./coef_norm, 2)
# Error in derivatives: only look at where there's no control
vec_deriv, std_deriv = mean_and_std(all_err_deriv_subsample, 2)
vec_naive_deriv, std_naive_deriv = mean_and_std(all_naive_err_deriv, 2)

#####
##### Save Lotka-Volterra data
#####
this_dat_name = "dat_library_with_noise_lv_"

save_for_plotting(noise_vals,
                vec_naive, std_naive, vec_err, std_err,
                vec_deriv, std_deriv, vec_naive_deriv, std_naive_deriv,
                this_dat_name)


text:
# * Neon

NeonHF = PP"""! [He] 2s2 2p6
!  Q=8., MEFIT, HF, Ref 22.
ECP,Ne,2,3,0;
1; 2,1.000000,0.00000000;
2; 2,32.029855,112.52822460; 2,12.294136,28.44589604;
2; 2,21.525406,-10.93982872; 2,13.091176,3.39426049;
1; 2,0.850385,-0.16409678;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

NeonWB = PP"""! [He] 2s2 2p6
!  Q=8., MEFIT, WB, Ref 22.
ECP,Ne,2,3,0;
1; 2,1.000000,0.00000000;
2; 2,31.860162,112.52543566; 2,12.362219,28.30083454;
2; 2,21.508034,-11.12658543; 2,12.910447,3.38754919;
1; 2,0.850385,-0.18408921;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

# * Argon

ArgonHF = PP"""! [Ne] 3s2 3p6
!  Q=8., MEFIT, HF, Ref 22.
ECP,Ar,10,4,0;
1; 2,1.000000,0.00000000;
2; 2,10.291970,68.69327818; 2,3.947937,24.48814203;
2; 2,5.388907,27.59672144; 2,2.704463,4.13654627;
2; 2,8.084705,-8.14025335; 2,4.018469,-1.66241708;
1; 2,5.295301,-3.40922653;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

ArgonWB = PP"""! [Ne] 3s2 3p6
!  Q=8., MEFIT, WB, Ref 22.
ECP,Ar,10,4,0;
1; 2,1.000000,0.00000000;
2; 2,10.261721,68.66778801; 2,3.952725,24.04276629;
2; 2,5.392714,27.73076331; 2,2.699967,4.04545904;
2; 2,8.086235,-8.13747696; 2,4.016632,-1.66452808;
1; 2,5.208459,-3.40009845;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

# * Calcium

CalciumDF = PP"""! [Ne] 3s2 3p6 4s2
!  Q=10., MEFIT, MCDHF+Breit, Ref 41.
ECP,Ca,10,4,3;
1; 2,1.,0.;
2; 2,10.556058402,138.832728745; 2,5.293154907,16.754824599;
4; 2,12.714571445,27.661607611; 2,12.878516057,55.297862576; 2,4.147923622,4.210974307; 2,3.841287457,7.627482586;
4; 2,4.703264136,-0.514324716; 2,4.707975898,-0.765184799; 2,13.696151450,-6.490621584; 2,13.698165108,-9.735555558;
2; 2,14.820019137,-14.000148237; 2,22.000803239,-61.616066332;
4; 2,12.714571445,-55.323215221; 2,12.878516057,55.297862576; 2,4.147923622,-8.421948613; 2,3.841287457,7.627482586;
4; 2,4.703264136,0.514324716; 2,4.707975898,-0.510123200; 2,13.696151450,6.490621584; 2,13.698165108,-6.490370372;
2; 2,14.820019137,9.333432158; 2,22.000803239,-30.808033166;

! References:
! [41] I.S. Lim, H. Stoll, P. Schwerdtfeger, J. Chem. Phys. 124, 034107 (2006).
"""sso

# * Krypton

KryptonDF = PP"""! [Ne] 3s2 3p6 3d10 4s2 4p6
!  Q=26., MEFIT, MCDHF+Breit, Ref 36.
ECP,Kr,10,4,3;
1; 2,1.000000,0.000000;
3; 2,70.034410,49.953835; 2,31.895971,369.978238; 2,7.353728,10.054544;
4; 2,47.265183,99.101833; 2,46.664268,198.232524; 2,23.008133,28.204670; 2,22.100979,56.516792;
6; 2,50.768165,-18.588877; 2,50.782228,-27.875398; 2,15.479497,-0.293411; 2,15.559383,-0.469399; 2,2.877154,0.068881; 2,1.991196,0.132270;
2; 2,15.437567,-1.223389; 2,22.055742,-2.991233;
4; 2,47.265183,-198.203666; 2,46.664268,198.232524; 2,23.008133,-56.409339; 2,22.100979,56.516792;
6; 2,50.768165,18.588877; 2,50.782228,-18.583599; 2,15.479497,0.293411; 2,15.559383,-0.312933; 2,2.877154,-0.068881; 2,1.991196,0.088180;
2; 2,15.437567,0.815593; 2,22.055742,-1.495617;
! References:
! [36] K.A. Peterson, D. Figgen, E. Goll, H. Stoll, M. Dolg, J. Chem. Phys. 119, 11113 (2003).
"""sso

KryptonDF2c = PP"""! [Ne] 3s2 3p6 3d10 4s2 4p6
!  Q=26., MEFIT, MCDHF+Breit, Ref 36.
ECP,Kr,10,4,3;
3; 2,70.034410,49.953835; 2,31.895971,369.978238; 2,7.353728,10.054544;
2; 2,46.664268,297.348786; 2,22.100979,84.775187;
3; 2,50.782228,-46.458997; 2,15.559383,-0.782332; 2,1.991196,0.220450;
1; 2,22.055742,-5.234658;
2; 2,47.265183,297.305500; 2,23.008133,84.614009;
3; 2,50.768165,-46.472192; 2,15.479497,-0.733528; 2,2.877154,0.172203;
1; 2,15.437567,-2.854575;
! References:
! [36] K.A. Peterson, D. Figgen, E. Goll, H. Stoll, M. Dolg, J. Chem. Phys. 119, 11113 (2003)."""r

# * Xenon

XenonHF = PP"""! [Ar] 3d10c 4s2c 4p6c 4d10c 5s2 5p6
!  Q=8., MEFIT, HF, Ref 22.
ECP,Xe,46,4,0;
1; 2,1.000000,0.00000000;
2; 2,3.815600,122.76231371; 2,1.878604,8.20135456;
2; 2,3.003078,68.75731963; 2,1.283819,3.64084871;
2; 2,2.027610,23.08319830; 2,0.830435,3.17298823;
2; 2,6.246157,-48.30201724; 2,1.564998,-6.91572892;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

XenonWB = PP"""! [Ar] 3d10c 4s2c 4p6c 4d10c 5s2 5p6
!  Q=8., MEFIT, WB, Ref 22; CPP: alpha=0.8345;delta=0.97391;ncut=1.
ECP,Xe,46,5,0;
1; 2,1.000000,0.00000000;
2; 2,3.940263,122.76382934; 2,2.277264,8.30885115;
2; 2,3.028373,68.82300437; 2,1.394319,3.64674223;
2; 2,2.122605,23.65207854; 2,0.798669,3.25844113;
2; 2,6.164360,-47.70319876; 2,1.542374,-6.54113991;
1; 2,1.847892,-7.10585060;
! References:
! [22] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102, 8942 (1995)."""

XenonDF = PP"""! [Ar] 3d10c 4s2 4p6 4d10 5s2 5p6
!  Q=26., MEFIT, MCDHF+Breit, Ref 36.
ECP,Xe,28,4,3;
1; 2,1.000000,0.000000;
3; 2,40.005184,49.997962; 2,17.812214,281.013303; 2,9.304150,61.538255;
4; 2,15.701772,67.439142; 2,15.258608,134.874711; 2,9.292184,14.663300; 2,8.559003,29.354730;
6; 2,15.185600,35.436908; 2,14.284500,53.195772; 2,7.121889,9.046232; 2,6.991963,13.223681; 2,0.623946,0.084853; 2,0.647284,0.044155;
4; 2,20.881557,-23.089295; 2,20.783443,-30.074475; 2,5.253389,-0.288227; 2,5.361188,-0.386924;
4; 2,15.701772,-134.878283; 2,15.258608,134.874711; 2,9.292184,-29.326600; 2,8.559003,29.354730;
6; 2,15.185600,-35.436908; 2,14.284500,35.463848; 2,7.121889,-9.046232; 2,6.991963,8.815787; 2,0.623946,-0.084853; 2,0.647284,0.029436;
4; 2,20.881557,15.392863; 2,20.783443,-15.037237; 2,5.253389,0.192151; 2,5.361188,-0.193462;
! References:
! [36] K.A. Peterson, D. Figgen, E. Goll, H. Stoll, M. Dolg, J. Chem. Phys. 119, 11113 (2003)."""sso

XenonDF2c = PP"""! [Ar] 3d10c 4s2 4p6 4d10 5s2 5p6
!  Q=26., MEFIT, MCDHF+Breit, Ref 36.
ECP,Xe,28,4,3;
3; 2,40.005184,49.997962; 2,17.812214,281.013303; 2,9.304150,61.538255;
2; 2,15.258608,202.312067; 2,8.559003,44.032095;
3; 2,14.284500,88.659619; 2,6.991963,22.039468; 2,0.647284,0.073591;
2; 2,20.783443,-52.630331; 2,5.361188,-0.677116;
2; 2,15.701772,202.317425; 2,9.292184,43.989900;
3; 2,15.185600,88.592271; 2,7.121889,22.615580; 2,0.623946,0.212133;
2; 2,20.881557,-53.875022; 2,5.253389,-0.672529;
! References:
! [36] K.A. Peterson, D. Figgen, E. Goll, H. Stoll, M. Dolg, J. Chem. Phys. 119, 11113 (2003)."""r

# * Radon

RadonDF = PP"""! [Kr] 4d10c 4f14c 5s2 5p6 5d10 6s2 6p6
!  Q=26., MEFIT, MCDHF+Breit, Ref 36.
ECP,Rn,60,5,4;
1; 2,1.000000,0.000000;
3; 2,30.151242,49.965551; 2,14.521241,283.070000; 2,8.052038,62.002870;
4; 2,11.009942,71.969119; 2,9.617625,143.860559; 2,7.336008,4.714761; 2,6.406253,9.013065;
4; 2,8.369220,36.368365; 2,8.116975,54.551761; 2,5.353656,9.634487; 2,5.097212,14.387902;
4; 2,6.348571,21.797290; 2,6.295949,28.946805; 2,2.882118,1.447365; 2,2.908048,2.177964;
4; 2,11.015205,-25.228095; 2,10.909853,-30.566121; 2,3.482761,-0.574800; 2,3.600418,-0.779538;
4; 2,11.009942,-143.938238; 2,9.617625,143.860559; 2,7.336008,-9.429523; 2,6.406253,9.013065;
4; 2,8.369220,-36.368365; 2,8.116975,36.367841; 2,5.353656,-9.634487; 2,5.097212,9.591935;
4; 2,6.348571,-14.531526; 2,6.295949,14.473402; 2,2.882118,-0.964910; 2,2.908048,1.088982;
4; 2,11.015205,12.614048; 2,10.909853,-12.226448; 2,3.482761,0.287400; 2,3.600418,-0.311815;
! References:
! [36] K.A. Peterson, D. Figgen, E. Goll, H. Stoll, M. Dolg, J. Chem. Phys. 119, 11113 (2003)."""sso


text:
export Metric, AccuracyMetric, RankingMetric
export measure, count_true_positive

abstract type Metric end

abstract type AccuracyMetric <: Metric end

function measure(metric::AccuracyMetric, truth::AbstractVector, pred::AbstractVector)
    error("measure is not implemented for metric type $(typeof(metric))")
end

abstract type RankingMetric <: Metric end

function measure(metric::RankingMetric, truth::Array{T}, pred::Array{T}, k::Int) where T
    error("measure is not implemented for metric type $(typeof(metric))")
end

function count_true_positive(truth::Array{T}, pred::Array{T}) where T
    sum(in(truth), pred)
end


text:
using Test, Flux, NeuralNetDiffEq
using DiffEqDevTools

# Run a solve on scalars
linear = (u,p,t) -> cos(2pi*t)
tspan = (0.0f0, 1.0f0)
u0 = 0.0f0
prob = ODEProblem(linear, u0 ,tspan)
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE(chain,opt), dt=1/20f0, verbose = true,
            abstol=1e-10, maxiters = 200)

# Run a solve on vectors
linear = (u,p,t) -> [cos(2pi*t)]
tspan = (0.0f0, 1.0f0)
u0 = [0.0f0]
prob = ODEProblem(linear, u0 ,tspan)
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE(chain,opt), dt=1/20f0, abstol=1e-10,
            verbose = true, maxiters=200)

#Example 1
linear = (u,p,t) -> @. t^3 + 2*t + (t^2)*((1+3*(t^2))/(1+t+(t^3))) - u*(t + ((1+3*(t^2))/(1+t+t^3)))
linear_analytic = (u0,p,t) -> [exp(-(t^2)/2)/(1+t+t^3) + t^2]
prob = ODEProblem(ODEFunction(linear,analytic=linear_analytic),[1f0],(0.0f0,1.0f0))
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),verbose = true, dt=1/5f0)
err = sol.errors[:l2]
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),verbose = true, dt=1/20f0)
sol.errors[:l2]/err < 0.5

#=
dts = 1f0 ./ 2f0 .^ (6:-1:2)
sim = test_convergence(dts, prob, NeuralNetDiffEq.NNODE(chain, opt))
@test abs(sim.𝒪est[:l2]) < 0.1
@test minimum(sim.errors[:l2]) < 0.5
=#

#Example 2
linear = (u,p,t) -> -u/5 + exp(-t/5).*cos(t)
linear_analytic = (u0,p,t) ->  exp(-t/5)*(u0 + sin(t))
prob = ODEProblem(ODEFunction(linear,analytic=linear_analytic),0.0f0,(0.0f0,1.0f0))
chain = Flux.Chain(Dense(1,5,σ),Dense(5,1))
opt = Flux.ADAM(0.1, (0.9, 0.95))
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),verbose = true, dt=1/5f0)
err = sol.errors[:l2]
sol  = solve(prob,NeuralNetDiffEq.NNODE(chain,opt),verbose = true, dt=1/20f0)
sol.errors[:l2]/err < 0.5

#=
dts = 1f0 ./ 2f0 .^ (6:-1:2)
sim = test_convergence(dts, prob, NeuralNetDiffEq.NNODE(chain, opt))
@test abs(sim.𝒪est[:l2]) < 0.5
@test minimum(sim.errors[:l2]) < 0.1
=#


text:
abstract type AxesLayout end

struct HasFast{M} <: AxesLayout
    fast::Int
    slows::NTuple{M, Int}
end

struct NoFast{N} <: AxesLayout end
struct NoAxes <: AxesLayout end

NoFast(N::Int) = NoFast{N}()

NoSlow(_N::Int) = HasFast(1, ())  # maybe create a new type?

# Base.ndims(::T) where {T <: AxesLayout} = ndims(T)
# Base.ndims(::Type{HasFast{M}}) where M = M + 1
# Base.ndims(::Type{NoFast{N}}) where N = N
# Base.ndims(::Type{NoAxes}) = 0

# fastaxis(al::HasFast) = al.fast
# fastaxis(::AxesLayout) = nothing
slowaxes(al::HasFast) = al.slows
slowaxes(::NoFast{N}) where N = ntuple(identity, N)
slowaxes(::NoAxes) = ()

slowaxes(::T) where {T} = slowaxes(AxesLayout(T))

#=
AxesLayout(::Nothing, slow::NTuple) = NoFast(slow)
AxesLayout(fast::Int, slow::NTuple) = HasFast(fast, slow)
=#

const ColumnMajor{N} = Union{
    # List taken from `strides`:
    DenseArray{<:Any, N},
    StridedReshapedArray{<:Any, N},
    StridedReinterpretArray{<:Any, N},
}

AxesLayout(::T) where T = AxesLayout(T)
AxesLayout(T::Type) = error("AxesLayout is unknown for: ", T)
AxesLayout(::Type{<:ColumnMajor{N}}) where N =
    HasFast(1, Base.tail(ntuple(identity, N)))
AxesLayout(::Type{<:Union{
    <:Adjoint{<:Any, <:ColumnMajor{1}},
    <:Adjoint{<:Any, <:ColumnMajor{2}},
    <:Transpose{<:Any, <:ColumnMajor{1}},
    <:Transpose{<:Any, <:ColumnMajor{2}},
}}) = HasFast(2, (1,))
AxesLayout(::Type{A}) where {A <: ReferenceableArray} = AxesLayout(parenttype(A))

AxesLayout(::Type{<:AbstractArray{<:Any, N}}) where N = NoFast{N}()
AxesLayout(::Type{<:Ref}) = NoAxes()

AxesLayout(::Type{BC}) where {BC <: Broadcasted} = axeslayoutargs(argstype(BC))
argstype(::Type{<:Broadcasted{<:Any, <:Any, <:Any, Args}}) where Args = Args
axeslayoutargs(::Type{Tuple{}}) = NoAxes()
axeslayoutargs(::Type{Args}) where {Args <: Tuple} =
    bcaxeslayout(AxesLayout(tuple_type_head(Args)),
                 axeslayoutargs(tuple_type_tail(Args)))

# "Computed/virtual" arrays has no slow axes.
AxesLayout(::Type{<:AbstractRange}) = NoSlow(1)  # is it?
AxesLayout(::Type{<:Union{
    CartesianIndices{N},
    LinearIndices{N},
}}) where N =
    NoSlow(N)

# Merge two `AxesLayout` in a way compatible with broadcasting
bcaxeslayout(x::NoFast, ::AxesLayout) = x
bcaxeslayout(::NoAxes, y::AxesLayout) = y
bcaxeslayout(::HasFast, y::NoFast) = y
bcaxeslayout(x::HasFast, ::NoAxes) = x
bcaxeslayout(x::HasFast{N}, y::HasFast{M}) where {N, M} =
    if x.fast !== y.fast
        NoFast{N}()
    else
        slows = foldlargs((), ntuple(identity, max(N, M))...) do slows, i
            if i === x.fast
                slows
            else
                (slows..., i)
            end
        end
        return HasFast(x.fast, slows)
    end

struct Index{x} end
const Indices{N} = NTuple{N, Index}

Index(x::Index) = x
Index(x) = Index{x}()

indexname(::Index{x}) where x = x

struct AccessPattern{A, T <: Indices}
    indexable::A
    indices::T
end

Base.ndims(p::AccessPattern) = ndims(p.indexable)

AccessPattern(pair::Pair) = AccessPattern(pair[1], map(Index, pair[2]))

struct AccessPlan{T <: Indices}
    indices::T
end

"""
    plan(args::AccessPattern...) :: AccessPlan
    plan(pairs::Pair...) :: AccessPlan

Try to find an index that is categorized as "fast" for all `args` with
such index.

# Examples
```jldoctest
julia> using NDReducibles: plan

julia> A1 = ones(0)
       A2 = ones(0, 0)
       A3 = ones(0, 0, 0)
       A4 = ones(0, 0, 0, 0);

julia> plan(
           A1 => (:i,),
           A2 => (:i, :j),
       )
AccessPlan: j → i

julia> plan(
           A1 => (:i,),
           A2 => (:j, :i),
           A2 => (:j, :k),
       )
AccessPlan: i → k → j

julia> plan(
           A1 => (:i,),
           A2' => (:j, :i),
           A2 => (:j, :k),
       )
AccessPlan: j → k → i

julia> plan(
           A1 => (:i,),
           A2' => (:i, :j),
           A2 => (:i, :j),
           A3 => (:k, :i, :j),
           A4 => (:k, :j, :l, :m),
       )
AccessPlan: i → j → l → m → k
```
"""
plan(pairs::Pair...) = plan(map(AccessPattern, pairs)...)

function plan(args::AccessPattern...)
    allindices = foldlargs((), args...) do indices, pattern
        tupleunion(indices, pattern.indices)
    end :: Indices
    fastindices = foldlargs(allindices, args...) do indices, pattern
        slows = map(slowaxes(pattern.indexable)) do a
            pattern.indices[a]
        end
        tuplediff(indices::Indices, slows::Indices)
    end :: Indices
    # TODO: Sort `fastindices` by the number of arrays using it?
    return AccessPlan((tuplediff(allindices, fastindices)..., fastindices...))
end


text:
"""
Run all callbacks
"""
function _run_callbacks(callbacks::C, scenes::Union{EntityQueueRecord{S,D,I}, Vector{Frame{Entity{S,D,I}}}}, actions::Union{Nothing, Vector{Frame{A}}}, roadway::R, models::Dict{I,M}, tick::Int) where {S,D,I,A<:EntityAction,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
    isdone = false
    for callback in callbacks
        isdone |= run_callback(callback, scenes, actions, roadway, models, tick)
    end
    return isdone
end

function simulate!(
    ::Type{A},
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    callbacks::C,
    ) where {S,D,I,A,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
    Base.depwarn(
"`simulate!` using `EntityQueueRecord`s is deprecated since v0.7.10 and may be removed in future versions.
 You should pass a pre-allocated vector of entitites `scenes::Vector{Frame{Entity{S,D,I}}}` to `simulate!`
 or use the convenience function `simulate` without pre-allocation instead.",
        :simulate_rec
    )

    empty!(rec)
    update!(rec, scene)

    # potential early out right off the bat
    if _run_callbacks(callbacks, rec, nothing, roadway, models, 0)
        return rec
    end

    actions = Array{A}(undef, length(scene))
    for tick in 1 : nticks
        get_actions!(actions, scene, roadway, models)
        tick!(scene, roadway, actions, get_timestep(rec))
        update!(rec, scene)
        if _run_callbacks(callbacks, rec, nothing, roadway, models, tick)
            break
        end
    end

    return rec
end
function simulate!(
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    callbacks::C,
    ) where {S,D,I,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}

    return simulate!(Any, rec, scene, roadway, models, nticks, callbacks)
end

## Implementations of useful callbacks

"""
    CollisionCallback

Terminates the simulation once a collision occurs
"""
@with_kw struct CollisionCallback
    mem::CPAMemory=CPAMemory()
end

function run_callback(
    callback::CollisionCallback,
    rec::EntityQueueRecord{S,D,I},
    roadway::R,
    models::Dict{I,M},
    tick::Int
    ) where {S,D,I,R,M<:DriverModel}

    return !is_collision_free(rec[0], callback.mem)
end

function run_callback(
    callback::CollisionCallback,
    scenes::Vector{Frame{E}},
    actions::Union{Nothing, Vector{Frame{A}}},
    roadway::R,
    models::Dict{I,M},
    tick::Int
) where {E<:Entity,A<:EntityAction,R,I,M<:DriverModel}
    return !is_collision_free(scenes[tick], callback.mem)
end


text:
using ResTime
using Test

@testset "ResTime.jl" begin

    # Write your tests here.
end


text:
mutable struct LineNumberingIO
    # the number of the line last read, first line is line 1
    line::Int
    io::IO

    LineNumberingIO(io::IO) = new(0, io)
end

function LineNumberingIO(f, io::IO)
    lnio = LineNumberingIO(io)
    try
        f(lnio)
    finally
        close(lnio)
    end
end

Base.close(in::LineNumberingIO ) = close(in.io)

Base.position(in::LineNumberingIO) = position(in.io)

function Base.readline(in::LineNumberingIO; keep::Bool=false)
    in.line += 1
    readline(in.io; keep=keep)
end


text:
function get_price_data_single(currency::String)

    date = Dates.today()

    filename = "$(currency)_EUR_data_$(date).csv"
    filepath = joinpath("data", filename)

    raw_df = DataFrame()

    # Look for present day's CSV file, if not found, download and save data to a new file
    if isfile(filepath)
        @info "Reading $(currency) price/vol data from CSV file on disk"
        raw_df = CSV.File(filepath) |> DataFrame
    else
        try
            @info "Fetching $(currency) price/vol data from Alpha Vantage"
            raw = AlphaVantage.digital_currency_daily(currency, "EUR", datatype="csv")
            raw_df = raw_to_df(raw)
            CSV.write(filepath, raw_df)
        catch
            @info "Could not fetch data, try again later!"
        end
    end

    # Return processed DataFrame only when raw data has been fetched successfully
    if ~isempty(raw_df)
        df_out_price, df_out_candle = average_price_df(currency, raw_df)
        df_out_vol = vol_df(currency, raw_df)

        return df_out_price, df_out_candle, df_out_vol
    else
        return DataFrame[], DataFrame[], DataFrame[]
    end
end

function get_ratings_data(currency::String)

    date = Dates.today()

    filename = "$(currency)_metrics_data_$(date).csv"
    filepath = joinpath("data", filename)

    metrics_df = DataFrame()

    # Look for present day's CSV file, if not found, download and save data to a new file
    if isfile(filepath)
        @info "Reading $(currency) FCAS data from CSV file on disk"
        metrics_df = CSV.File(filepath) |> DataFrame
    else
        try
            @info "Fetching $(currency) FCAS data from Alpha Vantage"
            rating = AlphaVantage.crypto_rating(currency)
            scores = rating["Crypto Rating (FCAS)"]
            CSV.write(filepath, scores)
            metrics_df = CSV.File(filepath) |> DataFrame
        catch err
            if isa(err, KeyError)
                @info "Could not retrieve data. Something wrong with the API, try again later!"
            else
                @info "This is a new error: $(err)"
            end
        end
    end

    ratings = ["utility", "fcas score", "developer", "market", "fcas rating"]
    index = Array{Int64}(undef,0)

    # Return scores only when metrics data has been fetched successfully
    if ~isempty(metrics_df)

        # Find the row(index) of a string match
        for rating in ratings
            i = findall(occursin.(rating, metrics_df[!,1]))
            push!(index, i[1])
        end

        # Variables should be assigned in the same order as the list of ratings above
        utility_score, fcas_score, dev_score = metrics_df[!,2][index[1]], metrics_df[!,2][index[2]], metrics_df[!,2][index[3]]
        mark_score, fcas_rating = metrics_df[!,2][index[4]], metrics_df[!,2][index[5]]

        return utility_score, fcas_score, dev_score, mark_score, fcas_rating
    else
        return Number[], Number[], Number[], Number[], Number[]
    end

end

################# Functions for CoinGecko API #################

function get_API_response(params::String, url::String = URL)

	CG_request = HTTP.request("GET", url * params; verbose = 0, retries = 2)
	response_text = String(CG_request.body)
	response_dict = JSON.parse(response_text)

	return response_dict
end

function get_coin_id(currency::String)

    date = Dates.today()

    filename = "List_of_all_coins_$(date).csv"
    filepath = joinpath("data", filename)

    df_coins, df_filter, df_filter_1 = [DataFrame() for i = 1:3]

    # Look for present day's CSV file, if not found, download and save data to a new file
    if isfile(filepath)
        @info "Reading list of coins from CSV file on disk"
        df_coins = CSV.File(filepath) |> DataFrame
    else
        try
            @info "Fetching list of coins from CoinGecko"
            coins_dict = get_API_response("/coins/list")
            df_coins = vcat(DataFrame.(coins_dict)...)
            CSV.write(filepath, df_coins)
        catch
            @info "Could not fetch data, try again later!"
        end
    end

    # Return valid coin id only when list of coins is available
    if ~isempty(df_coins)

        # Filter on matching currency
        df_filter = df_coins |> @filter(_.symbol == currency) |> DataFrame

        try
            # For multiple matches, first filter on coin ids and then on names,
            # which do not have "-" in them
            if size(df_filter)[1] > 1

                df_filter_1 = df_filter |>
                            @filter(~occursin("-", _.id)) |> DataFrame

                if isempty(df_filter_1)
                    df_filter_1 = df_filter |>
                            @filter(~occursin("-", _.name)) |> DataFrame
                end

                return df_filter_1[!, :id][1]
            end

            return df_filter[!, :id][1]

        catch err
            if isa(err, BoundsError)
                @info "Could not find an id for the given currency"
            else
                @info "Something went wrong, check this error: $(err)"
            end
        end

    else
        return ""
    end
end

function dict_to_df(data_dict::Dict, df::DataFrame)

    # Collect only the key-value data which is suitable for plotting
    for key in collect(keys(data_dict))
        if ~isnothing(data_dict[key]) && length(data_dict[key]) == 1
            push!(df, [key Float64(data_dict[key])])
        end
    end

    return df
end

function get_dev_comm_data(currency::String)

    coin_id = get_coin_id(currency)

    coin_dict, dev_dict, comm_dict = [Dict() for i = 1:3]

    try
        @info "Fetching coin data from CoinGecko"
        coin_dict = get_API_response("/coins/$(coin_id)")
    catch
        @info "Could not fetch data, try again later!"
    end

    # Get developer data
    try
        dev_dict = coin_dict["developer_data"]
    catch err
        if isa(err, KeyError)
            @info "Could not find developer data!"
        else
            @info "This is a new error: $(err)"
        end
    end

    # Get community data
    try
        comm_dict = coin_dict["community_data"]
    catch err
        if isa(err, KeyError)
            @info "Could not find community data!"
        else
            @info "This is a new error: $(err)"
        end
    end

    # Convert dict to DataFrame
    df_dev, df_comm = [DataFrame(Metric = String[], Value = Float64[]) for i = 1:2]

    if ~isempty(dev_dict)
        df_dev = dict_to_df(dev_dict, df_dev)
    end

    if ~isempty(comm_dict)
        df_comm = dict_to_df(comm_dict, df_comm)
    end

    return df_dev, df_comm
end

function get_list_of_exchanges(num_exchanges::Int64)

    date = Dates.today()

    filename = "List_of_top_$(num_exchanges)_exchanges_$(date).csv"
    filepath = joinpath("data", filename)

    df_all_ex, df_ex = [DataFrame() for i = 1:2]

    # Look for present day's CSV file, if not found, download and save data to a new file
    if isfile(filepath)
        @info "Reading list of exchanges from CSV file on disk"
        df_ex = CSV.File(filepath) |> DataFrame
    else
        try
            @info "Fetching list of exchanges from CoinGecko"
            ex_dict = get_API_response("/exchanges")
            df_all_ex = vcat(DataFrame.(ex_dict)...)

            # Keep only name and id columns
            df_ex = DataFrame(name = df_all_ex[!, :name], id = df_all_ex[!, :id])
            CSV.write(filepath, df_ex)
        catch err
            @info "Could not fetch data, try again later!"
            @info "$(err)"
        end
    end

    # Return filtered list of exchanges only when data is available
    if ~isempty(df_ex)
        df_ex_vol = DataFrame(Name = df_ex[!, :name][1:num_exchanges],
                               ID  = df_ex[!, :id][1:num_exchanges])

        return df_ex_vol
    else
        return DataFrame()
    end
end

function get_exchange_vol_data(currency::String, num_exchanges::Int64)

    coin_id = get_coin_id(currency)

    df_ex_vol = get_list_of_exchanges(num_exchanges)
    allowmissing!(df_ex_vol)

    exchange_coin_vol, exchange_usd_vol = [Union{Missing, Float64}[] for i = 1:2]

    # Extract volume data only when a list of exchanges is available
    if ~isempty(df_ex_vol)

        for exchange in df_ex_vol[!, :ID]

            coin_vol, usd_vol = [Float64[] for i = 1:2]
            coin_vol_tickers_dict = Dict()

            try
                coin_vol_tickers_dict = get_API_response("/exchanges/$(exchange)/tickers?coin_ids=$(coin_id)")
            catch err
                @info "Could not find $(coin_id) volume data on $(exchange)!"
                @info "$(err)"
            end

            if ~isempty(coin_vol_tickers_dict)

                for i = 1:length(coin_vol_tickers_dict["tickers"])
                    push!(coin_vol, coin_vol_tickers_dict["tickers"][i]["volume"])
                    push!(usd_vol, coin_vol_tickers_dict["tickers"][i]["converted_volume"]["usd"])
                end

                push!(exchange_coin_vol, sum(coin_vol))
                push!(exchange_usd_vol, sum(usd_vol))

            else
                push!(exchange_coin_vol, missing)
                push!(exchange_usd_vol, missing)
            end

        end
    end

    insertcols!(df_ex_vol, 3, :Coin_volume => exchange_coin_vol, :USD_volume => exchange_usd_vol)

    return df_ex_vol
end

function get_vol_chart(exchange::String)

    date = Dates.today()

    filename = "$(exchange)_vol_data_$(date).txt"
    filepath = joinpath("data", filename)

    ex_vol_chart = Vector{Any}[]

    # Historical data is fetched and saved for these many days. Takes more time for longer duration,
    # and queries for many exchanges don't even return more points.
    days = 365

    if isfile(filepath)
        @info "Reading $(exchange) vol data from file on disk"
        ex_vol_chart = readdlm(filepath, ';')

        # Drop the dimension and convert to Float64
        ex_vol_chart = Float64.(ex_vol_chart[:])
    else
        try
            # Fetch and save data for 365 days
            ex_vol_chart = get_API_response("/exchanges/$(exchange)/volume_chart?days=$(days)")

            open(filepath, "w") do f
                for i = 1:length(ex_vol_chart)
                    writedlm(f, [ex_vol_chart[i]], ";")
                end
            end

        catch
            @info "Could not find volume data for $(exchange)"
            @info "$(err)"
        end
    end

    if length(ex_vol_chart) > days
        start_index = Int64(length(ex_vol_chart)/2) + 1
        return ex_vol_chart[start_index:end]

    else
        ex_vol = Union{Missing, Float64}[]

        for i = 1:length(ex_vol_chart)
            try
                push!(ex_vol, round(parse(Float64, ex_vol_chart[i][2]); digits = 2))
            catch
                push!(ex_vol, missing)
            end
        end

        return ex_vol
    end
end

function get_overall_vol_data(duration::Int64, num_exchanges::Int64)

    df_ex_list = get_list_of_exchanges(num_exchanges)

    # Check on duration, current maximum is set to 365
    if duration > 365
        duration = 365
    end

    # Create a column with dates
    f_day = Dates.today()
    i_day = f_day - Dates.Day(duration-1)

    time = collect(i_day:Dates.Day(1):f_day)

    df_ex_vol = DataFrame(Time = time)
    allowmissing!(df_ex_vol)

    # Extract volume data only when a list of exchanges is available
    if ~isempty(df_ex_list)

        for exchange in df_ex_list[!, :ID]

            ex_vol = Vector{Any}[]

            try
                ex_vol = get_vol_chart(exchange)
            catch err
                @info "Could not find volume data for $(exchange), will continue to next!"
                @info "$(err)"

                # Skip next part of the code, and continue to next exchange
                continue
            end

            # Check and filter on duration
            if duration > length(ex_vol)
                duration = length(ex_vol)
            end
            ex_vol = ex_vol[end-duration+1:end]

            # Find name of the exchange corresponding to its ID
            df_row = df_ex_list |> @filter(_.ID == exchange) |> DataFrame
            name = df_row[!, :Name][1]

            try
                insertcols!(df_ex_vol, 2, Symbol(name) => ex_vol)
            catch err
                if isa(err, DimensionMismatch)
                    @info "Data is missing for $(exchange) for the requested duration"
                else
                    @info "Something went wrong, check this error: $(err)"
                end
            end

        end

        return df_ex_vol
    else
        return DataFrame()
    end

end


text:
struct MeshGenerationOptions
    distance_fcn
	node_spacing_fcn
	initial_spacing::Float64
	bounding_box::Rectangle
	fixed_points
	max_iteration::Int

    delta_point_tolerance::Float64 	   # = 0.001;
	retriangulation_threshold::Float64 # = 0.1;
	F_scale::Float64 				   # = 1.2;
	time_step_size::Float64 		   # = 0.2;
	geometrical_tolerance::Float64 	   # = 0.001 * initial_spacing;
	gradient_step::Float64 			   # = sqrt(eps)*initial_spacing;
end
function MeshGenerationOptions(
    distance_fcn,
	node_spacing_fcn,
	initial_spacing,
	bounding_box,
	fixed_points;
	max_iteration             = 10,
    delta_point_tolerance 	  = 0.001,
	retriangulation_threshold = 0.1,
	F_scale 				  = 1.2,
	time_step_size 			  = 0.1,
	geometrical_tolerance 	  = 0.001 * initial_spacing,
	gradient_step 			  = 10*sqrt(eps(1.0))*initial_spacing
    )

    return MeshGenerationOptions(
        distance_fcn,
        node_spacing_fcn,
        initial_spacing,
        bounding_box,
        fixed_points,
        max_iteration,
        delta_point_tolerance,
        retriangulation_threshold,
        F_scale,
        time_step_size,
        geometrical_tolerance,
        gradient_step)

end

function create_mesh(options::MeshGenerationOptions; visualize = false)

    # 1. create initial rough grid
    (xs, ys) = create_initial_meshgrid_distribution(options)
    (xs, ys) = shift_even_rows!(xs, ys, options)

    # 2. filter points
    (xs, ys) = remove_points_outside_region!(xs, ys, options)
    (xs, ys) = keep_interior_points!(xs, ys, options)

    # create the point array
    point_array = create_point_array(xs, ys, options)
    n_points = size(point_array,1)

    if visualize
        PyPlot.figure()
        PyPlot.grid()
        PyPlot.plot(xs, ys,
            marker = :., linestyle = :none, markersize = 10)
    end

    # For first iteration
	point_array_old = Inf # .* ones(n_points, 2);

    # out params
    triangle_index = Matrix{Int}(undef, 0, 3)
    delmesh = DelaunayMesh()

    bars = Matrix{Int}(undef, 0, 2)

    # main loop
    iter = 0
    for outer iter = 1:options.max_iteration
        println("At iter: $(iter)")

        # 3. retriangluation
        dp_mat = point_array .- point_array_old
        # dp_vec = map( x -> sum(x.^2), dp_mat )
        dp_vec = sum(dp_mat.^2; dims=2)
        max_dp = maximum(dp_vec)
        # println(max_dp)

        if iter == 1 || max_dp >= options.retriangulation_threshold * options.initial_spacing
            println("Retriangulation at iter $(iter)")
            # Save current positions
            point_array_old = copy(point_array)

            # create delaunay triangulation
            pointtyped_point_array = matrix_to_point_array(point_array)
            delmesh = build_delaunay_mesh!(pointtyped_point_array)
            triangle_index = triangle_indices(delmesh)

            # compute centroids
            p1 = point_array[triangle_index[:, 1], :];
			p2 = point_array[triangle_index[:, 2], :];
			p3 = point_array[triangle_index[:, 3], :];
			triangle_centroids = (p1 .+ p2 .+ p3) ./ 3;

            # Keep interior triangles
			dist_at_triange_centroids = vec(mapslices( options.distance_fcn, triangle_centroids; dims = 2))
			bool_interior_triangle    = dist_at_triange_centroids .< (-options.geometrical_tolerance);
			triangle_index 			  = triangle_index[bool_interior_triangle, :];

            if visualize
                px = point_array[:,1]
                py = point_array[:,2]
                PyPlot.figure()
                PyPlot.grid()

                PyPlot.plot( px, py,
                    marker = :., linestyle = :none, markersize = 10)

                for tt = 1:size(triangle_index,1)
                    t = triangle_index[tt,:]
                    tx = px[ [t[1], t[2], t[3], t[1]]]
                    ty = py[ [t[1], t[2], t[3], t[1]]]
                    PyPlot.plot(tx, ty)
                end
            end

            # 4. Describe each bar by a unique pair of nodes
			# Interior bars duplicated
			bars = [	triangle_index[:, [1, 2]];
						triangle_index[:, [1, 3]];
						triangle_index[:, [2, 3]]
					]

			# Bars as node pairs
            sort!(bars; dims = 2)
			bars = unique(bars; dims = 2);

        end # if

        # 6. Move mesh points based on bar lengths bar_lengths and forces F_bar
		bar_vector 			= point_array[bars[:, 1], :] - point_array[bars[:, 2], :];
		# bar_lengths 		= sqrt.( map( x -> sum(x.^2), bar_vector ) );
        bar_lengths         = vec(sqrt.( sum( bar_vector.^2; dims = 2) ))
		bar_midpoints 		= 0.5 .* ( point_array[bars[:, 1], :] + point_array[bars[:, 2], :] );
		hbars 				= vec(mapslices( options.node_spacing_fcn, bar_midpoints; dims = 2))
		desired_bar_length 	= hbars .* options.F_scale .* sqrt.( sum(bar_lengths.^2) ./ sum(hbars.^2) );

        # Bar forces (scalars)
		F_bar = max.(desired_bar_length - bar_lengths, 0.0)
        # println("F_bar = $(F_bar[1:5])")

		# Bar forces (x, y components)
		Fvec = F_bar .* bar_vector ./ bar_lengths

		n_bars = size(bar_vector, 1)
		F_total = zeros(n_points, 2)

		for bb = 1:n_bars

			start_index = bars[bb,1];
			end_index   = bars[bb,2];

			force_from_start_to_end = Fvec[bb,:]
			F_total[start_index,:]  = F_total[start_index,:] + force_from_start_to_end;
			F_total[end_index,:]    = F_total[end_index,:]   - force_from_start_to_end;

		end

        # Force  =  0 at fixed points
		F_total[1:size(options.fixed_points, 1), :] .= 0.0

		# Update node positions
		point_array = point_array .+ options.time_step_size * F_total;

		# 7. Bring outside points back to the boundary
		distance_at_points = vec(mapslices( options.distance_fcn, point_array; dims = 2))

		# Find points outside (distance_at_points > 0)
		outside = distance_at_points .> 0;

        dist_out = distance_at_points[outside]
		tmp = point_array[outside,:]

        # arg = map( p -> [p[1] + options.gradient_step, p[2]], tmp )
        arg = copy(tmp)
        arg[:,1] .+= options.gradient_step
		dx = vec(mapslices( options.distance_fcn, arg; dims = 2))
		grad_x_of_distfun = ( dx - dist_out ) ./ options.gradient_step

        # arg = map( p -> [p[1], p[2] + options.gradient_step], tmp )
		arg = copy(tmp)
        arg[:,2] .+= options.gradient_step
        dy = vec(mapslices( options.distance_fcn, arg; dims = 2))
		grad_y_of_distfun = ( dy - dist_out ) ./ options.gradient_step;

        # Project back to boundary
		point_array[outside,:] .-= [dist_out .* grad_x_of_distfun     dist_out .* grad_y_of_distfun]

        # 8. Termination criterion: All interior nodes move less than dptol (scaled)
		delta_p_for_interior = options.time_step_size .* F_total[distance_at_points .< (-options.geometrical_tolerance), :];
		# delta_p_for_interior = sqrt.( map( x -> sum(x.^2), delta_p_for_interior ) );
        delta_p_for_interior = vec(sqrt.( sum( delta_p_for_interior.^2; dims = 2) ))
		if maximum(delta_p_for_interior) < options.initial_spacing * options.delta_point_tolerance
			break
		end

    end # for

    # final retriangluation

    # create delaunay triangulation
    pointtyped_point_array = matrix_to_point_array(point_array)
    delmesh = build_delaunay_mesh!(pointtyped_point_array)
    triangle_index = triangle_indices(delmesh)

    # compute centroids
    p1 = point_array[triangle_index[:, 1], :];
    p2 = point_array[triangle_index[:, 2], :];
    p3 = point_array[triangle_index[:, 3], :];
    triangle_centroids = (p1 .+ p2 .+ p3) ./ 3;

    # Keep interior triangles
    dist_at_triange_centroids = vec(mapslices( options.distance_fcn, triangle_centroids; dims = 2))
    bool_interior_triangle    = dist_at_triange_centroids .< (-options.geometrical_tolerance);
    triangle_index 			  = triangle_index[bool_interior_triangle, :];

    # set the valid in geometry indices
    delmesh.triangle_index = triangle_index

    return (point_array, triangle_index, delmesh)

end

function create_initial_meshgrid_distribution(options::MeshGenerationOptions)

    x_min = options.bounding_box.p1.x
    y_min = options.bounding_box.p1.y

    x_max = options.bounding_box.p2.x
    y_max = options.bounding_box.p2.y

    x_min = min(x_min, x_max)
    x_max = max(x_min, x_max)
    y_min = min(y_min, y_max)
    y_max = max(y_min, y_max)

    dx = x_max - x_min
    dy = y_max - y_min
    Nx = Int64( max(2, ceil(dx / options.initial_spacing) ) )
    Ny = Int64( max(2, ceil(dy / (sqrt(3)/2 * options.initial_spacing)) ) )

    xr = range(x_min, x_max, Nx)
    yr = range(y_min, y_max, Ny)

    (xs, ys) = util.mesh_grid(xr, yr)
    return (xs, ys)

end

function shift_even_rows!(xs, ys, options::MeshGenerationOptions)
    xs[:,2:2:end] .+=  options.initial_spacing / 2
    xs = vec(xs)
    ys = vec(ys)
    return (xs, ys)
end

function remove_points_outside_region!(xs, ys, options::MeshGenerationOptions)

    L = length(xs)
    dist = zeros(L)

    for ii = 1:L
        dist[ii] = options.distance_fcn( (xs[ii], ys[ii]) )
    end

    interior = dist .<= options.geometrical_tolerance

    xs = xs[interior]
    ys = ys[interior]
    return (xs, ys)
end

function keep_interior_points!(xs, ys, options::MeshGenerationOptions)

    L = length(xs)

    # Probability to keep point
    node_spacings_at_points = zeros(L)
    for ii = 1:L
        node_spacings_at_points[ii] = options.node_spacing_fcn( (xs[ii], ys[ii]) );
    end
	r0 = 1 ./ node_spacings_at_points.^2;
	max_r0 = maximum(r0)

	# Rejection method
	probability = rand(L);
	keep_point  = probability .< (r0 ./ max_r0);

    xs = xs[keep_point]
    ys = ys[keep_point]

    return (xs, ys)
end

function create_point_array(xs, ys, options::MeshGenerationOptions)

    #=
    Lf = length(options.fixed_points)
    Lx = length(xs)
    L = Lx + Lf
    point_array = Vector{MutablePoint}(undef, L)

    for ii = 1:Lf
        point_array[ii] = options.fixed_points[ii]
    end
    for ii = (Lf+1):L
        point_array[ii] = MutablePoint(xs[ii], ys[ii])
    end
    =#

    point_array = vcat(options.fixed_points, [xs ys])

    return point_array
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia