text:
Base.@kwdef mutable struct TinkerforgeConnection
  host::IPAddr = ip"127.0.0.1"
  port::Integer = 4223
  connection

  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedGaussMeterParams
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}
end

self.ip_con = IPConnection()
            self.ip_con.connect(self.HOST, self.PORT)


text:
using Documenter, FiniteDifferences

# Copy the README.
repo_root = dirname(dirname(@__FILE__))
cp(
    joinpath(repo_root, "README.md"),
    joinpath(repo_root, "docs", "src", "index.md");
    # `index.md` may already exist if `make.jl` was run previously, so we set `force=true`
    # to overwrite it.
    force=true
)

makedocs(
    modules=[FiniteDifferences],
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages=["Home" => "index.md", "API" => "pages/api.md"],
    authors="Invenia Labs",
    checkdocs=:exports,
    repo="https://github.com/JuliaDiff/FiniteDifferences.jl/blob/{commit}{path}#L{line}",
    sitename="FiniteDifferences.jl",
    strict=true,
)

deploydocs(repo="github.com/JuliaDiff/FiniteDifferences.jl.git")

ext: jl


text:
# Autogenerated wrapper script for HDF5_jll for x86_64-w64-mingw32
export libhdf5, libhdf5_hl

using Zlib_jll
using OpenSSL_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("HDF5")
JLLWrappers.@declare_library_product(libhdf5, "libhdf5-0.dll")
JLLWrappers.@declare_library_product(libhdf5_hl, "libhdf5_hl-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libhdf5,
        "bin\\libhdf5-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhdf5_hl,
        "bin\\libhdf5_hl-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using NeRCA
using Test

@testset "most_frequent()" begin
    a = [1, 1, 2, 3, 1, 5]
    @test 1 == most_frequent(x->x, a)

    a = [[1, 2], [1, 2, 3], [1, 2], [1, 2], [1], [1]]
    @test 3 == most_frequent(sum, a)
end


text:
module MPIDistributedArrays
    using MPI
    export MArray,mzeros,distribute,free
# Write your package code here.

    const masterrank = 0

    function initialize()
        MPI.Init()
    end

    function Comm_size()
        return  MPI.Comm_size(MPI.COMM_WORLD)
    end

    struct MArray{T,N,A} <: AbstractArray{T,N}
        localpart::A#Union{Nothing,A}
        localindices::Array{UnitRange{Int64},1}
        myrank::Int64
        nprocs::Int64
        comm::MPI.Comm
        arrayrank::Array{Int64,1}
        win::MPI.Win
        localindices_set::Array{Array{UnitRange{Int64},1},1}
        size::NTuple{N,Int64}
        localsize::NTuple{N,Int64}

        MArray{T}(dims::Dims,parallel;debug= false) where {T} = MArray{T}(MPI.COMM_WORLD,dims,parallel,debug= debug)
        MArray(::Type{T},dims::Dims,parallel;debug= false) where {T} = MArray{T}(MPI.COMM_WORLD,dims,parallel,debug= debug)
        MArray(::Type{T},comm::MPI.Comm,dims::Dims,parallel;debug= false) where {T} = MArray{T}(comm,dims,parallel,debug= debug)

        function MArray{T}(comm::MPI.Comm,dims::Dims,parallel;debug= false) where {T}
            @assert length(dims) == length(parallel) "dimension mismatch! length(dims) = $(length(dims)) but length(parallel) = $(length(parallel))"
            N = length(dims)
            nprocs = MPI.Comm_size(comm)
            myrank = MPI.Comm_rank(comm)

            #println(typeof(comm))
            #exit()

            if prod(parallel) != nprocs
                if myrank == masterrank
                    println("number of MPI processes is mismatched. nprocs = $nprocs and parallel = $parallel")
                end
                MPI.Barrier(comm)
            end
            arrayrank = get_arrayrank(myrank,parallel)
            localindices = Array{UnitRange{Int64},1}(undef,N)
            localdims = zeros(Int64,N)

            for idim =  1:length(dims)
                @assert dims[idim] % parallel[idim] == 0 "$idim-th dimension: dims[$idim] % parallel[$idim] should be 0! now dims[$idim] = $(dims[idim]) and parallel[$idim] = $(parallel[idim])"
                ndim = dims[idim] ÷ parallel[idim]
                istart = arrayrank[idim]*ndim + 1
                iend = istart + ndim-1
                localindices[idim] = istart:iend
                localdims[idim] = ndim
            end

            localindices_set= Array{Array{UnitRange{Int64},1},1}(undef,nprocs)

            for rank=0:nprocs-1
                localindices_set[rank+1] = copy(localindices)
                MPI.Bcast!(localindices_set[rank+1], rank, comm)

            end

            if debug
                for rank=0:nprocs-1
                    #isinside = get_local_or_not(localindices_set[rank+1],1,1,1)

                    if rank == myrank
                        println("myrank : $myrank localindices = ",localindices)
                        #println("isinside = $isinside")
                        #println("localindices_set: $localindices_set")
                    end
                    MPI.Barrier(comm)
                end
            end

            localpart = Array{T,N}(undef,localdims...)
            win = MPI.Win_create(localpart,comm)
            return new{T,N,Array{T,N}}(localpart,localindices,myrank,nprocs,comm,arrayrank,win,
            localindices_set,dims,size(localpart))

        end

    end

    function free(a::MArray{T,N}) where {T,N}
        MPI.Barrier(a.comm)
        MPI.free(a.win)
    end

    Base.size(a::MArray) = a.size

    function Base.getindex(a::MArray{T,N,A},i...) where {T,N,A}
        isinside,lenind = get_local_or_not(a.localindices,i...)
        if isinside
            index = [i[idim] - a.localindices[idim][1]+1 for idim=1:lenind]
            return a.localpart[index...]
        else
            result = Ref{T}()
            target_rank,local_index = get_targetrank_and_index(a,i...)
            #println("target_rank ",target_rank)
            if target_rank == a.myrank
                #println("target $local_index",a.localpart[local_index])
                result[] = a.localpart[local_index]
            else
                #println("localindex = $local_index")
                MPI.Win_lock(MPI.LOCK_SHARED, target_rank, 0, a.win)
                MPI.Get(result,  target_rank, local_index-1 , a.win)
                MPI.Win_unlock(target_rank, a.win)
            end
            return result[]
        end
    end

    function get_local_or_not(indices,i...)
        #println("indices, $indices i = $i")
        lenind = length(indices)
        isinside = true
        for idim = 1:lenind
            isinside *= i[idim] in indices[idim]
        end
        return isinside,lenind
    end

    function get_targetrank_and_index(A::MArray,i...)
        for rank=0:A.nprocs
            isinside,lenind = get_local_or_not(A.localindices_set[rank+1],i...)
            if isinside
                index = [i[idim] - A.localindices_set[rank+1][idim][1]+1 for idim=1:lenind]
                #println("index = $index")
                prod = 1
                index_1d = 1
                for idim=1:lenind
                    index_1d += prod*(i[idim] - A.localindices_set[rank+1][idim][1]+1-1)
                    prod *= length(A.localindices_set[rank+1][idim])
                end
                return rank,index_1d
            end
        end
    end

    mzeros(comm::MPI.Comm,dims::Dims,parallel;debug=false)  = mzeros(Float64,comm,dims,parallel,debug=debug)
    mzeros(dims::Dims,parallel;debug=false)  = mzeros(Float64,MPI.COMM_WORLD,dims,parallel,debug=debug)
    mzeros(::Type{T},dims::Dims,parallel;debug=false) where {T} = mzeros(T,MPI.COMM_WORLD,dims,parallel,debug=debug)

    function mzeros(::Type{T},comm::MPI.Comm,dims::Dims,parallel;debug=false) where {T}
        A = MArray{T}(comm,dims,parallel,debug= debug)
        A.localpart .= zero(T)
        return A
    end

    distribute(rank::I,A::T,parallel;debug = false) where {T <: AbstractArray, I <: Int} = distribute(rank,MPI.COMM_WORLD,A,parallel,debug = debug)
    distribute(A::T,parallel;debug = false) where {T <: AbstractArray} = distribute(masterrank,MPI.COMM_WORLD,A,parallel,debug = debug)
    distribute(comm::MPI.Comm,A::T,parallel;debug = false) where {T <: AbstractArray} =distribute(masterrank,comm,A,parallel,debug = debug)

    function distribute(rank::I,comm::MPI.Comm,A::T,parallel;debug = false) where {T <: AbstractArray, I <:Int}
        MPI.Bcast!(A, rank, comm)
        A_local = mzeros(eltype(A),size(A),parallel,debug = debug)
        local_length = length(A_local.localpart)
        localindex =  zero(A_local.arrayrank)

        globalindex =  zero(A_local.arrayrank)
        for i=1:local_length
            get_localindex!(localindex,A_local ,i)
            get_globalindex!(globalindex,A_local,localindex,)
            A_local.localpart[localindex...] = A[globalindex...]
            #println("local $localindex",A_local.localpart[localindex...])
            #println("global $globalindex",A[globalindex...])
        end
        return A_local
    end

    function get_arrayrank(myrank,parallel)
        arrayrank = zero(parallel)
        psize = 1
        myranktemp = myrank
        for idim = 1:length(parallel)
            psize = parallel[idim]
            arrayrank[idim] = myranktemp % psize
            myranktemp = (myranktemp-arrayrank[idim]) ÷ psize
        end
        return arrayrank
    end

    function get_myrankfromarray(arrayrank,parallel)
        myrank = 0
        psize = 1
        for idim = 1:length(parallel)
            @assert arrayrank[idim] < parallel[idim]
            myrank += arrayrank[idim]*psize
            psize *= parallel[idim]
        end
        return myrank
    end

    function get_globalindex!(globalindex,A::MArray,localindex)
        for idim = 1:length(localindex)
            #println(A.localindices[idim])
            globalindex[idim] = A.localindices[idim][localindex[idim]]
        end
        return
    end

    function get_localindex!(localindex,A::MArray,i)
        #localindex =  zero(A.arrayrank)
        localsize = size(A.localpart)
        psize = 1
        itemp = i-1
        for idim = 1:length(localindex)
            psize = localsize[idim]
            #println(psize)
            localindex[idim] = itemp % psize
            itemp = (itemp-localindex[idim]) ÷ psize
        end
        localindex .+= 1
        return
        #return localindex .+ 1

    end

    function get_localindex(A::MArray,i)
        localindex =  zero(A.arrayrank)
        get_localindex!(localindex,A::MArray,i)
        return localindex
    end
end


text:
# Copyright (c) 2015 Ishita Takeshi
# License is MIT

# replace multiple whitespaces with single whitespace
strip_line(line) = replace(strip(line), r"\s+", " ")

parsefloat(x) = parse(Float64, x)
parseint(x) = parse(Int64, x)

isnumeric(s::AbstractString) = ismatch(r"[0-9]", s)
iscomment(line) = startswith(line, "#")

"""
Extract a sparse vector of `ndim` dimensions and its label
from the given string line.
If `ndim` is not passed, the vector dimension is automatically determined by
the contents of the given string line.
"""
function line_to_data(line)
    line = strip_line(line)
    splitted = split(line, " ")

    if iscomment(line) || length(line) == 0
        #the line contains only a comment, newline or whitespaces
        throw(NoDataException())
    end

    label = 0  # define label once here
    try
        label = parsefloat(splitted[1])
    catch error
        throw(InvalidFormatError(error.msg))
    end

    if length(splitted) < 2 || startswith(splitted[2], "#")
        # no vector per line or the case such as line = "-1 #comment"
        return (Int64[], Float64[]), label
    end

    indices = Int64[]
    values = Float64[]
    for element in splitted[2:end]
        if startswith(element, "#")
            # regard the remaining characters as a comment
            break
        end

        pair = split(element, ":")

        try
            index = parseint(pair[1])
            value = parsefloat(pair[2])
            push!(indices, index)
            push!(values, value)
        catch error
            throw(InvalidFormatError(error.msg))
        end
    end

    return ((indices, values), label)
end

function load_svmlight_file(filename, ndim = -1)
    I = Int64[]
    J = Int64[]
    V = Float64[]
    y = Array(Float64, 0)

    i = 1
    for line in eachline(open(filename))
        try
            ((indices, values), label) = line_to_data(line)
            append!(I, indices)
            append!(J, ones(length(indices))*i)
            append!(V, values)

            y = push!(y, label)

            i += 1
        catch error
            if isa(error, NoDataException)
                # do nothing
                continue
            end
            throw(error)
        end
    end

    if ndim < 0
        X = sparse(I, J, V)
    else
        X = sparse(I, J, V, ndim, i-1)
    end
    return X, y
end


text:
"""
    MDCProblemJumpStart(cost, p0 dp0, momentum, tspan, jumpsize, reinitialise_dp0)
Do I want to add an option: (reinitialise_dp0 == true) && give option to recalculate hessian
"""
struct MDCProblemJumpStart{A,B,C,D,E,F} <: CurveProblem
    cost::A
    p0::B
    dp0::C
    momentum::D
    tspan::E
    jumpsize::F
    reinitialise_dp0::Bool
    ## reverse initial direction and signflip curve span if the latter is nonpositive
    function MDCProblemJumpStart(a::A, b::B, c::C, d::D, e::E, f::F, g::Bool) where A where B where C where D where E where F
        if max(e...) <= 0.
            e = map(x -> -x |> abs, e) |> reverse
            c = -c
        end
        new{A,B,C,D,E,F}(a, b, c, d, e, f, g)
    end
end
isjumped(c::MDCProblemJumpStart) = JumpStart(c.jumpsize)
whatdynamics(c::MDCProblemJumpStart) = MDCDynamics()

function jumpstart(m::MDCProblem, jumpsize, reinitialise_dp0=false)
    return MDCProblemJumpStart(m.cost, m.p0, m.dp0, m.momentum, m.tspan, jumpsize, reinitialise_dp0)
end

function (j::JumpStart)(m::MDCProblem; reinitialise_dp0=false)
    return MDCProblemJumpStart(m.cost, m.p0, m.dp0, m.momentum, m.tspan, j.jumpsize, reinitialise_dp0)
end

function (c::MDCProblemJumpStart)()
    spans = make_spans(c, c.tspan, ZeroStart())
    cs = map(spans) do span
        mult = sign(span[end])
        return MDCProblemJumpStart(c.cost, c.p0, mult * c.dp0, c.momentum, abs.(span), c.jumpsize, c.reinitialise_dp0)
    end
    spans = map(x -> abs.(x), spans)
    u0s = initial_conditions.(cs)
    fs = dynamics.(cs)
    return ODEProblem.(fs, u0s, spans)
end

function make_spans(c::CurveProblem, span, j::JumpStart)
    spans = make_spans(c, span, ZeroStart())
    spans = map(spans) do span
        (span[1] +  0.1 * sign(span[2]), span[2])
    end
    return spans
end

"""
First need to move p0 in the jumpstart direction. BUT NOT change c.p0.
Then need to make ODEProblem tspan start at the jumpstart time instead of time. ie modify the spans.
Then need to make the costate: first find an initial u, potentially with reinitialise_dp0==true. Then its easy
Then dynamics are the same
"""
function initial_conditions(c::MDCProblemJumpStart)
    θ₀ = c.p0 + get_jump(c)
    λ₀ = initial_costate(c)
    return cat(θ₀, λ₀, dims=1)
end

function get_jump(c::CurveProblem)
    get_jump(c, isjumped(c))
end
get_jump(c, ::ZeroStart) = zero.(param_template(c))
function get_jump(c, j::JumpStart)
    return j.jumpsize * (c.dp0 / norm(c.dp0))
end

function initial_costate(c::MDCProblemJumpStart)
    u = get_initial_velocity(c::MDCProblemJumpStart)
    μ₂ = (-c.momentum + c.cost(c.p0)) / 2.
    λ₀ = -2. * μ₂ * u  # + μ₁*get_jump(c), but μ₁ = 0 by complementary slackness at this point
    return λ₀
end

"""
Algorithm:
Let y = θ - θ₀;
    f(x) = ∇Cᵀx; at θ
    g₁(x) = xᵀx - 1;
    g₂(x) = - xᵀy
Then optimisation problem is:
    minₓ f(x) subject to
    gᵢ(x) ≤ 0.
In other words, find a direction that maximally anticorrelates with the gradient, but has norm ≤ 1 and is pointing away from the curve origin. Norm = 1 would be ideal, but deconvexifies.
KKT conditions give:
    ∇C + 2μ₁x - μ₂y = 0;
    μᵢ ⩾ 0 + complementary slackness
Analytic solution:
    Case 1: ∇Cᵀy ≤ 0. Then μ₂ = 0 from complementary slackness and x ∝ - ∇C
    Case 2: ∇Cᵀ ≥ 0. then
            2μ₁x = μ₂y - ∇C
        ⇒   μ₂ = ∇Cᵀy / yᵀy
        ⇒   ∇Cᵀx ≥ 0 → x = 0
            ∇Cᵀx ≤ 0 → μ₁ s.t. norm(x) = 1
Issues:
- ∇C might be quite noisy close to the minimum. Might want the option of a second order condition involving hessian recalculation for cheaper problems. IE
- What do we do if x = 0? For now, keep the old dp0. Because at least that is in a shallow direction of the old Hessian. In future, could provide option to recalculate Hessian OR reuse Hessian at p0, which would be a good estimate of the new Hessian. If we have the Hessian, the optimisation problem would change:
    f(x) = ∇Cᵀx + 0.5 xᵀ∇²Cx ; at θ
    KKT:
    ∇C + ∇²Cx + 2μ₁x - μ₂y = 0
    (∇²C + 2μ₁𝕀)x = μ₂y - ∇C
And then go through the μᵢ = or ≂̸ 0 cases.
"""
function get_initial_velocity(c::MDCProblemJumpStart)
    (c.reinitialise_dp0 == false) && (return c.dp0)
    ∇C = param_template(c)
    y = get_jump(c)
    new_p0 = c.p0 + y
    c.cost(new_p0, ∇C)
    d = dot(∇C, y)
    if d ≤ 0
        new_dp0 = -∇C/(norm(∇C))
    else
        μ₂ = dot(∇C, y) / sum(abs2,y)
        x = μ₂*y - ∇C; nx = norm(x)
        new_dp0 = x/nx
    end
    if c.cost(new_p0 + 0.001new_dp0) > c.cost(new_p0 + 0.001c.dp0)
        @info("couldn't cheaply find a better initial curve direction dp0 for the jumpstarted problem. will not reinitialise dp0. Provide your own if you want")
        return c.dp0
    else
        @info("cheaply found a better initial curve direction dp0 for the jumpstarted problem. . re-initialising dp0")
        return x / nx
    end
end


text:
using NonstationaryProcesses
using Plots

function wallPlot(S, vars)
        # Square ones
        if vars[2] isa String
                vars, name = vars
        else
                name = getprocess(S)
        end
        mkpath(joinpath(dir, "Square"))
        p = plot(S; vars, colormode=:velocity, linecolor=colorgradient, background_color=:black, axis=nothing, framestyle=:none, size=(900, 900), aspect_ratio=:equal, margin=15Plots.mm, dpi=1000, linewidth=1.0);
        savefig(p, joinpath(dir, "Square", string(name)*".png"))

        # 16:9
        mkpath(joinpath(dir, "WallpaperSquare"))
        p = plot(S; vars, colormode=:velocity, linecolor=colorgradient, background_color=:black, axis=nothing, framestyle=:none, size=(1600, 900), aspect_ratio=:equal, margin=10Plots.mm, dpi=1000, linewidth=1.5);
        savefig(p, joinpath(dir, "WallpaperSquare", string(name)*".png"))

        # 16:9 fill
        mkpath(joinpath(dir, "WallpaperFill"))
        p = plot(S; vars, colormode=:velocity, linecolor=colorgradient, background_color=:black, axis=nothing, framestyle=:none, size=(1600, 900), left_margin=40Plots.mm, right_margin=40Plots.mm, dpi=1000, linewidth=1.5);
        savefig(p, joinpath(dir, "WallpaperFill", string(name)*".png"))
end

colorgradient = :turbo
dir = "../ProcessImages/"
mkpath(dir)

# Dict of system => variables to plot
𝒮 = Dict(
        simplestChaoticFlowVis()                => 2:3,
        thomasCyclicallySymmetricVis()          => 1:2,
        doubleScrollVis()                       => 1:3,
        diffusionlessLorenzVis()                => 1:3,
        piecewiseLinearHyperchaosVis()          => [1, 3],
        simplifiedLorenz4DVis()                 => [3, 2, 4],
        chensSystemVis()                        => [1, 2, 3],
        chensSystemVis(parameters=(46.0, 11.0, 29.0), tmax=500.0)
                                                => ([1, 2, 3], "chensSystem46"),
        cartesianDoublePendulumSim(profiles=(constant, constant, constant, constant),
                                   parameters=(0.5, 1.0, 1.0, 1.0),
                                   X0=[π/4, π, 0.0, 0.0],
                                   tmax=2000.0,
                                   savedt= 0.01)
                                                => [3, 4],
        lorenzVis()                             => 1:3,
)

for (S, vars) ∈ 𝒮
        wallPlot(S, vars)
end


text:
## Distribution Test
## =================
holesizes = logspace(-1,-10,3)

SplitBenchmarks = [Simplex(i,[:sample_distribution,
                                        :accumulative_KL,
                                        :total_time,],1,0.1)
                   for i = 2:6]

mh_captures = [:samples, :refine, :depth]
all_splits = [rand_partial_split]

SMTalgorithms = [SigmaSMT(mh_captures, solver, sampler, nprocs, split)
  for nprocs = [1],
      solver = [DRealSolverBinary],
      split = all_splits,
      sampler = [Sigma.pre_tlmh]][:]

Rejalgorithms = [SigmaRej(mh_captures)]

Bfsalgorithms = [SigmaBFS(mh_captures)]

function simplex_dim()
  @show "Got to start of dim"
  record(Bfsalgorithms, SplitBenchmarks;
         runname = "simplex_dim",prefix=benchdir,savedb=false,exceptions=true)
end

function count(ndims)
  results = DynamicAnalysis.all_records() |> r->DynamicAnalysis.where(r,:runname,j->j=="kl") |> r->DynamicAnalysis.where(r,:problem,j->j.ndims == ndims)
  @show size(results,1)

end

# function accumulative_kl(i, ndims)
#   records = DynamicAnalysis.all_records() |> r->DynamicAnalysis.where(r,:runname,j->j=="kl") |> r->DynamicAnalysis.where(r,:problem,j->j.ndims == ndims)
#   results = records[:results]
#   @show length(results)
#   allsamples = Lens.get(results[i];lensname=:samples)
#   kls = Float64[]
#   for i = 1:length(allsamples)
#     samples = allsamples[1:i]
#     v = vertex_distribution(samples,ndims,0.01)
#     truth = groundtruth(ndims)
#     push!(kls,KLsafe(truth,v))
#   end
#   kls
# end

# function accumulative_mean_kls(dims)
#   trials = [accumulative_kl(i, dims) for i = 1:count(dims)]
#   mean_kls = Array(Float64, length(trials[1]))
#   stds_kls = Float64[]
#   for i = 1:length(trials[1])
#     total = 0.0
#     for j = 1:length(trials)
#       total += trials[j][i]
#     end
#     mean_kls[i] = total/length(trials)
#   end
#   return mean_kls
# end

# function plot_kls(dims::Vector{Int})
#   dfs = DataFrame[]
#   for dim in dims
#     akl = accumulative_mean_kls(dim)
#     df = DataFrame(x=1:length(akl), y = akl, label = "Dim:$dim")
#     push!(dfs,df)
#   end
#   Gadfly.plot(vcat(dfs...), x="x", y="y", color="label",Gadfly.Geom.line,
#      Guide.xlabel("Number of Samples"),
#      Guide.ylabel("KL Divergence"))
# end

# splitkey = [rand_partial_split => "rand", weighted_mid_split => "mid", weighted_partial_split => "partial"]

# # # Analysis
# # ## KL
# using Lens
# using SQLite
# using Gadfly
# using DataFrames
# q = all_records() |> r->Lens.where(r,"runname",j->j=="kl")
# g = group(q)

# plot_db(y = tree, x = :iterate, groupby=)

# function plot_kl2(groups)
#   dfs = DataFrame[]
#   i = 1
#   for (k,v) in groups
#     results = Result[x[7] for x in v]
#     collated = Lens.collate(results, :samples, :x1)
#     @show length(collated)
#     samples = collated
#     ndims = k[2].ndims
#     holesize = k[2].holesize
#     gt = groundtruth(ndims)
#     akl = accumulative_KL(samples, ndims, gt,holesize)
#     df = DataFrame(x=1:length(akl), y = akl, label = "SMT$i")
#     push!(dfs,df)
#     i += 1
#   end
# #   return  dfs
#   plot(vcat(dfs...), x="x", y="y", color="label",Geom.line,
#        Guide.xlabel("Number of Samples"),
#        Guide.ylabel("KL Divergence"))
# end

# plot_kl2(g)

# draw(PDF("/home/zenna/Dropbox/writing/conferences/icml2015/figures/smt1.pdf",4*1inch,3*1inch), smtplot)


text:
module GeometryReferenceGenericSolver

using LinearAlgebra: norm
using LinearAlgebra: kron, I
using SparseArrays: sparse

using Parametron
using Parametron:setdirty!
using OSQP
using OSQP.MathOptInterfaceOSQP: OSQPSettings

import MathOptInterface
const MOI = MathOptInterface

function _constrain_convex_combination_coefficients(model, coeffs)
    M = length(coeffs)
    for i = 1:M
        @constraint(model, coeffs[i] >= 0)
    end
    s = sum(coeffs)
    @constraint(model, s == 1)
end

function _convex_combination(coeffs, points, m)
    s = @expression coeffs[1] * points[:, 1]
    for i = 2:m
        s = @expression s + coeffs[i] * points[:, i]
    end
    return s
end

function squared_distance(p1, p2, n)
    s = @expression (p1[1] - p2[1])^2
    for i = 2:n
        s = @expression s + (p1[i] - p2[i])^2
    end
    return s
end

struct SolverError{T} <: Exception
    status::T
end

function find_minimum_distance_parametron(convex1::NTuple{M1, PointT}, convex2::NTuple{M2, PointT}) where {M1, M2, PointT}
    model = Model(OSQP.Optimizer(
        max_iter=4000,
        warm_start=false,   # for deterministic results
        polish=true,
        eps_prim_inf=1e-8,
        eps_dual_inf=1e-8,
        eps_abs=1e-8,
        eps_rel=1e-6,
        scaling=0)
    )

    # from https://github.com/tkoolen/Parametron.jl/blob/d824c9f9ff4ba7bc38a9407129504bb8ecca31ee/test/model.jl#L22
    #MOI.set(model.optimizer, OSQPSettings.AdaptiveRhoInterval(), 25) # required for deterministic behavior
    #MOI.set(model.optimizer, OSQPSettings.Verbose(), false)
    #MOI.set(model.optimizer, OSQPSettings.Scaling(), false)
    N = length(PointT)      # TODO works for GeometryType.Point, general trait?
    # points
    c1 = Parameter(model, val=zeros(N, M1))
    c2 = Parameter(model, val=zeros(N, M2))

    # TODO use broadcast machinery to set ]up problem
    # coefficients of convex combination
    a1 = [Variable(model) for _ in 1:M1]
    a2 = [Variable(model) for _ in 1:M2]

    _constrain_convex_combination_coefficients(model, a1)
    _constrain_convex_combination_coefficients(model, a2)

    p1 = _convex_combination(a1, c1, M1)
    p2 = _convex_combination(a2, c2, M2)
    d_sq = squared_distance(p1, p2, N)

    @objective(model, Minimize, d_sq)

    # set parameters
    for i = 1:M1
        c1.val[][:,i] = convex1[i]
    end
    for i = 1:M2
        c2.val[][:,i] = convex2[i]
    end

    setdirty!(c1)
    setdirty!(c2)
    #=data_before_solve = deepcopy(
     unsafe_load_data(
        unsafe_load(
            unsafe_load(model.optimizer.inner.workspace).data))
    )
    =#
    data_before_solve = "syke"
    solve!(model)
    status =  MOI.get(model.optimizer, MOI.TerminationStatus())

    if status != MOI.OPTIMAL
        #throw(SolverError(status))
    end

    # TODO remove after done debugging
    g = value(model, a1[1])

    # extract decision varibles and re-compute expressions
    # TODO avoid expressing convex combinations and norm twice
    v_a1 = value.(model, a1)
    v_a2 = value.(model, a2)

    v_p1 = sum(v_a1 .* convex1)
    v_p2 = sum(v_a2 .* convex2)
    v_d_sq = sum((v_p1 - v_p2).^2)
    return (d_sq = v_d_sq, p1=v_p1, p2=v_p2, model=model,
        model_parts = (
            a1=a1, a2=a2, c1=c1, c2=c2,
        ),
        data_before_solve=data_before_solve
    )
end

"""
produce U from A such that `x' * U * x == x' * A * x` ∀ x
overwrite A in-place
"""
function upper_triangularize_quadratic!(A)
  s = size(A)
  n = first(s)
  @assert s == (n, n)

  # iterate below diagonal
  for i=2:n, j=1:(i-1)
    A[j,i] += A[i,j]
    A[i,j] = 0
  end
  return nothing
end

"""
Make QP matrices for the problem of finding the minimum distance between two convex hulls.
"""
function setup_qp_data(c1, c2)
  # order all a1 and then all a2
  m1 = length(c1)
  m2 = length(c2)
  n = length(c1[1]) # assume all same dimension

  M = [m1, m2]
  M_offset = [0, m1]

  _C = [c1, c2]

  C = zeros(2*n, m1 + m2)
  for point_set_i = 1:2
    for _j = 1:M[point_set_i]
      j = _j + M_offset[point_set_i]
      n_offset = (point_set_i-1) * n
      C[(1:n) .+ n_offset, j] .= _C[point_set_i][_j]
    end
  end

  p1_minus_p2_op = kron([1, -1]', I(n))

  # maps convex combination coefficients to the difference
  p1_minus_p2__Ta = p1_minus_p2_op * C

  QP_P = 2 * (p1_minus_p2__Ta' * p1_minus_p2__Ta)

  n_constraints = (m1 + m2 + length(M))
  n_variables = m1 + m2

  A = zeros(n_constraints, n_variables)
  l = zeros(n_constraints)
  u = zeros(n_constraints)

  A[1:n_variables, 1:n_variables] = I(n_variables)
  l[1:n_variables] .= 0
  u[1:n_variables] .= Inf

  A[n_variables+1, (1:m1) ] .= ones(m1)
  A[n_variables+2, (1:m2) .+ m1 ] .= ones(m2)
  l[(n_variables+1):end] .= 1
  u[(n_variables+1):end] .= 1
  #upper_triangularize_quadratic!(QP_P)
  return (
    A=A, l=l, u=u, P=QP_P, q=zeros(n_variables),
    parts=(
        m1=m1,
        m2=m2,
        p1_minus_p2__Ta=p1_minus_p2__Ta,
        C = C,
    )
  )
end

quad(x) = x' * x

function setup_model(data)
  model = OSQP.Optimizer().inner

  #model = OSQP.Model()
  #MOI.set(model.optimizer, OSQPSettings.Scaling(), false)
  #OSQP.update_settings!(model; scaling=0)

  OSQP.setup!(model;
          P=sparse(data.P),
          q=data.q,
          A=sparse(data.A),
          l=data.l,
          u=data.u,

          verbose=false,
          max_iter=4000,
          warm_start=false,   # for deterministic results
          polish=true,
          eps_prim_inf=1e-8,
          eps_dual_inf=1e-8,
          eps_abs=1e-8,
          eps_rel=1e-6,
          scaling=10,
          )
  model
end

function find_minimum_distance(c1, c2)
  qp_problem_data = setup_qp_data(c1, c2)

  model = setup_model(qp_problem_data)
  result = OSQP.solve!(model)
  coeffs = result.x
  a1 = coeffs[1:qp_problem_data.parts.m1 .+ 0]
  a2 = coeffs[(1:qp_problem_data.parts.m2) .+ qp_problem_data.parts.m1]
  p1 = sum(c1 .* a1)
  p2 = sum(c2 .* a2)

  return (
    d_sq=result.info.obj_val, p1=p1, p2=p2,
    model=model, result=result)
end

end # module

max_stars_repo_name: goretkin/GeometryReferenceGenericSolver

max_issues_repo_name: goretkin/GeometryReferenceGenericSolver

max_forks_repo_name: goretkin/GeometryReferenceGenericSolver


text:
using XGBoost
using Random: bitrand
using SparseArrays: sparse
using Test

include("utils.jl")

@testset "XGBoost" begin

@testset "Sparse matrices" begin
    X = sparse(randn(100,10) .* bitrand(100,10))
    y = randn(100)
    DMatrix(X, label = y)

    X = sparse(convert(Matrix{Float32}, randn(10,100) .* bitrand(10,100)))
    y = randn(100)
    DMatrix(X, true)

    X = sparse(randn(100,10) .* bitrand(100,10))
    y = randn(100)
    DMatrix(X)
end

@testset "DMatrix loading" begin
    dtrain = DMatrix("../data/agaricus.txt.train")
    train_X, train_Y = readlibsvm("../data/agaricus.txt.train", (6513, 126))
    @test dtrain != nothing

    labels = get_info(dtrain, "label")

    @test train_Y == labels
end

@testset "Agaricus training" begin
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = xgboost(dtrain, 2, watchlist=watchlist, eta = 1, max_depth = 2,
                  objective = "binary:logistic", silent = 1)
    @test bst != nothing

    preds = XGBoost.predict(bst, dtest)

    labels = get_info(dtest, "label")
    @test size(preds) == size(labels)

    err = count(!iszero, (preds .> 0.5) .!= labels) / length(preds)
    @test err < 0.1
end

@testset "Cross validation" begin
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = nfold_cv(dtrain, 5, 3, eta = 1, max_depth = 2, objective = "binary:logistic", silent = 1,
                   seed = 12345)
    # important_features = importance(bst)
    #
    # @test startswith(important_features[1].fname, "f28")
    # @test important_features[1].fname == "f28" # pending
end

@testset "Feature importance" begin
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = xgboost(dtrain, 5, watchlist = watchlist, eta = 1, max_depth = 2,
                  objective = "binary:logistic", silent = 1, seed = 12345)
    important_features = importance(bst)

    @test startswith(important_features[1].fname, "f28")
    # @test important_features[1].fname --> "f28" # pending
end

@testset "Example" begin
    include("example.jl")
end

end


text:
export
    sqrtmod,
    multiplicativeorder,
    linearmod,
    primitiveroot

Nemo.sqrtmod(x::Integer, m::T) where T<:Integer = convert(T, sqrtmod(fmpz(x), fmpz(m)))

#=
http://rosettacode.org/wiki/Multiplicative_order

One possible algorithm that is efficient also for large numbers is the
following: By the Chinese Remainder Theorem, it's enough to calculate the
multiplicative order for each prime exponent p^k of m, and combine the
results with the least common multiple operation. Now the order of a wrt.
to p^k must divide Φ(p^k). Call this number t, and determine it's factors
q^e. Since each multiple of the order will also yield 1 when used as
exponent for a, it's enough to find the least d such that (q^d)*(t/(q^e))
yields 1 when used as exponent.
=#
multiplicativeorder(a::Integer, m::Integer) = multiplicativeorder(promote(a, m)...)
multiplicativeorder(a::T, m::T) where {T<:Integer} = begin
    iscoprime(a, m) || error("$a and $m are not coprime")

    res = one(m)
    for (p, k) in factor(m)
        m = p^k
        t = (m ÷ p) * (p-1)
        for f in factors(t)
            if powermod(a, f, m) == 1
                res = lcm(res, f)
                break
            end
        end
    end
    res
end

##  Solve linear modulo equation a * x = b mod m
# https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl
linearmod(a::Integer, b::Integer, m::Integer) = linearmod(promote(a, b, m)...)
linearmod(a::T, b::T, m::T) where {T<:Integer} = begin
    m > 1 || error("Argument 'm' must be an integer greater 1")

    g, u, v = gcdx(a, m)
    mod(b, g) == 0 || return T[]

    x0 = mod(u * (b ÷ g), m)
    [mod(x0 + i * (m ÷ g), m) for i in 0:(g-1)]
end

##  Find a primitive root modulo m
# https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl
primitiveroot(m::T) where {T<:Integer} = begin
    isprime(m) || error("Argument 'm' must be a prime number")
    m == 2 && return one(T)

    P = primefactors(m-1)
    for r in 2:(m-1)
        not_found = true
        for p in P
            if powermod(r, (m-1) ÷ p, m) == 1
                not_found = false
                break
            end
        end
        not_found && return r
    end
    zero(T)
end


text:
using MathOptInterface
const MOI = MathOptInterface

"""MOI utility to create a square matrix from an upper-triangle vector"""
function square_a_triangle(f::MOI.AbstractVectorFunction, dimension::Integer)
        x = MOIU.eachscalar(f)
        S = MOIU.scalar_type(typeof(f))
        matrix = Matrix{S}(undef, dimension, dimension)
        k = 0
        for j in 1:dimension, i in 1:j
                k += 1
                matrix[j, i] = matrix[i, j] = x[k]
        end
        matrix
end


text:
# this script will precompile some commonly used functions
using Econometrics, StatsPlots, CSV, DataFrames
ols()
mleresults()
gmmresults()
samin()
fminunc()
fmincon()
mcmc()
npreg()
dstats(rand(100,3))
plot(1:10)
histogram(rand(100))
scatter(1:10)
d = DataFrame(rand(2,2), :auto)
CSV.write("/tmp/junk.csv", d)
d = CSV.read("/tmp/junk.csv", DataFrame)
clc()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ bab76e0b-52ce-413c-a0e5-3d2d9722d78b
begin
	using CairoMakie, Makie.GeometryBasics
end

# ╔═╡ 378bca1c-d4b1-11ec-144c-ef6ccaf37a63
md"# cairomakie circle"

# ╔═╡ 2cbbc96a-2dc0-4bc5-b22b-20e5be985745
md"used for test"

# ╔═╡ 64d89c8c-8627-4559-9eda-5bb10d24d291
let
	x = range(0, 10, length=100)
	y1 = sin.(x)
	y2 = cos.(x)

	lines(x, y1, color = (:red,0.5))
	lines!(x, y2, color = (:blue,0.7))
	current_figure()
end

# ╔═╡ Cell order:
# ╠═378bca1c-d4b1-11ec-144c-ef6ccaf37a63
# ╠═bab76e0b-52ce-413c-a0e5-3d2d9722d78b
# ╠═2cbbc96a-2dc0-4bc5-b22b-20e5be985745
# ╠═64d89c8c-8627-4559-9eda-5bb10d24d291

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_name: phpsmarter/plutonote


text:
module GenomicRegion
    include ("region.jl")

    include ("bed.jl")
    include ("bedgraph.jl")

    include ("tabix.jl")
end


text:
"""
	coords(X; columns=nothing)

Reshape and filter dataframe or matrix and return the matrix in right format
for unfolding.

## Parameters:

* `X`       - the input Dataframe or matrix with the coordinates.
* `columns` - if the input is a Dataframe with more than 3 columns, the column
  names of the three coordinates must be informed here (e.g. ["X","Y","Z"]).
"""
function coords(x; columns=nothing)
	!isnothing(columns) && (x = x[:,Symbol.(columns)])
	shape = size(x)
	@assert (3 in shape) "3-D coordinate matrix should be informed"
	!(x isa Matrix) && (x = Matrix{Float64}(x))
	(shape[2]==3 && shape[1]!=3) && (x = permutedims(x))
	x
end

"""
	to_vtk(outname, coordinates, props)

Export the `coordinates` as VTK points named `outname`.vtu. Optionally, extra
properties are passed via `props` as `NamedTuple` or `Dict`, where the first item
is the property name and the second is an array with the property values.
"""
function to_vtk(outname::String, coords::AbstractMatrix, props=nothing)
	verts = [MeshCell(VTKCellTypes.VTK_VERTEX, [i]) for i in 1:size(coords,2)]
	outfiles = vtk_grid(outname, coords, (verts)) do vtk
		if !isnothing(props)
			for (name,prop) in zip(keys(props), values(props))
				vtk[string(name)] = prop
			end
		end
	end
end

# get neighbors
function get_neighbors(origin::AbstractMatrix, target::AbstractMatrix,
	                   nhood::Symbol, neigh_val, calcdists=false)
	@assert nhood in [:knn,:radius] "invalid neighborhood type"
	tree = KDTree(origin)
	if nhood==:knn
		nneigh = minimum([neigh_val, size(origin,2)])
		idxs, dists = knn(tree, target, nneigh, true)
		idxs, dists
	elseif nhood==:radius
		idxs = inrange(tree, target, neigh_val)
		if calcdists
			dists = Vector{Vector{Float64}}(undef, 0)
			for i in 1:length(idxs)
				p1 = view(target,:,i)
				d  = [euclidean(p1, view(origin,:,j)) for j in idxs[i]]
				push!(dists,d)
			end
		else
			dists = nothing
		end
		idxs, dists
	end
end

get_neighbors(ref_coords::AbstractMatrix, nhood, neigh_val, calcdists=false) =
	get_neighbors(ref_coords, ref_coords, nhood, neigh_val, calcdists)


text:
function _generate_apt_install_command(config::Config)
    apt = config.apt
    pkgs = config.pkgs
    if isempty(config.apt)
        return ""
    else
        apt_list = join(apt, " ")
        return "apt-get install -yq --no-install-recommends $(apt_list)"
    end
end

function _generate_tests_must_pass_commands(config::Config)::Vector{String}
    tests_must_pass = config.tests_must_pass
    all_commands = String[]
    for pkg_name in tests_must_pass
        push!(all_commands,
              "cd /tmp; JULIA_DEBUG=all /usr/bin/julia -e 'println(:current_stage_is_the_tests_must_pass_stage); println(:now_testing_$(pkg_name)); import Pkg; Pkg.test(string(:$(pkg_name)))'")
    end
    return all_commands
end

function Template(config::Config)
    apt_install = _generate_apt_install_command(config)
    julia_url, asc_url = _get_julia_url(config)
    parent_image = config.parent_image
    tests_must_pass_commands = _generate_tests_must_pass_commands(config)

    file_list_vector = File[
        File(_generate_backupsofsimplecontainergenerator_content_1(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "backups_of_simplecontainergenerator_1.jl"),
        File(_generate_backupsofsimplecontainergenerator_content_2(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "backups_of_simplecontainergenerator_2.jl"),
        File(_generate_install_packages_content(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "install_packages.jl"),
        File(_generate_use_sysimage_julia_script_content(config),
             "/tmp/staging/copyfiles/usr/bin",
             "julia"),
        File(_generate_do_not_use_sysimage_julia_script_content(config),
             "/tmp/staging/copyfiles/usr/bin",
             "no_sysimage_julia"),
        File(_generate_packagecompiler_install_content(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "packagecompiler_install.jl"),
        File(_generate_packagecompiler_run_content(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "packagecompiler_run.jl"),
        File(_generate_precompile_execution_content(config),
             "/tmp/staging/copyfiles/opt/simplecontainergenerator_containers",
             "precompile_execution.jl"),
        File(_generate_global_startup_file_content(config),
             "/tmp/staging/copyfiles/opt/etc/julia",
             "startup.jl"),
    ]

    step_list_vector = Vector{AbstractStep}(undef, 0)
    #
    push!(step_list_vector, DockerOnlyLine("FROM $(parent_image)"))
    #
    push!(step_list_vector, EnvStep("DEBIAN_FRONTEND", "noninteractive"))
    #
    push!(step_list_vector, DockerOnlyLine("RUN apt-get update"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get -yq dist-upgrade"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get update"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get -yq dist-upgrade"))
    push!(step_list_vector, DockerOnlyLine("RUN $(apt_install)"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get update"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get -yq dist-upgrade"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get update"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get -yq dist-upgrade"))
    push!(step_list_vector, DockerOnlyLine("RUN apt-get clean"))
    #
    push!(step_list_vector, PlaceholderYumStep())
    #
    push!(step_list_vector, RunStep("rm -rf /var/lib/apt/lists/*"))
    push!(step_list_vector, DockerOnlyLine("RUN echo \"en_US.UTF-8 UTF-8\" > /etc/locale.gen && locale-gen"))
    #
    push!(step_list_vector, EnvStep("SHELL", "/bin/bash"))
    push!(step_list_vector, EnvStep("LC_ALL", "en_US.UTF-8"))
    push!(step_list_vector, EnvStep("LANG", "en_US.UTF-8"))
    push!(step_list_vector, EnvStep("LANGUAGE", "en_US.UTF-8"))
    push!(step_list_vector, EnvStep("PROMPT_COMMAND", "\"echo -n [\$()]\""))
    #
    push!(step_list_vector, RunStep("mkdir -p /tmp/simplecontainergenerator_containers-julia-gpg-key"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-julia-gpg-key && curl https://julialang.org/assets/juliareleases.asc --output juliareleases.asc"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-julia-gpg-key && gpg --import juliareleases.asc"))
    push!(step_list_vector, RunStep("rm -rf /tmp/simplecontainergenerator_containers-julia-gpg-key"))
    #
    push!(step_list_vector, RunStep("mkdir -p /tmp/simplecontainergenerator_containers-download-julia"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-download-julia && wget -O julia.tar.gz $(julia_url)"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-download-julia && wget -O julia.tar.gz.asc $(asc_url)"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-download-julia && gpg --verify julia.tar.gz.asc"))
    push!(step_list_vector, RunStep("mkdir -p /opt"))
    push!(step_list_vector, RunStep("cd /tmp/simplecontainergenerator_containers-download-julia && tar xzf julia.tar.gz -C /opt --strip-components=1"))
    push!(step_list_vector, RunStep("rm -rf /tmp/simplecontainergenerator_containers-download-julia"))
    #
    push!(step_list_vector, RunStep("rm -rf /tmp/simplecontainergenerator-containers-build-depot"))
    push!(step_list_vector, RunStep("mkdir -p /tmp/simplecontainergenerator-containers-build-depot"))
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEPOT_PATH=/tmp/simplecontainergenerator-containers-build-depot /opt/bin/julia -e 'import InteractiveUtils; InteractiveUtils.versioninfo(; verbose=true)'"))
    push!(step_list_vector, RunStep("rm -rf /tmp/simplecontainergenerator-containers-build-depot"))
    #
    push!(step_list_vector, RunStep("rm -rf /opt/etc/julia/startup.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/share/julia/compiled"))
    push!(step_list_vector, RunStep("rm -rf /usr/bin/julia"))
    push!(step_list_vector, RunStep("mkdir -p /opt/etc/julia"))
    #
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/backup_of_dockerfile"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/depot_backup_simplecontainergenerator"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/depot_backup_packagecompiler"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/packagecompiler_depot"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/julia_depot"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/julia_project"))
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers/sysimage"))
    #
    push!(step_list_vector, RunStep("mkdir -p /usr/bin"))
    #
    push!(step_list_vector, PlaceholderDockerCopyFiles())
    #
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_1.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_2.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/install_packages.jl"))
    push!(step_list_vector, RunStep("rm -rf /usr/bin/julia"))
    push!(step_list_vector, RunStep("rm -rf /usr/bin/no_sysimage_julia"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/packagecompiler_install.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/packagecompiler_run.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/precompile_execution.jl"))
    push!(step_list_vector, RunStep("rm -rf /opt/etc/julia/startup.jl"))
    #
    push!(step_list_vector, RunStep("mkdir -p /opt/simplecontainergenerator_containers"))
    push!(step_list_vector, RunStep("mkdir -p /usr/bin"))
    push!(step_list_vector, RunStep("mkdir -p /opt/etc/julia"))
    #
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_1.jl /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_1.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_2.jl /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_2.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/install_packages.jl /opt/simplecontainergenerator_containers/install_packages.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/usr/bin/julia /usr/bin/julia"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/usr/bin/no_sysimage_julia /usr/bin/no_sysimage_julia"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/packagecompiler_install.jl /opt/simplecontainergenerator_containers/packagecompiler_install.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/packagecompiler_run.jl /opt/simplecontainergenerator_containers/packagecompiler_run.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/simplecontainergenerator_containers/precompile_execution.jl /opt/simplecontainergenerator_containers/precompile_execution.jl"))
    push!(step_list_vector, RunStep("cp /tmp/staging/copyfiles/opt/etc/julia/startup.jl /opt/etc/julia/startup.jl"))
    #
    push!(step_list_vector, RunStep("rm -rf /tmp/staging"))
    #
    push!(step_list_vector, RunStep("chmod 444 /opt/etc/julia/startup.jl"))
    push!(step_list_vector, RunStep("chmod 555 /usr/bin/julia"))
    push!(step_list_vector, RunStep("chmod 555 /usr/bin/no_sysimage_julia"))
    #
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all SIMPLECONTAINERGENERATOR_CONTAINER_NO_TEMP_DEPOT=\"true\" /usr/bin/no_sysimage_julia /opt/simplecontainergenerator_containers/install_packages.jl"))
    #
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all /opt/bin/julia /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_1.jl"))
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all /opt/bin/julia /opt/simplecontainergenerator_containers/backups_of_simplecontainergenerator_2.jl"))
    #
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all /opt/bin/julia /opt/simplecontainergenerator_containers/packagecompiler_install.jl"))
    push!(step_list_vector, RunStep("cd /tmp; SIMPLECONTAINERGENERATOR_CONTAINER_NO_TEMP_DEPOT=\"true\" /usr/bin/no_sysimage_julia /opt/simplecontainergenerator_containers/packagecompiler_run.jl"))
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all /usr/bin/no_sysimage_julia -e 'import InteractiveUtils; InteractiveUtils.versioninfo(; verbose=true)'"))
    push!(step_list_vector, RunStep("cd /tmp; JULIA_DEBUG=all /usr/bin/julia -e 'import InteractiveUtils; InteractiveUtils.versioninfo(; verbose=true)'"))
    #
    for cmd in tests_must_pass_commands
        push!(step_list_vector, RunStep(cmd))
    end
    #
    push!(step_list_vector, RunStep("rm -rf /opt/bin/julia/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/etc/julia/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/share/julia/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/depot_backup_simplecontainergenerator/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/depot_backup_packagecompiler/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/packagecompiler_depot/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/julia_depot/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/julia_project/compiled"))
    push!(step_list_vector, RunStep("rm -rf /opt/simplecontainergenerator_containers/sysimage/compiled"))
    #
    push!(step_list_vector, RunStep("find /opt -type d -print0 | xargs -0 chmod a+rx"))
    push!(step_list_vector, RunStep("find /opt -type f -print0 | xargs -0 chmod a+r"))
    push!(step_list_vector, RunStep("chmod a+rx /opt/bin/julia"))
    push!(step_list_vector, RunStep("chmod a+rx /usr/bin/julia"))
    push!(step_list_vector, RunStep("chmod a+rx /usr/bin/no_sysimage_julia"))
    #
    push!(step_list_vector, DockerOnlyLine("RUN mkdir -p /opt/simplecontainergenerator_containers/backup_of_dockerfile"))
    push!(step_list_vector, DockerOnlyLine("COPY Dockerfile /opt/simplecontainergenerator_containers/backup_of_dockerfile/Dockerfile"))
    #
    push!(step_list_vector, DockerOnlyLine("ENTRYPOINT [\"/bin/bash\", \"-c\"]"))

    result = Template(FileList(file_list_vector),
                      StepList(step_list_vector))
    return result
end


text:
"""
CTPG_train(dynamics_plant::Function, dynamics_controller::Function, cost_running::Function, cost_terminal::Function, cost_regularisor::Function, policy_NN, scenario;
solve_alg = Tsit5(), sense_alg = InterpolatingAdjoint(autojacvec = ZygoteVJP()), ensemble_alg = EnsembleThreads(), opt_1 = ADAM(0.01), opt_2 = LBFGS(), maxiters_1 = 100, maxiters_2 = 100, progress_plot = true, solve_kwargs...)

`CTPG_train()` provides a high-level interface for optimisation of the neural networks inside an ODE-represented dynamics based on Continuous-Time Policy Gradient (CTPG) methods that belong to the adjoint sensitivity analysis techniques. The code implemented and the default values for keyword arguments are specified considering training of a neural controller as the main application. In the context herein, a neural controller refers to a dynamic controller that incorporates neural-network-represented components at some points in its mathematical description.

The code utilises the functionalities provided by the [DiffEqFlux.jl](https://github.com/SciML/DiffEqFlux.jl) and [DiffEqSensitivity.jl](https://github.com/SciML/DiffEqSensitivity.jl) packages, and the Automatic Differentiation (AD) capabilities provided by the [Zygote.jl](https://github.com/FluxML/Zygote.jl) package that is integrated in DiffEqFlux.jl. `CTPG_train()` presumes the consistency of the functions provided as its input arguments with the AD tool, hence, the dynamics and cost functions should maintain their transparence against AD tools.

The optimisation (training) problem minimises the cost function defined over deterministic samples of the initial plant state `x₀` and the reference `r` by performing ensemble simulation based on parallelised computation.

The signals are defined as described below:

- `t`: time
- `x`: plant state
- `y`: plant output (= sensor output)
- `x_c`: controller state
- `u`: plant input (= controller output)
- `r`: exogenous reference
- `x_aug`: augmented forward dynamics state (= `[x; x_c; ∫cost_running]`)
- `p_NN`: neural network parameter

The arguments should be provided as explained below:

- `dynamics_plant`: Describes the dynamics of the plant to be controlled. Input arguments `x` and `u` should be of Vector type.
- `dynamics_controller`: Describes the dynamics of the controller that includes neural networks components. Input arguments `x_c`, `y`, `r`, and `p_NN` should be of Vector type.
- `dynamics_sensor`: Describes the dynamics of the sensor that measures output variables fed to the controller. Input arguments `x` should be of Vector type:
- `cost_running`: Describes the running cost defined as the integrand of the Lagrange-form continuous functional. Input arguments `x`, `y`, `u`, and `r` should be of Vector type.
- `cost_terminal`: Describes the terminal cost defined as the Mayer-form problem cost function. Defines a Bolza-form problem along with `cost_running`. Input arguments `x_f` and `r` should be of Vector type.
- `cost_regularisor`: Describes the regularisation term appended to the cost (loss) function. Input argument `p_NN` should be of Vector type.
- `policy_NN`: The neural networks entering into the controller dynamics. DiffEqFlux-based FastChain is recommended for its construction.
- `scenario`: Contains the parameters related with the ensemble-based training scenarios.
    - `ensemble`: A vector of the initial plant state `x₀` and the reference `r` constituting the trajectory realisations.
    - `t_span`: Time span for forward-pass integration
    - `t_save`: Array of time points to be saved while solving ODE. Typically defined as `t_save = t_span[1]:Δt_save:t_span[2]`
    - `dim_x`: `length(x)`
    - `dim_x_c`: `length(x_c)`

The keyword arguments should be provided as explained below:

- `solve_alg`: The algorithm used for solving ODEs. Default value is `Tsit5()`
- `sense_alg`: The algorithm used for adjoint sensitivity analysis. Default value is `InterpolatingAdjoint(autojacvec = ZygoteVJP())`, because the control problems usually render the `BacksolveAdjoint()` unstable. The vjp choice `autojacvec = ReverseDiffVJP(true)` is usually faster than `ZygoteVJP()`, when the ODE function does not have any branching inside. Please refer to the [DiffEqFlux documentation](https://diffeqflux.sciml.ai/dev/ControllingAdjoints/) for further details.
- `ensemble_alg`: The algorithm used for handling ensemble of ODEs. Default value is `EnsembleThreads()` for multi-threaded computation in CPU.
- `opt_1`: The algorithm used for the first phase of optimisation which rapidly delivers the parameter to a favourable region around a local minimum. Default value is `ADAM(0.01)`.
- `opt_2`: The algorithm used for the second phase of opitmisaiton. Defalut value is `LBFGS()` which refines the result of the first phase to find a more precise minimum. Please refer to the [DiffEqFlux documentation](https://diffeqflux.sciml.ai/dev/sciml_train/) for further details about two-phase composition of optimisers.
- `maxiters_1`: The maximum number of iterations allowed for the first phase of optimisation with `opt_1`. Defalut value is `100`.
- `maxiters_2`: The maximum number of iterations allowed for the second phase of optimisation with `opt_2`. Defalut value is `100`.
- `progress_plot`: The indicator to plot the state history for a nominal condition among the ensemble during the learning process. Default value is `true`.
- `i_nominal`: The index to select the case to plot using `progress_plot` during optimisation process from the `ensemble` defined in `scenario`. Defalut value is `nothing`.
- `p_NN_0`: Initial value of the NN parameters supplied by the user to bypass random initialisation of `p_NN` or to continue optimisation from the previous result. Defalut value is `nothing`.
- `solve_kwargs...`: Additional keyword arguments that are passed onto the ODE solver.

`CTPG_train()` returns the following outputs:

- `result`: The final result of parameter optimisation.
- `fwd_ensemble_sol`: The ensemble solution of forward simulation using the final neural network parameters.
- `loss_history`: The history of loss function evaluated at each iteration.
"""

function CTPG_train(dynamics_plant::Function, dynamics_controller::Function, dynamics_sensor::Function, cost_running::Function, cost_terminal::Function, cost_regularisor::Function, policy_NN, scenario; solve_alg = Tsit5(), sense_alg = InterpolatingAdjoint(autojacvec = ZygoteVJP()), ensemble_alg = EnsembleThreads(), opt_1 = ADAM(0.01), opt_2 = LBFGS(), maxiters_1 = 100, maxiters_2 = 100, progress_plot = true, i_nominal = nothing, p_NN_0 = nothing, solve_kwargs...)

    # scenario parameters
    @unpack ensemble, t_span, t_save, dim_x, dim_x_c = scenario
    dim_ensemble = length(ensemble)
    mean_factor  = Float32(1 / dim_ensemble)
    dim_t_save   = length(t_save)
    if isnothing(i_nominal)
        i_nominal = max(round(Int, dim_ensemble / 2), 1)
    end

    # NN parameters initialisation
    if isnothing(p_NN_0)
        p_NN_0 = initial_params(policy_NN)
    end

    # augmented dynamics
    function fwd_dynamics(r)
        return function (x_aug, p_NN, t)
            x   = x_aug[1:dim_x]
            x_c = x_aug[dim_x+1:end-1]
            # ∫cost_running = x_aug[end]

            y            = dynamics_sensor(t, x)
            (dx_c, u, _) = dynamics_controller(t, x_c, y, r, p_NN, policy_NN)
            dx           = dynamics_plant(t, x, u)

            return [dx; dx_c; cost_running(t, x, y, u, r)]
        end
    end

    # ODE problem construction
    prob_base = ODEProblem(fwd_dynamics(ensemble[i_nominal].r), [ensemble[i_nominal].x₀; zeros(Float32, dim_x_c + 1)], t_span, p_NN_0)

    function generate_probs(p_NN)
        return function (prob, i, repeat)
            remake(prob, f = fwd_dynamics(ensemble[i].r), u0 = [ensemble[i].x₀; zeros(Float32, dim_x_c + 1)], p = p_NN)
        end
    end

    # if dim_ensemble == 1
    #     prob_mtk = modelingtoolkitize(prob_base)
    #     prob_base = ODEProblem(prob_mtk, [], t_span, jac = true)
    #     ensemble_alg = EnsembleSerial()
    # end

    # loss function definition
    function loss(p_NN)
        ensemble_prob = EnsembleProblem(prob_base, prob_func = generate_probs(p_NN))

        fwd_ensemble_sol_full = solve(ensemble_prob, solve_alg, ensemble_alg, saveat = t_save, trajectories = dim_ensemble, sensealg = sense_alg; solve_kwargs...)

        # version 1: mean(sol[end])
        sol_length = Float32.([max(1,length(fwd_ensemble_sol_full[i])) for i in 1:dim_ensemble])
        sol_length_ratio = maximum(sol_length)./sol_length

        # if fwd_ensemble_sol_full[1].retcode == :Success
        if size(Array(fwd_ensemble_sol_full[1]), 2) > 10
            loss_val = mean([fwd_ensemble_sol_full[i][end][end] + cost_terminal(fwd_ensemble_sol_full[i][end][1:dim_x], ensemble[i].r) for i in 1:dim_ensemble] .* sol_length_ratio) + cost_regularisor(p_NN)
        else
            loss_val = Inf
        end

        # version 2: mean(Array[end])
        # loss_val = mean([(Array(fwd_ensemble_sol_full[i])[end,end] + cost_terminal(Array(fwd_ensemble_sol_full[i])[1:dim_x,end], ensemble[i].r)) for i in 1:dim_ensemble] .* sol_length_ratio) + cost_regularisor(p_NN)

        # version 3: scalar operation (best for code robustness as it can handle divergent case in ensemble)
        # loss_val = 0.0f0
        # for i in 1:dim_ensemble
        #     fwd_sol = Array(fwd_ensemble_sol_full[i])

        #     if size(fwd_sol,2) > dim_t_save
        #         fwd_sol = fwd_sol[:,1:dim_t_save]
        #     end

        #     if size(fwd_sol,2) == dim_t_save
        #         x_aug_f = fwd_sol[:,end]
        #         loss_val += (x_aug_f[end] + cost_terminal(x_aug_f[1:dim_x], ensemble[i].r)) * mean_factor
        #     else
        #         loss_val += 1000.0f0
        #     end
        # end
        # loss_val += cost_regularisor(p_NN)

        return loss_val, fwd_ensemble_sol_full
    end

    # learning progress callback setup
    loss_history = fill(NaN32, maxiters_1 + maxiters_2 + 2)
    iterator_learning = 1
    cb_progress = function (p_NN_val, loss_val, fwd_ensemble_sol_full; plot_val = progress_plot)
        @show (loss_val, iterator_learning);
        loss_history[iterator_learning] = loss_val
        if plot_val
            fwd_ensemble_sol = Array(fwd_ensemble_sol_full[i_nominal])
            display(scatter(fwd_ensemble_sol[1:dim_x, :]', label = :false, plot_title = "System State: Learning Iteration $(iterator_learning)", layout = (dim_x, 1), size = (700, 200 * dim_x)))
        end
        iterator_learning += 1
        return false
    end

    # NN training
    if maxiters_1 <= 0
        result_coarse = (; u = p_NN_0)
    else
        result_coarse = DiffEqFlux.sciml_train(loss, p_NN_0, opt_1; cb = cb_progress, maxiters = maxiters_1)
    end

    if maxiters_2 <= 0
        result = result_coarse
    else
        result = DiffEqFlux.sciml_train(loss, result_coarse.u, opt_2; cb = cb_progress, maxiters = maxiters_2)
    end

    # Forward solution for optimised p_NN
    function eval_IO(r, p_NN)
        # Evaluation of input u, output y, and NN output y_NN at each t
        return function (t, x_aug)
            x = x_aug[1:dim_x]
            x_c = x_aug[dim_x+1:end-1]

            y = dynamics_sensor(t, x)
            (_, u, y_NN) = dynamics_controller(t, x_c, y, r, p_NN, policy_NN)

            return (; u = u, y = y, y_NN = y_NN)
        end
    end

    function generate_savedata(p_NN)
        return function (sol, i)
            IO = eval_IO(ensemble[i].r, p_NN).(sol.t, sol.u)
            return ((; sol = sol, u = [u for (u, y, y_NN) in IO], y = [y for (u, y, y_NN) in IO], y_NN = [y_NN for (u, y, y_NN) in IO]), false)
        end
    end

    fwd_ensemble_sol = solve(EnsembleProblem(prob_base, prob_func = generate_probs(result.u), output_func = generate_savedata(result.u)), solve_alg, ensemble_alg, saveat = t_save, trajectories = dim_ensemble, sensealg = sense_alg; solve_kwargs...)

    loss_history = filter(!isnan, loss_history)

    return result, fwd_ensemble_sol, loss_history
end

function view_result(i_plot_list, fwd_ensemble_sol, loss_history, save_fig = true;
    vars_x = nothing, x_names = [], vars_u = nothing, u_names = [], vars_y = nothing, y_names = [], vars_y_NN = nothing, y_NN_names = [], plot_kwargs...)

    # labels, variables to be plotted, and layouts setting
    (dim_x_aug, dim_u, dim_y, dim_y_NN) = (length(fwd_ensemble_sol[1].sol[1]), length(fwd_ensemble_sol[1].u[1]), length(fwd_ensemble_sol[1].y[1]), length(fwd_ensemble_sol[1].y_NN[1]))
    (xlabel_t, ylabel_x, ylabel_u, ylabel_y, ylabel_y_NN) = ("\$t\$", "state", "input", "output", "NN output")

    if ~isnothing(vars_x)
        layout_x = (length(vars_x), 1)
        if ~isempty(x_names)
            ylabel_x = x_names[vars_x']
        end
    else
        layout_x = (dim_x_aug - 1, 1)
        vars_x = 1:dim_x_aug-1
    end

    if ~isnothing(vars_u)
        layout_u = (length(vars_u), 1)
        if ~isempty(u_names)
            ylabel_u = u_names[vars_u']
        end
    else
        layout_u = (dim_u, 1)
        vars_u = 1:dim_u
    end

    if ~isnothing(vars_y)
        layout_y = (length(vars_y), 1)
        if ~isempty(y_names)
            ylabel_y = y_names[vars_y']
        end
    else
        layout_y = (dim_y, 1)
        vars_y = 1:dim_y
    end

    if ~isnothing(vars_y_NN)
        layout_y_NN = (length(vars_y_NN), 1)
        if ~isempty(y_NN_names)
            ylabel_y_NN = y_NN_names[vars_y_NN']
        end
    else
        layout_y_NN = (dim_y_NN, 1)
        vars_y_NN = 1:dim_y_NN
    end

    # plotting
    if isempty(i_plot_list)
        i_plot_list = 1:length(fwd_ensemble_sol)
    end

    (f_x, f_u, f_y, f_y_NN) = (plot(), plot(), plot(), plot())
    for i in i_plot_list
        # data preprocessing
        @unpack sol, u, y, y_NN = fwd_ensemble_sol[i]
        @unpack t = sol

        u = reduce(vcat, u')     # hcat(u...)'
        y = reduce(vcat, y')     # hcat(y...)'
        y_NN = reduce(vcat, y_NN')  # hcat(y_NN...)'

        if i == first(i_plot_list)
            f_x = plot(sol, vars = vars_x, layout = layout_x, label = :false, xlabel = xlabel_t, ylabel = ylabel_x, size = (800, 160 * length(vars_x)); plot_kwargs...)

            f_u = plot(t, u[:, vars_u], layout = layout_u, label = :false, xlabel = xlabel_t, ylabel = ylabel_u, size = (800, 160 * length(vars_u)); plot_kwargs...)

            f_y = plot(t, y[:, vars_y], layout = layout_y, label = :false, xlabel = xlabel_t, ylabel = ylabel_y, size = (800, 160 * length(vars_y)); plot_kwargs...)

            f_y_NN = plot(t, y_NN[:, vars_y_NN], layout = layout_y_NN, label = :false, xlabel = xlabel_t, ylabel = ylabel_y_NN, size = (800, 160 * length(vars_y_NN)); plot_kwargs...)
        else
            plot!(f_x, sol, vars = vars_x, layout = layout_x, label = :false, xlabel = xlabel_t, ylabel = ylabel_x, size = (800, 160 * length(vars_x)); plot_kwargs...)

            plot!(f_u, t, u[:, vars_u], layout = layout_u, label = :false, xlabel = xlabel_t, ylabel = ylabel_u, size = (800, 160 * length(vars_u)); plot_kwargs...)

            plot!(f_y, t, y[:, vars_y], layout = layout_y, label = :false, xlabel = xlabel_t, ylabel = ylabel_y, size = (800, 160 * length(vars_y)); plot_kwargs...)

            plot!(f_y_NN, t, y_NN[:, vars_y_NN], layout = layout_y_NN, label = :false, xlabel = xlabel_t, ylabel = ylabel_y_NN, size = (800, 160 * length(vars_y_NN)); plot_kwargs...)
        end
    end
    f_L = plot(loss_history, label = :false, xlabel = "iteration", ylabel = "\$L\$"; plot_kwargs...)

    display(f_x)
    display(f_u)
    display(f_y)
    display(f_y_NN)
    display(f_L)

    if save_fig
        savefig(f_x,    "f_x.pdf")
        savefig(f_u,    "f_u.pdf")
        savefig(f_y,    "f_y.pdf")
        savefig(f_y_NN, "f_y_NN.pdf")
        savefig(f_L,    "f_L.pdf")
    end

    return f_x, f_u, f_y, f_y_NN, f_L
end


text:
module MLJEnsembles

using MLJModelInterface
import MLJModelInterface: predict, fit
import MLJBase # still needed for aggregating measures in oob-estimates of error
using Random
using CategoricalArrays
using CategoricalDistributions
using ComputationalResources
using Distributed
import Distributions
using ProgressMeter
import StatsBase

export EnsembleModel

const MMI = MLJModelInterface

include("ensembles.jl")

end # module


text:
function update_sigma!(state::State)
  state.sigma .= Util.scalefromaltskewt.(sqrt.(state.omega), state.psi)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
mutable struct bin_BB_Tree{BB_Type}
    id::Int64
    box::BB_Type
    node_1::bin_BB_Tree
    node_2::bin_BB_Tree
    function bin_BB_Tree(id::Int64, bb::BB_Type) where {BB_Type <: BoundingBox}
        return new{BB_Type}(id, bb)
    end
    function bin_BB_Tree(node_1::bin_BB_Tree{BB_Type}, node_2::bin_BB_Tree{BB_Type}) where {BB_Type <: BoundingBox}
        aabb = BB_Type(node_1.box, node_2.box)
        return new{BB_Type}(-9999, aabb, node_1, node_2)
    end
    function bin_BB_Tree(bb::BB_Type, node_1::bin_BB_Tree{BB_Type}, node_2::bin_BB_Tree{BB_Type}) where {BB_Type <: BoundingBox}
        return new{BB_Type}(-9999, bb, node_1, node_2)
    end
end

function Base.show(io::IO, bin_tree::bin_BB_Tree)
    n_leaf = leafNumber(bin_tree)
    n_depth = treeDepth(bin_tree)
    println(io, "bin_BB_Tree with $n_leaf leaves in $n_depth layers:")
    if n_depth == 1
        frac_high = n_leaf / (2^n_depth - 1)
        frac_low = n_leaf / (2^(n_depth - 1) - 1)
        println("    fill between $(frac_high) and $(frac_low) of optimal")
        if (frac_high < 0.1) || (frac_low < 0.2)
            @warn("tree may be unbalanced")
        end
    end
end

mutable struct TT_Cache
    vc::VectorCache{Tuple{Int64, Int64}}
    t_a_b::SVector{3,Float64}
    R_a_b::SMatrix{3,3,Float64,9}
    abs_R_a_b::SMatrix{3,3,Float64,9}
    function TT_Cache()
        vc = VectorCache{NTuple{2,Int64}}()
        return new(vc)
    end
end

update_TT_Cache!(tt::TT_Cache, trans::SVector{3,Float64}, R::RotMatrix{3,Float64,9}) = update_TT_Cache!(tt, trans, SMatrix{3,3,Float64,9}(R))
function update_TT_Cache!(tt::TT_Cache, trans::SVector{3,Float64}, R::SMatrix{3,3,Float64,9})
    tt.t_a_b = trans
    tt.R_a_b = R
    tt.abs_R_a_b = abs.(R) .+ 1.0e-14
    empty!(tt.vc)
    return nothing
end

Base.length(tt::TT_Cache) = length(tt.vc)
Base.isempty(tt::TT_Cache) = isempty(tt.vc)
Base.@propagate_inbounds Base.getindex(tt::TT_Cache, i::Int) = tt.vc[i]

is_leaf(tree::bin_BB_Tree)     = (tree.id != -9999)
is_not_leaf(tree::bin_BB_Tree) = (tree.id == -9999)

treeDepth(t::bin_BB_Tree) = leafNumberDepth(t, 0, 0)[1]
leafNumber(t::bin_BB_Tree) = leafNumberDepth(t, 0, 0)[2]

function leafNumberDepth(t::bin_BB_Tree, k_depth::Int64, k_leaf::Int64)
    if is_leaf(t)
        k_leaf = 1
    else is_not_leaf(t)
        k_depth_1, k_leaf_1 = leafNumberDepth(t.node_1, k_depth + 1, k_leaf)
        k_depth_2, k_leaf_2 = leafNumberDepth(t.node_2, k_depth + 1, k_leaf)
        k_depth = max(k_depth_1, k_depth_2)
        k_leaf = k_leaf_1 + k_leaf_2
    end
    return k_depth, k_leaf
end

function extractData!(tree::bin_BB_Tree, v::Vector{Int64})
    if is_leaf(tree)
        push!(v, tree.id)
    else
        extractData!(tree.node_1, v)
        extractData!(tree.node_2, v)
    end
end
function extractData(tree::bin_BB_Tree)
    v = Vector{Int64}()
    extractData!(tree, v)
    return v
end

function tree_tree_intersect(ttCache::TT_Cache, tree_1::bin_BB_Tree, tree_2::bin_BB_Tree)
    if BB_BB_intersect(ttCache, tree_1.box, tree_2.box)
        is_leaf_1 = is_leaf(tree_1)
        is_leaf_2 = is_leaf(tree_2)
        if is_leaf_1
            if is_leaf_2
                addCacheItem!(ttCache.vc, (tree_1.id, tree_2.id))
            else
                tree_tree_intersect(ttCache, tree_1, tree_2.node_1)
                tree_tree_intersect(ttCache, tree_1, tree_2.node_2)
            end
        else
            if is_leaf_2
                tree_tree_intersect(ttCache, tree_1.node_1, tree_2)
                tree_tree_intersect(ttCache, tree_1.node_2, tree_2)
            else
                tree_tree_intersect(ttCache, tree_1.node_1, tree_2.node_1)
                tree_tree_intersect(ttCache, tree_1.node_2, tree_2.node_1)
                tree_tree_intersect(ttCache, tree_1.node_1, tree_2.node_2)
                tree_tree_intersect(ttCache, tree_1.node_2, tree_2.node_2)
            end
        end
    end
end


text:
using Test

@testset "ME" begin
    include("ModelEnumeration.jl")
end

@testset "DB" begin
    include("DB.jl")
end

@testset "FLS" begin
    include("FLS.jl")
end


text:
using ForwardDiff, Distributions, FiniteDifferences, Tracker, Random, LinearAlgebra, PDMats
using Turing: Turing, gradient_logp_reverse, invlink, link, SampleFromPrior
using Turing.Core.RandomVariables: getval
using Turing.Core: TuringMvNormal, TuringDiagNormal
using ForwardDiff: Dual
using StatsFuns: binomlogpdf, logsumexp
using Test, LinearAlgebra
const FDM = FiniteDifferences

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

_to_cov(B) = B * B' + Matrix(I, size(B)...)
@testset "ad.jl" begin
    @turing_testset "adr" begin
        ad_test_f = gdemo_default
        vi = Turing.VarInfo()
        ad_test_f(vi, SampleFromPrior())
        svn = vi.vns[1]
        mvn = vi.vns[2]
        _s = getval(vi, svn)[1]
        _m = getval(vi, mvn)[1]

        x = map(x->Float64(x), vi[SampleFromPrior()])
        ∇E = gradient_logp_reverse(x, vi, ad_test_f)[2]
        grad_Turing = sort(∇E)

        dist_s = InverseGamma(2,3)

        # Hand-written logp
        function logp(x::Vector)
          s = x[2]
          # s = invlink(dist_s, s)
          m = x[1]
          lik_dist = Normal(m, sqrt(s))
          lp = logpdf(dist_s, s) + logpdf(Normal(0,sqrt(s)), m)
          lp += logpdf(lik_dist, 1.5) + logpdf(lik_dist, 2.0)
          lp
        end

        # Call ForwardDiff's AD
        g = x -> ForwardDiff.gradient(logp, x);
        # _s = link(dist_s, _s)
        _x = [_m, _s]
        grad_FWAD = sort(g(_x))

        # Compare result
        @test grad_Turing ≈ grad_FWAD atol=1e-9
    end
    @turing_testset "passing duals to distributions" begin
        float1 = 1.1
        float2 = 2.3
        d1 = Dual(1.1)
        d2 = Dual(2.3)

        @test logpdf(Normal(0, 1), d1).value ≈ logpdf(Normal(0, 1), float1) atol=0.001
        @test logpdf(Gamma(2, 3), d2).value ≈ logpdf(Gamma(2, 3), float2) atol=0.001
        @test logpdf(Beta(2, 3), (d2 - d1) / 2).value ≈ logpdf(Beta(2, 3), (float2 - float1) / 2) atol=0.001

        @test pdf(Normal(0, 1), d1).value ≈ pdf(Normal(0, 1), float1) atol=0.001
        @test pdf(Gamma(2, 3), d2).value ≈ pdf(Gamma(2, 3), float2) atol=0.001
        @test pdf(Beta(2, 3), (d2 - d1) / 2).value ≈ pdf(Beta(2, 3), (float2 - float1) / 2) atol=0.001
    end
    @turing_testset "general AD tests" begin
        # Tests gdemo gradient.
        function logp1(x::Vector)
            dist_s = InverseGamma(2, 3)
            s = x[2]
            m = x[1]
            lik_dist = Normal(m, sqrt(s))
            lp = Turing.logpdf_with_trans(dist_s, s, false) + Turing.logpdf_with_trans(Normal(0,sqrt(s)), m, false)
            lp += logpdf(lik_dist, 1.5) + logpdf(lik_dist, 2.0)
            return lp
        end

        test_model_ad(gdemo_default, logp1, [:m, :s])

        # Test Wishart AD.
        @model wishart_ad() = begin
            v ~ Wishart(7, [1 0.5; 0.5 1])
            v
        end

        # Hand-written logp
        function logp3(x)
            dist_v = Wishart(7, [1 0.5; 0.5 1])
            v = [x[1] x[3]; x[2] x[4]]
            lp = Turing.logpdf_with_trans(dist_v, v, false)
            return lp
        end

        test_model_ad(wishart_ad(), logp3, [:v])
    end
    @turing_testset "Tracker + logdet" begin
        rng, N = MersenneTwister(123456), 7
        ȳ, B = randn(rng), randn(rng, N, N)
        test_tracker_ad(B->logdet(cholesky(_to_cov(B))), ȳ, B; rtol=1e-8, atol=1e-8)
    end
    @turing_testset "Tracker + fill" begin
        rng = MersenneTwister(123456)
        test_tracker_ad(x->fill(x, 7), randn(rng, 7), randn(rng))
        test_tracker_ad(x->fill(x, 7, 11), randn(rng, 7, 11), randn(rng))
        test_tracker_ad(x->fill(x, 7, 11, 13), rand(rng, 7, 11, 13), randn(rng))
    end
    @turing_testset "Tracker + MvNormal" begin
        rng, N = MersenneTwister(123456), 11
        B = randn(rng, N, N)
        m, A = randn(rng, N), B' * B + I

        # Generate from the TuringMvNormal
        d, back = Tracker.forward(TuringMvNormal, m, A)
        x = Tracker.data(rand(d))

        # Check that the logpdf agrees with MvNormal.
        d_ref = MvNormal(m, PDMat(A))
        @test logpdf(d, x) ≈ logpdf(d_ref, x)

        test_tracker_ad((m, B, x)->logpdf(MvNormal(m, _to_cov(B)), x), randn(rng), m, B, x)
    end
    @turing_testset "Tracker + Diagonal Normal" begin
        rng, N = MersenneTwister(123456), 11
        m, σ = randn(rng, N), exp.(0.1 .* randn(rng, N)) .+ 1

        d = TuringDiagNormal(m, σ)
        x = rand(d)

        # Check that the logpdf agrees with MvNormal.
        d_ref = MvNormal(m, σ)
        @test logpdf(d, x) ≈ logpdf(d_ref, x)

        test_tracker_ad((m, σ, x)->logpdf(MvNormal(m, σ), x), randn(rng), m, σ, x)
    end
    @turing_testset "Tracker + MvNormal Interface" begin
        # Note that we only test methods where the `MvNormal` ctor actually constructs
        # a TuringMvNormal.

        rng, N = MersenneTwister(123456), 7
        m, b, B, x = randn(rng, N), randn(rng, N), randn(rng, N, N), randn(rng, N)
        ȳ = randn(rng)

        # zero mean, dense covariance
        test_tracker_ad((B, x)->logpdf(MvNormal(_to_cov(B)), x), randn(rng), B, x)
        test_tracker_ad(B->logpdf(MvNormal(_to_cov(B)), x), randn(rng), B)

        # zero mean, diagonal covariance
        test_tracker_ad((b, x)->logpdf(MvNormal(exp.(b)), x), randn(rng), b, x)
        test_tracker_ad(b->logpdf(MvNormal(exp.(b)), x), randn(rng), b)

        # dense mean, dense covariance
        test_tracker_ad((m, B, x)->logpdf(MvNormal(m, _to_cov(B)), x),
            randn(rng),
            randn(rng, N), randn(rng, N, N), randn(rng, N),
        )
        test_tracker_ad((m, B)->logpdf(MvNormal(m, _to_cov(B)), x),
            randn(rng),
            randn(rng, N), randn(rng, N, N),
        )
        test_tracker_ad((m, x)->logpdf(MvNormal(m, _to_cov(B)), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((B, x)->logpdf(MvNormal(m, _to_cov(B)), x),
            randn(rng),
            randn(rng, N, N), randn(rng, N),
        )
        test_tracker_ad(m->logpdf(MvNormal(m, _to_cov(B)), x), randn(rng), randn(rng, N))
        test_tracker_ad(B->logpdf(MvNormal(m, _to_cov(B)), x), randn(rng), randn(rng, N, N))

        # dense mean, diagonal covariance
        test_tracker_ad((m, b, x)->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N), randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((m, b)->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((m, x)->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((b, x)->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad(m->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N),
        )
        test_tracker_ad(b->logpdf(MvNormal(m, Diagonal(exp.(b))), x),
            randn(rng),
            randn(rng, N),
        )

        # dense mean, diagonal variance
        test_tracker_ad((m, b, x)->logpdf(MvNormal(m, exp.(b)), x),
            randn(rng),
            randn(rng, N), randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((m, b)->logpdf(MvNormal(m, exp.(b)), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((m, x)->logpdf(MvNormal(m, exp.(b)), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((b, x)->logpdf(MvNormal(m, exp.(b)), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad(m->logpdf(MvNormal(m, exp.(b)), x), randn(rng), randn(rng, N))
        test_tracker_ad(b->logpdf(MvNormal(m, exp.(b)), x), randn(rng), randn(rng, N))

        # dense mean, constant covariance
        b_s = randn(rng)
        test_tracker_ad((m, b, x)->logpdf(MvNormal(m, exp(b)), x),
            randn(rng),
            randn(rng, N), randn(rng), randn(rng, N),
        )
        test_tracker_ad((m, b)->logpdf(MvNormal(m, exp(b)), x),
            randn(rng),
            randn(rng, N), randn(rng),
        )
        test_tracker_ad((m, x)->logpdf(MvNormal(m, exp(b_s)), x),
            randn(rng),
            randn(rng, N), randn(rng, N)
        )
        test_tracker_ad((b, x)->logpdf(MvNormal(m, exp(b)), x),
            randn(rng),
            randn(rng), randn(rng, N),
        )
        test_tracker_ad(m->logpdf(MvNormal(m, exp(b_s)), x), randn(rng), randn(rng, N))
        test_tracker_ad(b->logpdf(MvNormal(m, exp(b)), x), randn(rng), randn(rng))

        # dense mean, constant variance
        b_s = randn(rng)
        test_tracker_ad((m, b, x)->logpdf(MvNormal(m, exp(b) * I), x),
            randn(rng),
            randn(rng, N), randn(rng), randn(rng, N),
        )
        test_tracker_ad((m, b)->logpdf(MvNormal(m, exp(b) * I), x),
            randn(rng),
            randn(rng, N), randn(rng),
        )
        test_tracker_ad((m, x)->logpdf(MvNormal(m, exp(b_s) * I), x),
            randn(rng),
            randn(rng, N), randn(rng, N),
        )
        test_tracker_ad((b, x)->logpdf(MvNormal(m, exp(b) * I), x),
            randn(rng),
            randn(rng), randn(rng, N),
        )
        test_tracker_ad(m->logpdf(MvNormal(m, exp(b_s) * I), x), randn(rng), randn(rng, N))
        test_tracker_ad(b->logpdf(MvNormal(m, exp(b) * I), x), randn(rng), randn(rng))

        # zero mean, constant variance
        test_tracker_ad((b, x)->logpdf(MvNormal(N, exp(b)), x),
            randn(rng),
            randn(rng), randn(rng, N),
        )
        test_tracker_ad(b->logpdf(MvNormal(N, exp(b)), x), randn(rng), randn(rng))
    end
    @testset "Simplex Tracker AD" begin
        @model dir() = begin
            theta ~ Dirichlet(1 ./ fill(4, 4))
        end
        Turing.setadbackend(:reverse_diff)
        sample(dir(), HMC(0.01, 1), 1000);
    end
    @testset "PDMatDistribution Tracker AD" begin
        @model wishart() = begin
            theta ~ Wishart(4, Matrix{Float64}(I, 4, 4))
        end
        Turing.setadbackend(:reverse_diff)
        sample(wishart(), HMC(0.01, 1), 1000);

        @model invwishart() = begin
            theta ~ InverseWishart(4, Matrix{Float64}(I, 4, 4))
        end
        Turing.setadbackend(:reverse_diff)
        sample(invwishart(), HMC(0.01, 1), 1000);
    end
end


text:
using Test
# using ProgressMeter

function constructSet(line)
    path = []
    pathWithSteps = []
    position = (0, 0)
    paths = (0, 0, 0)
    instructions = split(line, ",")
    totalSteps = 1
    for inst in instructions
        # println(inst)
        dir = inst[1]
        number = parse(Int, inst[2:end])
        for i in 1:number
            if dir == 'R'
                position = position[1], position[2] + 1
            elseif dir == 'L'
                position = position[1], position[2] - 1
            elseif dir == 'U'
                position = position[1] + 1, position[2]
            elseif dir == 'D'
                position = position[1] - 1, position[2]
            end
            append!(path, [position])
            append!(pathWithSteps, [(position[1], position[2], totalSteps)])
            totalSteps += 1
        end
    end
    return path, pathWithSteps
end

function getClosestIntersection(set1, set2)
    inter = intersect(set1, set2)
    best = minimum(pos -> abs(pos[1]) + abs(pos[2]), inter)
    return best
end

function part1(input)
    set1, _ = constructSet(input[1])
    set2, _ = constructSet(input[2])
    return getClosestIntersection(set1, set2)
end

function part2(input)
    path1, path1withSteps = constructSet(input[1])
    path2, path2withSteps = constructSet(input[2])
    inter = intersect(path1, path2)

    distances1 = Dict()
    # println(path1withSteps)
    for (x, y, steps) in path1withSteps
        if !haskey(distances1, (x, y))
            distances1[(x, y)] = steps
        end
    end

    distances2 = Dict()
    for (x, y, steps) in path2withSteps
        # println(x, y, steps)
        if !haskey(distances2, (x, y))
            distances2[(x, y)] = steps
        end
    end

    return minimum(pos -> distances1[pos] + distances2[pos], inter)
end

@test part2(["R75,D30,R83,U83,L12,D49,R71,U7,L72", "U62,R66,U55,R34,D71,R55,D58,R83"]) == 610
@test part2(["R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51", "U98,R91,D20,R16,D67,R40,U7,R15,U6,R7"]) == 410

# println("part1")
# println(part1(lines))
# println("part2")
# println(part2(lines))

function run(s)
    return part2(split(s, "\n"))
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
using QuadGK
export BlackBody,
    spectral_radiance,
    spectral_radiance_frequency,
    spectral_band_radiance,
    spectral_band_radiance_frequency,
    radiance,
    disk_nt_rel_factors,
    gravity_radius,
    disk_flux,
    compute_disk_temperature,
    compute_disk_luminosity,
    compute_disk_spectral_luminosity,
    compute_disk_spectral_flux,
    compute_disk_spectral_flux_at_radius,
    compute_disk_flux_at_radius

struct BlackBody
    T::Float64
end

kev_to_hz(energy) = energy / ERG_TO_KEV / H_PLANCK
lambda_peak(T) = 0.2897771955 / T
nu_peak(T) = C / lambda_peak(T)
nu_peak(bb::BlackBody) = nu_peak(bb.T)

"""
Spectral radiance of a black body bb for a frequency. Input frequency is
in Hz, and output units are erg / (s cm^2 sr Hz)
"""
function spectral_radiance_frequency(bb::BlackBody, frequency)
    a = 2 * H_PLANCK * frequency^3 / C^2
    b = 1 / (exp((H_PLANCK * frequency) / (K_B * bb.T)) - 1)
    return a * b
end

"""
Spectral radiance of a black body at a given energy. The input energy
has units of keV, and the result output is in units of
keV / (s cm² keV sr)
"""
function spectral_radiance(bb::BlackBody, energy)
    frequency = kev_to_hz(energy)
    rad = spectral_radiance_frequency(bb, frequency)
    return rad
end

"""
    spectral_band_radiance_frequency(bb::BlackBody, low, high)
Computes the radiance on a frequency band from low to high. low and high
have to be in Hz, return units are erg / sr / m^2
"""
function spectral_band_radiance_frequency(bb::BlackBody, low, high)
    low = max(low, nu_peak(bb) / 1e4)
    high = min(high, nu_peak(bb) * 1e4)
    integral, err = quadgk(x -> spectral_radiance_frequency(bb, x), low, high, rtol = 1e-8)
    return integral
end

"""
Computes the radiance on an energetic band from low to high.
"""
function spectral_band_radiance(bb::BlackBody, low, high)
    low = kev_to_hz(low)
    high = kev_to_hz(high)
    low = max(low, nu_peak(bb) / 1e3)
    high = min(high, nu_peak(bb) * 1e3)
    integral, err = quadgk(x -> spectral_radiance_frequency(bb, x), low, high, rtol = 1e-8)
    return integral
end

radiance(bb::BlackBody) = SIGMA_SB * bb.T^4 / π

"""
Novikov-Thorne relativistic factors for the AD spectrum.
If the passed radius is smaller than ISCO, it returns 0.

Parameters
----------
radius
    disc radius measured in R_g
spin
    normalized black hole spin (-1, 1)
isco
    innermost stable circular orbit measured in Rg
"""
function disk_nt_rel_factors(radius, spin, isco)
    if radius <= isco
        return 0.0
    end
    yms = sqrt(isco)
    y1 = 2 * cos((acos(spin) - pi) / 3)
    y2 = 2 * cos((acos(spin) + pi) / 3)
    y3 = -2 * cos(acos(spin) / 3)
    y = sqrt(radius)
    C = 1 - 3 / radius + 2 * spin / radius^1.5
    B = 3 * (y1 - spin)^2 * log((y - y1) / (yms - y1)) / (y * y1 * (y1 - y2) * (y1 - y3))
    B += 3 * (y2 - spin)^2 * log((y - y2) / (yms - y2)) / (y * y2 * (y2 - y1) * (y2 - y3))
    B += 3 * (y3 - spin)^2 * log((y - y3) / (yms - y3)) / (y * y3 * (y3 - y1) * (y3 - y2))
    A = 1 - yms / y - 3 * spin * log(y / yms) / (2 * y)
    factor = (A - B) / C
    return factor
end

disk_nt_rel_factors(bh, r) = disk_nt_rel_factors(r, bh.spin, bh.isco)

"""
Self-gravity radius as described by Laor & Netzer (1989).
"""
function gravity_radius(bh)
    mass = (bh.M / (1e9 * M_SUN))
    alpha = 0.1 # assumption
    r_sg = 2150 * mass^(-2 / 9) * bh.mdot^(4 / 9) * alpha^(2 / 9)# / bh.Rg
    return r_sg
end

"""
Computes the disk flux of the accretion disc of black hole ``bh``
at a given radius ``r`` (in Rg). See equation 5 of https://arxiv.org/pdf/2001.04720.pdf
"""
function disk_flux(bh, r)
    Mdot = compute_mass_accretion_rate(bh)
    NT_factors = disk_nt_rel_factors(bh, r)
    return 3 * G * bh.M * Mdot / (8 * π * (r * bh.Rg)^3) * NT_factors
end

"""
Computes the temperature of the disk at a radius r (in units of Rg).
"""
function compute_disk_temperature(bh, r)
    flux = disk_flux(bh, r)
    return (flux / SIGMA_SB)^(1 / 4)
end

"""
    compute_disk_luminosity(bh, r_min, r_max)
computes the luminosity of the disk between r_min and r_max
"""
function compute_disk_luminosity(bh, r_min, r_max)
    constant = SIGMA_SB * 4 * π * bh.Rg^2
    integ, _ = quadgk(
        r -> r * compute_disk_temperature(bh, r)^4,
        r_min,
        r_max,
        atol = 0,
        rtol = 1e-2,
    )
    return integ * constant
end

"""
    compute_disk_spectral_flux_at_radius(bh, radius, energy_range)
Returns an array with the photon flux at each bin. (# photons / cm^2)
"""
function compute_disk_spectral_flux_at_radius(bh, radius, energy)
    ret = Float64[]
    bb = BlackBody(compute_disk_temperature(bh, radius))
    freq = kev_to_hz(energy)
    radiance = spectral_radiance_frequency(bb, freq)
    return radiance * π / HZ_TO_KEV
end

"""
    compute_disk_flux_at_radius(bh, radius, energy_range)
Returns an array with the photon flux at each bin. (# photons / cm^2)
"""
function compute_disk_flux_at_radius(bh, radius, low = 1e-5, high = 1e2)
    ret = Float64[]
    bb = BlackBody(compute_disk_temperature(bh, radius))
    radiance = spectral_band_radiance(bb, low, high)
    return radiance * π
end

function compute_disk_spectral_luminosity(bh, energy_range; r_min=6, r_max=1600)
    ret = Float64[]
    for energy in energy_range
        integ, _ = quadgk(
            r -> r * compute_disk_spectral_flux_at_radius(bh, r, energy),
            r_min,
            r_max,
            rtol = 1e-8,
            atol = 0,
        )
        push!(ret, integ)
    end
    return ret * 2 * 2π * bh.Rg^2
end