text:
# Data structures and methods for actin filaments
# Alex Tam, 12/10/2020

"Data structure for actin filament segments"
mutable struct Segment
    index::Int # Integer segment index
    L_eq::Float64 # [μm] Segment equilibrium length
    t::Vector{Vector{Int}} # Translations required by periodic BCs
end

"Mutable data structure for actin filaments"
mutable struct Actin_Filament
    index::Int # Integer filament index
    segments::Vector{Segment} # List of segments on actin filament
    "Inner constructor to generate actin filaments and update State"
    function Actin_Filament(s, ind, parA, Lxx, Lyy)
        # 1. Generate actin filament segments
        segments = Vector{Segment}(undef, parA.nSeg); # Pre-allocate
        for i = 1:parA.nSeg
            segments[i] = Segment(i, parA.LSeg, []); # Create segments
        end
        # 2. Generate centre position and orientation
        angle = pi/4 + (ind-1)*pi/2; # Prescribe actin filament orientation
        # 3. Calculate minus end position
        nodes = Vector{}(undef, parA.nSeg+1); # Pre-allocate
        nx = 0.5; ny = 0.5; # Dimensionless actin filament minus end position
        nodes[1] = [nx, ny]; # Store minus end position
        # 4. Calculate node positions
        for i = 1:parA.nSeg
            nx += segments[i].L_eq*cos(angle)/Lxx;
            ny += segments[i].L_eq*sin(angle)/Lyy;
            nodes[i+1] = [nx, ny];
        end
        push!(s.an, nodes); # Push nodes to State
        # 5. Apply periodic BCs to each segment
        for i = 1:parA.nSeg
            segments[i].t = periodic(s.an[ind][i+1], s.an[ind][i]); # Store required translations
        end
        return new(ind, segments)
    end
end

"Obtain dimensional segment lengths"
function get_segment_lengths(f::Actin_Filament, s::State{T}, Lxx, Lxy, Lyx, Lyy) where {T}
    Ls = Vector{T}(undef, length(f.segments)); # Pre-allocate segment lengths
    for seg in f.segments
        mx, my, px, py = get_segment_nodes(f, seg, s, Lxx, Lxy, Lyx, Lyy); # Obtain segment nodes (dimensional, un-translated)
        Ls[seg.index] = sqrt((px - mx)^2 + (py - my)^2); # Compute segment length
    end
    return Ls
end

"Obtain dimensional length of a particular segment"
function get_segment_length(f::Actin_Filament, s::State{T}, seg::Segment, Lxx, Lxy, Lyx, Lyy) where {T}
    mx, my, px, py = get_segment_nodes(f, seg, s, Lxx, Lxy, Lyx, Lyy); # Obtain segment nodes (dimensional, un-translated)
    Ls = sqrt((px - mx)^2 + (py - my)^2); # Compute segment length
    return Ls
end

"Obtain dimensional length of a particular segment"
function get_segment_length(f::Actin_Filament, s::State{T}, seg::Segment, Lxx, Lxy, Lyx, Lyy) where {T}
    mx, my, px, py = get_segment_nodes(f, seg, s, Lxx, Lxy, Lyx, Lyy); # Obtain segment nodes (dimensional, un-translated)
    Ls = sqrt((px - mx)^2 + (py - my)^2); # Compute segment length
    return Ls
end

"Obtain dimensionless, un-translated node positions"
function get_segment_nodes(f::Actin_Filament, seg::Segment, s::State{T}) where {T}
    mx = s.an[f.index][seg.index][1];
    my = s.an[f.index][seg.index][2];
    px = s.an[f.index][seg.index+1][1];
    py = s.an[f.index][seg.index+1][2];
    return mx, my, px, py
end

"Obtain dimensionless, translated node positions"
function get_segment_nodes(f::Actin_Filament, seg::Segment, s::State{T}, t::Vector{Int}) where {T}
    mx = s.an[f.index][seg.index][1] - t[1];
    my = s.an[f.index][seg.index][2] - t[2];
    px = s.an[f.index][seg.index+1][1] - t[1];
    py = s.an[f.index][seg.index+1][2] - t[2];
    return mx, my, px, py
end

"Obtain dimensional, un-translated node positions"
function get_segment_nodes(f::Actin_Filament, seg::Segment, s::State{T}, Lxx, Lxy, Lyx, Lyy) where {T}
    mx = s.an[f.index][seg.index][1]*Lxx + s.an[f.index][seg.index][2]*Lyx;
    my = s.an[f.index][seg.index][1]*Lxy + s.an[f.index][seg.index][2]*Lyy;
    px = s.an[f.index][seg.index+1][1]*Lxx + s.an[f.index][seg.index+1][2]*Lyx;
    py = s.an[f.index][seg.index+1][1]*Lxy + s.an[f.index][seg.index+1][2]*Lyy;
    return mx, my, px, py
end

"Obtain dimensional, translated node positions"
function get_segment_nodes(f::Actin_Filament, seg::Segment, s::State{T}, t::Vector{Int}, Lxx, Lxy, Lyx, Lyy) where {T}
    mx = (s.an[f.index][seg.index][1] - t[1])*Lxx + (s.an[f.index][seg.index][2] - t[2])*Lyx;
    my = (s.an[f.index][seg.index][1] - t[1])*Lxy + (s.an[f.index][seg.index][2] - t[2])*Lyy;
    px = (s.an[f.index][seg.index+1][1] - t[1])*Lxx + (s.an[f.index][seg.index+1][2] - t[2])*Lyx;
    py = (s.an[f.index][seg.index+1][1] - t[1])*Lxy + (s.an[f.index][seg.index+1][2] - t[2])*Lyy;
    return mx, my, px, py
end


text:
using TinyG

TinyG.download_firmware()

ext: jl


text:
@testset "Scaling and Distance Matrices" begin

    #### L2 = Euclidean
        @testset "Euclidean Distance Matrix" begin
            objs = [[1,1,1], [2,3,4]]
            expected = [0. sqrt(14) ; sqrt(14) 0.]
            A = permutedims(hcat(objs...)) # makes 2x3 array
            @test pairwise(Euclidean(), A, dims = 1) ≈ expected
        end

        @testset "Euclidean simple" begin
            objs = [[1,1,1], [2,3,4]]
            expected = float(sqrt(14))
            @test Euclidean()(objs...) ≈ expected
        end

    #### KL Divergence
        @testset "Same as scipy.stats.entropy" begin
            u,v = [1/2, 1/2], [9/10, 1/10]
            expected = 0.5108256237659907
            @test KLDivergence()(u,v) ≈ expected
        end

    ##### EMD or 1-d Wasserstein

        # EMD or 1-d Wasserstein
        @testset "EMD No Fail" begin
            objs = [[1,1,1], [2,3,4]]
            @test EMD(objs[1], objs[2]) !== nothing
        end

        @testset "EMD with grid provided, same result as python" begin
            objs = [[3.4, 3.9, 7.5, 7.8], [4.5, 1.4],
            [1.4, 0.9, 3.1, 7.2], [3.2, 3.5]]
            expected = 4.0781331438047861
            @test EMD(objs...) ≈ expected
        end

        @testset "EMD grid provided, different constructor, no fail" begin
            g = [0.5, 1.0, 1.5, 2.0]
            u = [3.4, 3.9, 7.5, 7.8]
            v = [1.4, 0.9, 3.1, 7.2]
            @test EMD((g,g))(u,v) !== nothing
        end

        @testset "EMD with grid, different constructor, same result as python" begin
            g1 = [3.4, 3.9, 7.5, 7.8]
            g2 = [4.5, 1.4]
            u = [1.4, 0.9, 3.1, 7.2]
            v = [3.2, 3.5]
            expected = 4.0781331438047861
            @test EMD((g1,g2))(u,v) ≈ expected
        end

    ##### Energy

        @testset "Energy() no fail" begin
            objs = [[1,1,1], [2,3,4]]
            @test Energy(objs[1], objs[2]) !== nothing
        end

        @testset "energy() not null" begin
            objs = [[3.4, 3.9, 7.5, 7.8], [4.5, 1.4],
            [1.4, 0.9, 3.1, 7.2], [3.2, 3.5]]
            expected = 4.0781331438047861
            @test energy(objs...) !== nothing
        end

        @testset "Energy() produces same result as python, simple." begin
            objs = ([0, 8], [0, 8], [3, 1], [2, 2])
            expected = 1.0
            @test Energy()(objs...) ≈ expected
        end

        @testset "Energy() produces same result as python, harder." begin
            objs = [[3.4, 3.9, 7.5, 7.8], [4.5, 1.4],
            [1.4, 0.9, 3.1, 7.2], [3.2, 3.5]]
            expected = 2.3674181638546394
            @test Energy()(objs...) ≈ expected
        end

        @testset "energy() w/ different u and v. Compare to python result." begin
            objs = ([0.7, 7.4, 2.4, 6.8], [1.4, 8. ],[2.1, 4.2, 7.4, 8. ], [7.6, 8.8])
            expected = Float32(0.8800334097615822) # from scipy.stats.energy_distance
            @test energy(objs...) ≈ expected
        end

        @testset "energy() w/ different u and v. Compare to python result. Float32 math." begin
            objs = (Float32[0.7, 7.4, 2.4, 6.8], Float32[1.4, 8. ],Float32[2.1, 4.2, 7.4, 8. ], Float32[7.6, 8.8])
            expected = 0.8800334097615822 # from scipy.stats.energy_distance
            @test Float32(energy(objs...)) ≈ expected
        end

        @testset "energy((g1,g2)) w/ different u and v. " begin
            g1 = [0.7, 7.4, 2.4, 6.8]
            g2 = [1.4, 8. ]
            u = [2.1, 4.2, 7.4, 8. ]
            v = [7.6, 8.8]
            expected = 0.8800334097615822 # from scipy.stats.energy_distance
            @test Energy((g1,g2))(u,v) ≈ expected
        end

end


text:
"""
    PoissonProcess{D,M}

Homogeneous temporal Poisson process with arbitrary mark distribution.

# Fields
- `λ::Float64`: event rate.
- `mark_dist::D`: mark distribution.
"""
struct PoissonProcess{D,M} <: TemporalPointProcess{M}
    λ::Float64
    mark_dist::D
end

function PoissonProcess(λ::Float64, mark_dist::D) where {D}
    return PoissonProcess{D,sampletype(mark_dist)}(λ, mark_dist)
end

## Access

function mark_distribution(pp::PoissonProcess)
    return pp.mark_dist
end

function intensity(pp::PoissonProcess, m = nothing)::Float64
    if isnothing(m)
        return pp.λ
    else
        return pp.λ * density(mark_distribution(pp), m)
    end
end

## Simulation

function Base.rand(rng::AbstractRNG, pp::PoissonProcess, tmin::Real, tmax::Real)
    N = rand(rng, Dists.Poisson(intensity(pp) * (tmax - tmin)))
    times = rand(rng, Dists.Uniform(tmin, tmax), N)
    marks = [rand(rng, mark_distribution(pp)) for _ = 1:N]
    return TemporalHistory(times, marks, tmin, tmax)
end

Base.rand(pp::PoissonProcess, tmin::Real, tmax::Real) = rand(GLOBAL_RNG, pp, tmin, tmax)

## Likelihood

function MeasureTheory.logdensity(pp::PoissonProcess, h::TemporalHistory)
    l = -intensity(pp) * duration(h)
    mark_dist = mark_distribution(pp)
    for m in event_marks(h)
        l += logdensity(mark_dist, m)
    end
    return l
end

## Prior

struct PoissonProcessPrior{D}
    λα::Float64
    λβ::Float64
    mark_dist_prior::D
end

## Prior likelihood

function MeasureTheory.logdensity(prior::PoissonProcessPrior, pp::PoissonProcess)
    lλ = logdensity(Dists.Gamma(prior.λα, 1 / prior.λβ), pp.λ)
    ld = logdensity(prior.mark_dist_prior, pp.mark_dist)
    return lλ + ld
end

## Fitting

function Dists.fit(::Type{<:PoissonProcess{D}}, history::TemporalHistory) where {D}
    λ = nb_events(history) / duration(history)
    mark_dist = fit(D, history.marks)
    return PoissonProcess(λ, mark_dist)
end

function Dists.fit(
    ::Type{<:PoissonProcess{D}},
    prior::PoissonProcessPrior,
    history::TemporalHistory,
) where {D}
    λ = (nb_events(history) + prior.λα - 1) / (duration(history) + prior.λβ)
    mark_dist = fit(D, prior.mark_dist_prior, history.marks)
    return PoissonProcess(λ, mark_dist)
end


text:
include("twoBeams2D.jl")
N = 20
xhats = Array{Float64}(12,1000,N,2)
x = 0
first = true
p1 = plot(title="Improved, N = 3000")
p2 = plot(title="Improved, N = 3000")
for i = 1:N
  x, xhat = runTest(3000, pf, twoBeams2DSettings(true), 1000)
  xhats[:,:,i,1] = xhat
  if first
    plot!(p1, x[1,:]', x[2,:]', width=3)
    plot!(p2, x[5:6,:]', width=3)
  end
  first = false
  plot!(p1, xhat[1,:]', xhat[2,:]')
  plot!(p2, xhat[5:6,:]')
  gui(p1)
  gui(p2)
  println(i)
end

first = true
p3 = plot(title="Old, N = 3000")
p4 = plot(title="Old, N = 3000")
for i = 1:N
  x, xhat = runTest(4096, pf, twoBeams2DSettings(false), 1000)
  xhats[:,:,i,2] = xhat
  if first
    plot!(p3, x[1,:]', x[2,:]', width=3)
    plot!(p4, x[5:6,:]', width=3)
  end
  first = false
  plot!(p3, xhat[1,:]', xhat[2,:]')
  plot!(p4, xhat[5:6,:]')
  gui(p3)
  gui(p4)
  println(i)
end


text:
#################################################################################
# Type Alias for long type signatures
const MinMax = NamedTuple{(:min, :max), NTuple{2, Float64}}
const NamedMinMax = Tuple{String, MinMax}
const UpDown = NamedTuple{(:up, :down), NTuple{2, Float64}}
const InOut = NamedTuple{(:in, :out), NTuple{2, Float64}}
const StartUpStages = NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}

const BUILD_PROBLEMS_TIMER = TimerOutputs.TimerOutput()
const RUN_OPERATION_MODEL_TIMER = TimerOutputs.TimerOutput()
const RUN_SIMULATION_TIMER = TimerOutputs.TimerOutput()

# Type Alias for JuMP and PJ containers
const PGAE = PJ.ParametrizedGenericAffExpr{Float64, JuMP.VariableRef}
const GAE = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}
const JuMPAffineExpressionArray = Matrix{GAE}
const JuMPAffineExpressionVector = Vector{GAE}
const JuMPConstraintArray = DenseAxisArray{JuMP.ConstraintRef}
const JuMPVariableArray = DenseAxisArray{JuMP.VariableRef}
const JuMPParamArray = DenseAxisArray{PJ.ParameterRef}
const JuMPFloatArray = DenseAxisArray{Float64}

# Settings constants
const UNSET_HORIZON = 0
const UNSET_INI_TIME = Dates.DateTime(0)

# Tolerance of comparisons
# MIP gap tolerances in most solvers are set to 1e-4
const ABSOLUTE_TOLERANCE = 1.0e-3
const BALANCE_SLACK_COST = 1e6
const SERVICES_SLACK_COST = 1e5
const COST_EPSILON = 1e-3
const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0
const SECONDS_IN_MINUTE = 60.0
const MINUTES_IN_HOUR = 60.0
const SECONDS_IN_HOUR = 3600.0
const MAX_START_STAGES = 3
const OBJECTIVE_FUNCTION_POSITIVE = 1.0
const OBJECTIVE_FUNCTION_NEGATIVE = -1.0
const INITIALIZATION_PROBLEM_HORIZON = 3
# The DEFAULT_RESERVE_COST value is used to avoid degeneracy of the solutions, reseve cost isn't provided.
const DEFAULT_RESERVE_COST = 1.0
const KiB = 1024
const MiB = KiB * KiB
const GiB = MiB * KiB

const UNSUPPORTED_POWERMODELS =
    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]

const PSI_NAME_DELIMITER = "__"

const M_VALUE = 1e6

const NO_SERVICE_NAME_PROVIDED = ""
const CONTAINER_KEY_EMPTY_META = ""
const UPPER_BOUND = "ub"
const LOWER_BOUND = "lb"

# File Names definitions
const PROBLEM_SERIALIZATION_FILENAME = "operation_problem.bin"
const PROBLEM_LOG_FILENAME = "operation_problem.log"
const HASH_FILENAME = "check.sha256"
const SIMULATION_SERIALIZATION_FILENAME = "simulation.bin"
const SIMULATION_LOG_FILENAME = "simulation.log"
const REQUIRED_RECORDERS = (:simulation_status, :execution)
const KNOWN_SIMULATION_PATHS = [
    "data_store",
    "logs",
    "models_json",
    "problems",
    "recorder",
    "results",
    "simulation_files",
]
const RESULTS_DIR = "results"

# Enums
IS.@scoped_enum(BuildStatus, IN_PROGRESS = -1, BUILT = 0, FAILED = 1, EMPTY = 2,)
IS.@scoped_enum(
    RunStatus,
    NOT_READY = -2,
    READY = -1,
    SUCCESSFUL = 0,
    RUNNING = 1,
    FAILED = 2,
)
IS.@scoped_enum(SOSStatusVariable, NO_VARIABLE = 1, PARAMETER = 2, VARIABLE = 3,)

const ENUMS = (BuildStatus, RunStatus, SOSStatusVariable)

const ENUM_MAPPINGS = Dict()

for enum in ENUMS
    ENUM_MAPPINGS[enum] = Dict()
    for value in instances(enum)
        ENUM_MAPPINGS[enum][lowercase(string(value))] = value
    end
end

"""
Get the enum value for the string. Case insensitive.
"""
function get_enum_value(enum, value::String)
    if !haskey(ENUM_MAPPINGS, enum)
        throw(ArgumentError("enum=$enum is not valid"))
    end

    val = lowercase(value)
    if !haskey(ENUM_MAPPINGS[enum], val)
        throw(ArgumentError("enum=$enum does not have value=$val"))
    end

    return ENUM_MAPPINGS[enum][val]
end

Base.convert(::Type{BuildStatus}, val::String) = get_enum_value(BuildStatus, val)
Base.convert(::Type{RunStatus}, val::String) = get_enum_value(RunStatus, val)
Base.convert(::Type{SOSStatusVariable}, x::String) = get_enum_value(SOSStatusVariable, x)

# Store const definitions
# Update src/simulation/simulation_store_common.jl with any changes.
const STORE_CONTAINER_DUALS = :duals
const STORE_CONTAINER_PARAMETERS = :parameters
const STORE_CONTAINER_VARIABLES = :variables
const STORE_CONTAINER_AUX_VARIABLES = :aux_variables
const STORE_CONTAINER_EXPRESSIONS = :expressions
const STORE_CONTAINERS = (
    STORE_CONTAINER_DUALS,
    STORE_CONTAINER_PARAMETERS,
    STORE_CONTAINER_VARIABLES,
    STORE_CONTAINER_AUX_VARIABLES,
    STORE_CONTAINER_EXPRESSIONS,
)


text:
# Parameterize by T so that way it can be Vector{Expression} which is defined after
struct Operation <: AbstractOperation
    op::Function
    args::Vector{Expression}
end

# Recursive ==
function Base.:(==)(x::Operation,y::Operation)
    x.op == y.op && length(x.args) == length(y.args) && all(isequal.(x.args,y.args))
end
Base.:(==)(::Operation, ::Number   ) = false
Base.:(==)(::Number   , ::Operation) = false
Base.:(==)(::Operation, ::Variable ) = false
Base.:(==)(::Variable , ::Operation) = false
Base.:(==)(::Operation, ::Constant ) = false
Base.:(==)(::Constant , ::Operation) = false

Base.convert(::Type{Expr}, O::Operation) =
    build_expr(:call, Any[Symbol(O.op); convert.(Expr, O.args)])
Base.show(io::IO, O::Operation) = print(io, convert(Expr, O))

"""
find_replace(O::Operation, x::Expression, y::Expression)

Finds the expression `x` in Operation `O` and replaces it with the Expression `y`
"""
function find_replace!(O::Operation, x::Expression, y::Expression)
    for i in eachindex(O.args)
        if isequal(O.args[i], x)
            O.args[i] = y
        elseif typeof(O.args[i]) <: Operation
            find_replace!(O.args[i],x,y)
        end
    end
end

# For inv
Base.convert(::Type{Operation}, x::Number) = Operation(identity, Expression[Constant(x)])
Base.convert(::Type{Operation}, x::Operation) = x
Base.convert(::Type{Operation}, x::Expression) = Operation(identity, Expression[x])
Operation(x) = convert(Operation, x)


text:
abstract type StandardBasisPolynomial{T} <: AbstractPolynomial{T} end

const ⟒ = constructorof #\upin

function showterm(io::IO, ::Type{<:StandardBasisPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T}
    if iszero(pj) return false end
    pj = printsign(io, pj, first, mimetype)
    if !(pj == one(T) && !(showone(T) || j == 0))
        printcoefficient(io, pj, j, mimetype)
    end
    printproductsign(io, pj, j, mimetype)
    printexponent(io, var, j, mimetype)
    return true
end

# allows  broadcast  issue #209
evalpoly(x, p::StandardBasisPolynomial) = p(x)

domain(::Type{<:StandardBasisPolynomial}) = Interval(-Inf, Inf)
mapdomain(::Type{<:StandardBasisPolynomial}, x::AbstractArray) = x

## generic test if polynomial `p` is a constant
isconstant(p::StandardBasisPolynomial) = degree(p) <= 0

Base.convert(P::Type{<:StandardBasisPolynomial}, q::StandardBasisPolynomial) = isa(q, P) ? q : P([q[i] for i in 0:degree(q)], q.var)

function fromroots(P::Type{<:StandardBasisPolynomial}, r::AbstractVector{T}; var::SymbolLike = :x) where {T <: Number}
    n = length(r)
    c = zeros(T, n + 1)
    c[1] = one(T)
    for j in 1:n, i in j:-1:1
        c[(i + 1)] = c[(i + 1)] - r[j] * c[i]
    end
    #return P(c, var)
    return P(reverse(c), var)
end

function Base.:+(p::P, c::S) where {T, P <: StandardBasisPolynomial{T}, S<:Number}
    R = promote_type(T,S)
    as = R[c  for c in coeffs(p)]
    as[1] += c
    ⟒(P)(as, p.var)
end

function derivative(p::P, order::Integer = 1) where {P <: StandardBasisPolynomial}
    order < 0 && error("Order of derivative must be non-negative")
    T = eltype(p)
    # we avoid usage like Base.promote_op(*, T, Int) here, say, as
    # Base.promote_op(*, Rational, Int) is Any, not Rational in analogy to
    # Base.promote_op(*, Complex, Int)
    R = eltype(one(T)*1)
    order == 0 && return p
    hasnan(p) && return P(R[NaN], p.var)
    order > length(p) && return P(R[0], p.var)
    d = degree(p)
    d <= 0 && return zero(p)
    n = d + 1
    a2 = Vector{R}(undef, n - order)
    @inbounds for i in order:n - 1
        a2[i - order + 1] = reduce(*, (i - order + 1):i, init = p[i])
    end
    return ⟒(P)(a2, p.var)
end

function integrate(p::P, k::S) where {P <: StandardBasisPolynomial, S<:Number}
    T = eltype(p)
    R = eltype((one(T)+one(S))/1)
    if hasnan(p) || isnan(k)
        return ⟒(P)([NaN])
    end
    n = length(p)
    a2 = Vector{R}(undef, n + 1)
    a2[1] = k
    @inbounds for i in 1:n
        a2[i + 1] = p[i - 1] / i
    end
    return ⟒(P)(a2, p.var)
end

function Base.divrem(num::P, den::Q) where {P <: StandardBasisPolynomial, Q <: StandardBasisPolynomial}
    num.var != den.var && error("Polynomials must have same variable")
    var = num.var

    n = degree(num)
    m = degree(den)

    m == -1 && throw(DivideError())
    if m == 0 && den[0] ≈ 0 throw(DivideError()) end

    T, S =  eltype(num), eltype(den)
    R = eltype(one(T)/one(S))

    deg = n - m + 1

    if deg ≤ 0
        return zero(P, var), num
    end

    q_coeff = zeros(R, deg)
    r_coeff = R[ num[i-1] for i in 1:n+1 ]

    @inbounds for i in n:-1:m
        q = r_coeff[i + 1] / den[m]
        q_coeff[i - m + 1] = q
        @inbounds for j in 0:m
            elem = den[j] * q
            r_coeff[i - m + j + 1] -= elem
        end
    end

    return ⟒(P)(q_coeff, var), ⟒(P)(r_coeff, var)

end

function companion(p::P) where {P <: StandardBasisPolynomial}
    d = length(p) - 1
    d < 1 && error("Series must have degree greater than 1")
    d == 1 && return diagm(0 => [-p[0] / p[1]])

    T = eltype(p)
    R = eltype(one(T)/one(T))

    comp = diagm(-1 => ones(R, d - 1))
    ani = 1 / p[end]
    for j in  0:(degree(p)-1)
        comp[1,(d-j)] = -p[j] * ani # along top row has smaller residual than down column
    end
    return comp
end

function  roots(p::P; kwargs...)  where  {P <: StandardBasisPolynomial}
    T = eltype(p)
    R = eltype(one(T)/one(T))
    d = degree(p)
    if d < 1
        return []
    end
    d == 1 && return R[-p[0] / p[1]]

    as = [p[i] for i in 0:d]
    K  = findlast(!iszero, as)
    if K == nothing
        return R[]
    end
    k =  findfirst(!iszero, as)

    k  == K && return zeros(R, k-1)

    # find eigenvalues of the  companion matrix
    comp  = companion(⟒(P)(as[k:K], p.var))
    L = eigvals(comp; kwargs...)
    append!(L, zeros(eltype(L), k-1))

    L
end

function vander(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, n::Integer) where {T <: Number}
    A = Matrix{T}(undef, length(x), n + 1)
    A[:, 1] .= one(T)
    @inbounds for i in 1:n
        A[:, i + 1] = A[:, i] .* x
    end
    return A
end


text:
function CompositeDataContext(arg0::TimeScales, arg1::Frames, arg2::CelestialBodies, arg3::GravityFields, arg4::GeoMagneticFields)
    return CompositeDataContext((TimeScales, Frames, CelestialBodies, GravityFields, GeoMagneticFields), arg0, arg1, arg2, arg3, arg4)
end

function get_celestial_bodies(obj::CompositeDataContext)
    return jcall(obj, "getCelestialBodies", CelestialBodies, ())
end

function get_frames(obj::CompositeDataContext)
    return jcall(obj, "getFrames", Frames, ())
end

function get_geo_magnetic_fields(obj::CompositeDataContext)
    return jcall(obj, "getGeoMagneticFields", GeoMagneticFields, ())
end

function get_gravity_fields(obj::CompositeDataContext)
    return jcall(obj, "getGravityFields", GravityFields, ())
end

function get_time_scales(obj::CompositeDataContext)
    return jcall(obj, "getTimeScales", TimeScales, ())
end


text:
export
  evolve_operator,
  evolve

"""
    evolve_operator(evo_gen, time)

Return an exponent of `time`×`evo_gen`. This function is useful in the case of
multiple initial states and fixed `evo_gen` and `time`, as it is faster to
compute the exponent once and use it for evolving on different initial states.

*Note:* Parameter `evo_gen` must be of type `Matrix`. For type `SparseMatrixCSC`
case different numerical approach is used. See function ```epmv``` in package
```Expokit```.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
([0 1; 1 0], Array{Int64,2}[[0 1; 0 0], [0 0; 1 0]])

julia> evolve_operator(evolve_generator(H, L), 4.0)
4×4 Array{Complex{Float64},2}:
 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im
```
"""
function evolve_operator(evo_gen::AbstractMatrix{<:Number}, time::Real)
  @argumentcheck time>= 0 "Time has to be nonnegative"
  @argumentcheck size(evo_gen, 1) == size(evo_gen, 2) "Argument evo_gen has to be a square matrix"

  exp(time*evo_gen)
end

"""

    evolve(evogen, state, time)
    evolve(evogen, state, tpoints)
    evolve(evosuper, state)

Simulate the GKSL master equation according to the equation

``|result⟩⟩ = exp(time*evogen)|state⟩⟩``

where ``|⋅⟩⟩`` denotes the vectorization.

*Note:* The function returns unvectorized `result`.

The evolution can be calculated using three different approaches.

In the simplest case the function accepts matrix `evogen` specifying the
generator of the evolution, `state` describing the starting point of the
evolution, and `time` specifying the time of the evolution.

*Note:* If `evogen` is of type `Matrix`, the exponent is calculated using
`exp` function. If `evogen` is of type `SparseMatrixCSC`, `expmv` from
`Expokit.jl` is used.

Alternatively, a list of point of time (`tpoints`) can be given. Points of time
needs to be non-negative (you cannot go back in time). In this case a list of
resulting states is returned.

*Note:* It is up to the user to provide `evogen` and `state` fulfilling
the appropriate conditions. For the procedure to work correctly `evogen` should
be generated by `evolve_generator` function and `state` should be a proper
density matrix.

The third approach can be used if the superoperator is known. In this case
argument `evosuper` can be specified. This argument can be generated by
`evolve_operator` function. This is useful to simulate a fixed model of
evolution in the case of multiple initial states and the same time point.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
([0 1; 1 0], Array{Int64,2}[[0 1; 0 0], [0 0; 1 0]])

julia> Matrix(evolve(evolve_generator(H, L), proj(1, 2), 4.))
2×2 Array{Complex{Float64},2}:
 0.499815+0.0im              0.0+0.00127256im
      0.0-0.00127256im  0.500185+0.0im

julia> Matrix.(evolve(evolve_generator(H, L), proj(1, 2), [1., 2., 3., 4.]))
4-element Array{Array{Complex{Float64},2},1}:
 [0.43320327312718954 + 0.0im 0.0 + 0.10760477583156115im; 0.0 - 0.10760477583156117im 0.5667967268728107 + 0.0im]
 [0.48576602987832085 + 0.0im 0.0 - 0.017171799503962696im; 0.0 + 0.017171799503962682im 0.5142339701216796 + 0.0im]
 [0.5055971005015641 + 0.0im 0.0 - 0.002617013961759585im; 0.0 + 0.002617013961759615im 0.49440289949843613 + 0.0im]
 [0.49981547041444147 + 0.0im 0.0 + 0.0012725622225037549im; 0.0 - 0.0012725622225037724im 0.5001845295855593 + 0.0im]

julia> ev_op = evolve_operator(evolve_generator(H, L), 4.)
4×4 Array{Complex{Float64},2}:
 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im

julia> Matrix(evolve(ev_op, proj(1, 2)))
2×2 Array{Complex{Float64},2}:
 0.499815+0.0im              0.0+0.00127256im
      0.0-0.00127256im  0.500185+0.0im
```
"""
function evolve(exp_evolve_generator::AbstractMatrix{<:Number},
                initial_state::AbstractMatrix{<:Number})
  @argumentcheck size(exp_evolve_generator, 1) == size(exp_evolve_generator, 2) "Argument exp_evolve_generator should be square"
  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Initial_state should be a square matrix"
  @assert size(exp_evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of exp_evolve_generator size"

  unres(exp_evolve_generator*res(initial_state))
end,

function evolve(evolve_generator::AbstractMatrix{<:Number},
                initial_state::AbstractMatrix{<:Number},
                timepoint::Real)
  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "Argument evolve_generator should be square"
  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Initial_state should be a square matrix"
  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
  @argumentcheck timepoint>= 0 "Time needs to be nonnegative"

  unres(exp(timepoint*evolve_generator)*res(initial_state))
end,

function evolve(evolve_generator::SparseMatrixCSC{<:Number},
                initial_state::AbstractMatrix{<:Number},
                timepoint::Real)
  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "Argument evolve_generator should be a square matrix"
  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Argument initial_state should be a square matrix"
  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
  @argumentcheck timepoint>= 0 "Time needs to be nonnegative"

  unres(expmv(timepoint, evolve_generator, Vector(res(initial_state))))
end,

function evolve(evolve_generator::AbstractMatrix{<:Number},
                initial_state::AbstractMatrix{<:Number},
                timepoints::AbstractVector{<:Real})
  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "evolve_generator should be a square matrix"
  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Argument initial_state should be a square matrix"
  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
  @argumentcheck all(timepoints.>= 0) "All time points need to be nonnegative"

  [evolve(evolve_generator, initial_state, t) for t = timepoints]
end


text:
module TestShapes

import Modia3D
using  Modia3D.Test
using  Modia3D.StaticArrays

material1 = Modia3D.VisualMaterial(;color="DarkRed", wireframe=true)
material2 = Modia3D.VisualMaterial(;color=[0,0,255], shininess=0.2)

sphere = Modia3D.Visual(shape=Modia3D.Sphere(diameter=2.0),
    visualMaterial=material1)
ellipsoid = Modia3D.Visual(shape=Modia3D.Ellipsoid(lengthX=1.0, lengthY=2.0,
    lengthZ=3.0))
box = Modia3D.Visual(shape=Modia3D.Box(lengthX=1.0, lengthY=2.0,
    lengthZ=3.0))
cylinder = Modia3D.Visual(shape=Modia3D.Cylinder(axis=3, diameter=3.0,
    length=1.0, innerDiameter=1.5))
cone = Modia3D.Visual(shape=Modia3D.Cone(diameter=3.0, length=1.0,
    topDiameter=0.6), visualMaterial=material1)
capsule = Modia3D.Visual(shape=Modia3D.Capsule(axis=3, diameter=3.0,
    length=1.0))
beam = Modia3D.Visual(shape=Modia3D.Beam(axis=1, length=1.0, width=2.0,
    thickness=3.0), visualMaterial=material1)
spring = Modia3D.Visual(shape=Modia3D.Spring(axis=3, length=1.0, diameter=3.0,
    windings=10, wireDiameter=0.2), visualMaterial=material2)
gearWheel = Modia3D.Visual(shape=Modia3D.GearWheel(axis=3, diameter=3.0,
    length=1.0, innerDiameter=2.4, teeth=30, angle=0.2), visualMaterial=material2)
coordsys = Modia3D.Visual(shape=Modia3D.CoordinateSystem(length=2.0))
grid= Modia3D.Visual(shape=Modia3D.Grid(axis=3, length=3.0, width=1.0, distance=0.1, lineWidth=2.0))

filename = joinpath(Modia3D.path, "objects", "engine", "crank", "crank.obj")
mesh1 = Modia3D.Visual(shape=Modia3D.FileMesh(filename=filename))
mesh2 = Modia3D.Visual(shape=Modia3D.FileMesh(filename=filename, scale=[2.0,2.0,2.0],
    useMaterialColor=true, smoothNormals=true))
mesh3 = Modia3D.Visual(shape=Modia3D.FileMesh(filename=filename, scale=[2.0,3.0,4.0],
    useMaterialColor=true, smoothNormals=true))

text1 = Modia3D.TextShape(text="This is a text line")
font1 = Modia3D.Font()
font2 = Modia3D.Font(fontFamily="FreeSans", bold=false, italic=false, charSize=0.1,
    color = "LightBlue", transparency=0.0)
text2 = Modia3D.TextShape(text="Another text line", font=font2, offset=[0.1,0.2,0.3],
    axisAlignment=Modia3D.XZ_Plane,alignment=Modia3D.Left)

@testset "Modia3D.Shapes: Test shape, color, fonts" begin
   @test material1.color     == MVector{3,Cint}(139,0,0)
   @test material1.wireframe == true
end

end


text: DataFrames.Index(sch::Data.Schema) = DataFrames.Index(Symbol.(Data.header(sch)))


text:
############################################################################################
"""
$(TYPEDEF)

MCMC Tuning container.

# Fields
$(TYPEDFIELDS)
"""
struct MCMCTune{
    A<:BaytesCore.UpdateBool,
    B<:BaytesCore.UpdateBool,
    F<:AbstractFloat,
    T<:Tagged,
    E<:Tuple,
    P<:Proposal} <: AbstractTune
    "Tagged Parameter."
    tagged::T
    "Current Phase in MCMC Cycle"
    phase::PhaseTune{E}
    "Stepsize container"
    stepsize::StepSizeTune{A,F}
    "Information for posterior covariance estimate"
    proposal::P
    "Boolean if generated quantities should be generated while sampling"
    generated::B
    "Current iteration number"
    iter::Iterator
    function MCMCTune(
        objective::Objective,
        phase::PhaseTune{E},
        stepsize::StepSizeTune{A,F},
        proposal::P,
        generate::B,
    ) where {A<:BaytesCore.UpdateBool,B<:BaytesCore.UpdateBool,F<:AbstractFloat,E<:Tuple,P<:Proposal}
        #!NOTE: Start with 0, so first proposal step will update iter to 1
        iter = Iterator(0)
        return new{A,B,F,typeof(objective.tagged),E,P}(
            objective.tagged, phase, stepsize, proposal, generate, iter
        )
    end
end

############################################################################################
"""
$(SIGNATURES)
Update MCMC tuning fields at current iteration.

# Examples
```julia
```

"""
function update!(
    tune::MCMCTune, result::S, acceptrate::T
) where {S<:ℓObjectiveResult,T<:Real}
    @unpack θᵤ = result
    ##  Update Current iteration counter
    update!(tune.iter)
    ## Update current sampling phase
    #!NOTE: Done AFTER tune.iter update; both iterations start from 0, so to check if phase needs to be updated, tune.iter has to be updated before
    update!(tune.phase, tune.iter.current)
    ## Update stepsize and proposal parameter given current sampling phase
    update!(tune.stepsize, acceptrate, tune.phase)
    update!(tune.proposal, θᵤ, tune.phase)
    ## Pack container
    return nothing
end

############################################################################################
#=
"""
$(SIGNATURES)
Some basic tuning output figures. Useful for debugging.

# Examples
```julia
```

"""
function printtune(tune::MCMCTune, diagparam::Integer=size(tune.proposal.chain, 2))
    ArgCheck.@argcheck tune.iter.current > 0 "Need to run a proposal step before accessing tuning configuration."
    println("#################################################")
    println("Tune iter: ", tune.iter.current)
    println("Phase iter: ", tune.phase.iter.current)
    println("Phase name: ", tune.phase.name[tune.phase.counter.current])
    println("Phase update: ", tune.phase.update.current)
    println("Proposal update: ", tune.proposal.adaption)
    println("Phase counter: ", tune.phase.counter.current)
    #!NOTE: This will be Dualaveraging Parameter for first update! in new samplingphase, hence in final sampling phase only the smoothed ϵ is returned, and the dualaverage has no impact cause it wont compute anymore when Exploration+update=false
    #!NOTE: Hence, the stepsize multiplier will only come into effect at NEXT iteration with un-smoothed epsilon, in case adaption is still performed
    println("Phase stepsize: ", tune.stepsize.ϵ)
    println(
        "Phase (avg) stepsize: ",
        exp(tune.stepsize.dualaverage.logϵ),
        " (",
        exp(tune.stepsize.dualaverage.logϵ̄),
        ")",
    )
    println(
        "Proposal chain current row: ", tune.proposal.chain[end, tune.phase.iter.current]
    )
    println("Proposal chain last element: ", tune.proposal.chain[end])
    println("Proposal covariance: ", tune.proposal.Σ.diag[1:diagparam])
    return nothing
end
=#
############################################################################################
# Export
export MCMCTune


text:
import MKL
using SafeTestsets: @safetestset

@time @safetestset "layers" begin include("layers.jl") end
@time @safetestset "binary" begin include("binary.jl") end
@time @safetestset "pcd" begin include("pcd.jl") end
@time @safetestset "centered_gradient" begin include("centered_gradient.jl") end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(problem::EstimationProblem)
  pdata   = data(problem)
  pdomain = domain(problem)
  vars    = join(name.(variables(pdata)), ", ")

  title --> "Estimation Problem"
  legend --> true

  @series begin
    label --> "domain"
    pdomain
  end
  @series begin
    seriescolor --> :blue
    label --> "data ($vars)"
    marker --> :xcross
    domain(pdata)
  end
end


text:
"""
    ability()

Return a valid ability value.
"""
function ability()
    # While the exercise doesn't explicitely require it, we simulate 4 dice rolls and discard the smallest roll.
    # A minimal solution would be ability() = rand(3:18)

    # Throw 4 dice
    dice = [rand(1:6) for i in 1:4]

    # Sort the dice so that [2:end] excludes the lowest roll
    sort!(dice)
    sum(dice[2:end])
end

"""
    modifier(ability)

Return the modifier based on an ability value.
"""
modifier(ability) = fld(ability - 10, 2)

Base.@kwdef struct DNDCharacter
    strength::Int = ability()
    dexterity::Int = ability()
    constitution::Int = ability()
    intelligence::Int = ability()
    wisdom::Int = ability()
    charisma::Int = ability()
    hitpoints::Int = 10 + modifier(constitution)
end


text:
struct SearchTrace{S, V}
    s_init::S
    v_init::V
    df::DataFrame
end

SearchTrace(v) = SearchTrace(nothing, v, basedf(eltype(v)))
SearchTrace(s, v) = SearchTrace(s, v, basedf(eltype(v)))

function init!(st::SearchTrace, s, v)
    st.s_init = copy(s)
    st.v_init = v
    empty!(st.df)
end

function init!(st::SearchTrace{Nothing, V}, s, v) where {V}
    st.v_init = v
    empty!(st.df)
end

function Base.push!(st::SearchTrace, data::Tuple)
    push!(st.df, data)
    return st
end

basedf(V) = DataFrame(basecols(V))

basecols(V) = [:it => Int[], :op => Function[], :params => Tuple[], :v_i => V[], :v_s => V[], :accept => Bool[]]

function accepted(st::SearchTrace; accept::Symbol = :accept)
    filterdf(st, accept)
end

function accepted!(st::SearchTrace; accept::Symbol = :accept)
    filterdf!(st, accept)
end

function improving(st::SearchTrace; output::Symbol = :output)
    filterdf(st, output)
end

function improving!(st::SearchTrace; output::Symbol = :output)
    filterdf!(st, output)
end

function filterdf(st::SearchTrace, symbol::Symbol)
    SearchTrace(st.s_init, st.v_init, st.df[st.df[:, symbol], :])
end

function filterdf!(st::SearchTrace, symbol::Symbol)
    st.df = st.df[st.df[:, symbol], :]
end

function filterdf(st::SearchTrace, filter)
    SearchTrace(st.s_init, st.v_init, st.df[filter, :])
end

function filterdf!(st::SearchTrace, filter)
    st.df = st.df[filter, :]
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz

max_forks_repo_name: Michiel-VL/PolyaViz


text:
using Serialization, DataFrames, Unitful, FileIO

isotopes = deserialize("isotopes_data.jls")

quantity_repr(x) = repr(x)
function quantity_repr(x::Quantity)
	n,u = repr(x)|>split
	return n*"*"*u
end

open("isotopes_data.jl", "w") do io
	println(io, "_isotopes_data = [")
	for istp in eachrow(isotopes)
		println(io, "    Isotope(")
		for field in istp[1:end-1]
			println(io, "        ",quantity_repr(field),',')
		end
		println(io, "        ",quantity_repr(istp[end]))
		println(io, "    ),")
	end
	println(io, "]")
end

max_stars_repo_name: Gregstrq/Isotope-data

max_issues_repo_name: Gregstrq/Isotope-data

max_forks_repo_name: Gregstrq/Isotope-data


text:
# Single-source shortest path in complete graphs.
# Only the upper part of the reward matrix is considered!

abstract type ElementaryPathSolver end
function build!(solver::ElementaryPathSolver, graph::SimpleDiGraph, source::Int, destination::Int)
  nothing
end

struct ElementaryPath <: CombinatorialInstance{Tuple{Int, Int}}
  # Mandatory properties.
  n_arms::Int
  optimal_average_reward::Float64

  # Probability distributions for the arm rewards.
  graph::SimpleDiGraph
  reward::Dict{Tuple{Int, Int}, Distribution} # Roughly a matrix, except that not all entries should be filled.
  # Even for a complete graph, the matrix has too many elements: the diagonal is not used.
  # For a very large number of reads into the data structure, the dictionary seems a good fit.

  # Characteristics of the path.
  source::Int
  destination::Int

  # Internal solver.
  solver::ElementaryPathSolver

  function ElementaryPath(graph::SimpleDiGraph, reward::Dict{Tuple{Int, Int}, Distribution}, source::Int, destination::Int, solver::ElementaryPathSolver)
    n = nv(graph)

    if source <= 0
      error("Source node is not an acceptable node index: $(source) is negative or zero. The first node is numbered 1.")
    end
    if source > n
      error("Source node is not an acceptable node index: $(source) is too large. The last node is numbered $n.")
    end
    if destination <= 0
      error("Destination node is not an acceptable node index: $(destination) negative or zero. The first node is numbered 1.")
    end
    if destination > n
      error("Destination node is not an acceptable node index: $(destination) is too large. The last node is numbered $n.")
    end

    if length(reward) != ne(graph)
      error("The edges and rewards do not perfectly match: there are $(length(reward)) rewards and $(ne(graph)) edges.")
    end
    for v in keys(reward)
      if ! has_edge(graph, v[1], v[2])
        error("The edges and rewards do not perfectly match: the edge $v has a reward, but does not exist in the graph.")
      end
    end

    if is_cyclic(graph)
      @warn "The input graph is cyclic, bandit algorithms are usually not useful in this case"
      # Just play a positive-reward cycle an infinite number of times at each round...
    end

    # Prepare the solver to be used (if required).
    build!(solver, graph, source, destination)

    avg_reward = Dict{Tuple{Int, Int}, Float64}(k => mean(v) for (k, v) in reward)
    opt_sol = solve_linear(solver, avg_reward)
    opt = sum(avg_reward[i] for i in opt_sol)

    # Done!
    return new(ne(graph), opt, graph, reward, source, destination, solver)
  end
end

function initial_state(instance::ElementaryPath) # TODO: share a lot of code among graph instances?
  zero_counts = Dict(k => 0 for (k, _) in instance.reward)
  zero_rewards = Dict(k => 0.0 for (k, _) in instance.reward)
  return State{Tuple{Int, Int}}(0, 0.0, 0.0, zero_counts, zero_rewards, copy(zero_rewards))
end

function is_feasible(instance::ElementaryPath, arms::Vector{Tuple{Int, Int}})
  # Many checks are similar for just a partially acceptable path.
  if ! is_partially_acceptable(instance, arms)
    return false
  end

  if length(arms) == 0
    return false
  end

  # Check whether the path goes to the destination. Think about the case where
  # the path is just one edge. The check about the source is already performed
  # by is_partially_acceptable.
  last_edge = arms[end]
  if last_edge[2] != instance.destination
    return false
  end

  return true
end

function is_partially_acceptable(instance::ElementaryPath, arms::Vector{Tuple{Int, Int}})
  path_length = length(arms)

  if path_length > instance.n_arms
    return false
  end

  if path_length < 1 # Major difference with is_feasible: an empty path is partially acceptable!
    return true
  end

  n = nv(instance.graph)

  # Check whether the path starts from the source and possibly ends at the destination.
  # Think about the case where the path is just one edge.
  for (i, arm) in enumerate(arms)
    # The first edge must start at the source.
    if i == 1 && arm[1] != instance.source
      return false
    end

    # No edge can end at the destination if it is not the last one in the path (otherwise, cycle).
    if i < length(arms) && arm[2] == instance.destination
      return false
    end

    # The last edge must end at the destination, but only if the path is as long as possible.
    if i == n && path_length == instance.n_arms && arm[2] != instance.destination
      return false
    end

    # If after the first edge, this edge must be adjacent to the same node as the previous one.
    if i > 1 && arms[i - 1][2] != arm[1]
      return false
    end
  end

  return true
end


text:
function newton_solve(F, ∇ₓF, x; Ftol=1e-10)
    while norm(F(x)) ≥ Ftol
        x .-= ∇ₓF(x) \ F(x)
    end
    return x
end

# Create a type for the solver's algorithm
struct MyAlg <: DiffEqBase.AbstractSteadyStateAlgorithm end

# Overload DiffEqBase's solve function
function DiffEqBase.solve(prob::DiffEqBase.AbstractSteadyStateProblem,
                          alg::MyAlg;
                          Ftol=1e-10)
    # Define the functions according to DiffEqBase.SteadyStateProblem type
    p = prob.p
    t = 0
    x0 = copy(prob.u0)
    dx, df = copy(x0), copy(x0)
    F(x) = prob.f(dx, x, p, t)
    ∇ₓF(x) = prob.f(df, dx, x, p, t)
    # Compute `u_steady` and `resid` as per DiffEqBase using my algorithm
    x_steady = newton_solve(F, ∇ₓF, x0, Ftol=Ftol)
    resid = F(x_steady)
    # Return the common DiffEqBase solution type
    DiffEqBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
end

# Overload DiffEqBase's SteadyStateProblem constructor
function DiffEqBase.SteadyStateProblem(F, ∇ₓF, x, p)
    f(dx, x, p, t) = F(x, p)
    f(df, dx, x, p, t) = ∇ₓF(x, p)
    return DiffEqBase.SteadyStateProblem(f, x, p)
end


text:
@testset "1561.maximum-number-of-coins-you-can-get.jl" begin
    @test max_coins([29, 24, 17, 44, 44, 39, 12, 31, 20, 36, 34, 59, 29, 59, 50, 38, 16, 67, 26, 76]) == 305
    @test max_coins([9,8,7,6,5,1,2,3,4]) == 18
    @test max_coins([2,4,1,2,7,8]) == 9
    @test max_coins([3, 4, 6]) == 4
end


text:
function steadystate!(m::HetDSGEGovDebt;
                      βlo::S = 0.5*exp(m[:γ].scaledvalue)/(1 + m[:r].scaledvalue),
                      βhi::S = exp(m[:γ].scaledvalue)/(1 + m[:r].scaledvalue),
                      excess::S = 5000.,
                      tol::S = 1e-4,
                      maxit::Int64 = 20,
                      βband::S = 1e-2) where {S<:AbstractFloat}
    # If we have already solved for βstar (i.e. it's not NaN) and we only want to
    # estimate the non steady state parameters, there's no need to recompute
    # zlo/zhi, etc.
    if !isnan(m[:βstar].value) && get_setting(m, :estimate_only_non_steady_state_parameters)
        return
    else
        reset_grids!(m)
        nx = get_setting(m, :nx)
        ns = get_setting(m, :ns)
        ni = get_setting(m, :ni)
        nz = get_setting(m, :nz)

        # Determines whether random or not
        us, zs = if get_setting(m, :fix_random_matrices)
            get_setting(m, :us), get_setting(m, :zs)
        else
            generate_us_and_zs(ni, nz)
        end

        # This is a test setting; TODO: Remove because now we can just fix the randomness
        if get_setting(m, :steady_state_only)
            find_steadystate!(m; βlo = βlo, βhi = βhi, excess = excess, tol = tol, maxit = maxit,
                              βband = βband)
            return
        else
            # Do you want to calibrate for matching the income moments?
            if get_setting(m, :calibrate_income_targets)

                target = get_setting(m, :calibration_targets)
                lower  = get_setting(m, :calibration_targets_lb)
                upper  = get_setting(m, :calibration_targets_ub)
                m[:sH_over_sL], m[:zlo], _, _ = best_fit(m[:pLH].value, m[:pHL].value,
                                                         target, lower, upper, us, zs)
            else
                m[:varlinc], m[:vardlinc] = skill_moments(m[:sH_over_sL].value, m[:zlo].value,
                                                          m[:pLH].value, m[:pHL].value,
                                                          us, zs, ni)
            end
            # Whether one calibrates or not, zhi needs to be updated here
            m[:zhi] = 2.0 - m[:zlo].value

            # Construct Markov transition matrix for skill
            f, sgrid, swts, sscale = persistent_skill_process(m[:sH_over_sL].value, m[:pLH].value,
                                                              m[:pHL].value, get_setting(m, :ns))
            m.grids[:sgrid] = Grid(sgrid, swts, sscale)
            m.grids[:fgrid] = f

            xgrid, xwts, xlo, xhi, xscale = cash_grid(sgrid, m[:ωstar].value, m[:H].value,
                                                      m[:r].scaledvalue, m[:η].value,
                                                      m[:γ].scaledvalue,
                                                      m[:Tstar].value, m[:zlo].value, nx)

            m.grids[:xgrid] = Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)
            m <= Setting(:xlo, xlo)
            m <= Setting(:xhi, xhi)
            m <= Setting(:xscale, xscale)
            xswts = kron(swts, xwts)
            m.grids[:weights_total] = xswts

            # Once have updated grids, can call steady state and compute other two moments
            find_steadystate!(m; βlo = βlo, βhi = βhi,
                              excess = excess, tol = tol, maxit = maxit,
                              βband = βband)
            m[:mpc] = ave_mpc(m[:μstar].value,   m[:cstar].value, xgrid, xswts, nx, ns)
            m[:pc0] = frac_zero(m[:μstar].value, m[:cstar].value, xgrid, xswts, ns)
        end
    end
end

function find_steadystate!(m::HetDSGEGovDebt;
                           βlo::S = 0.5*exp(m[:γ].scaledvalue)/(1 + m[:r].scaledvalue),
                           βhi::S = exp(m[:γ].scaledvalue)/(1 + m[:r].scaledvalue),
                           excess::S = 5000.,
                           tol::S = 1e-4,
                           maxit::Int64 = 20,
                           βband::S = 1e-2) where {S<:AbstractFloat}
    # Load settings
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)

    xlo = get_setting(m, :xlo)
    xhi = get_setting(m, :xhi)

    zlo = m[:zlo].value
    zhi = m[:zhi].value

    # Load parameters
    R  = 1 + m[:r].scaledvalue
    H  = m[:H].value
    η  = m[:η].value
    γ  = m[:γ].scaledvalue
    ω  = m[:ωstar].value
    T  = m[:Tstar].value
    bg = m[:bg].value

    # Load grids
    xgrid = m.grids[:xgrid].points
    sgrid = m.grids[:sgrid].points
    xswts = m.grids[:weights_total]
    f     = m.grids[:fgrid]

    βlo_temp = βlo
    βhi_temp = βhi

    reject  = false
    counter = 1

    n  = ns*nx
    c  = zeros(n)
    bp = zeros(n)
    KF = zeros(n,n)
    μ  = zeros(n)

    # Initial guess
    β   = 1.0
    Win = Vector{Float64}(undef, n)
    Win_guess = ones(n)

    if get_setting(m, :use_last_βstar) && !isnan(m[:βstar].value)

        # This short-circuits computation of the policy function
        βlo = βhi = m[:βstar].value
    elseif !isnan(m[:βstar].value)

        # If one has computed β* before, we first bisect into a neighborhood around it
        βlo_temp = m[:βstar].value - βband
        βhi_temp = m[:βstar].value + βband

        c, bp, Win, KF, reject = policy_hetdsgegovdebt(nx, ns, βlo_temp, R, ω, H, η, T, γ,
                                                       zhi, zlo, xgrid, sgrid, xswts, Win_guess,
                                                       f, damp = get_setting(m, :policy_damp), maxit = get_setting(m, :policy_maxit))
        excess_lo, μ = compute_excess(xswts, KF, bp, bg)

        if excess_lo < 0 && abs(excess_lo) > tol
            βlo = βlo_temp
            c, bp, Win, KF, reject = policy_hetdsgegovdebt(nx, ns, βhi_temp, R, ω, H, η, T, γ, zhi,
                                                           zlo, xgrid, sgrid, xswts, Win_guess, f,
                                                           damp = get_setting(m, :policy_damp), maxit = get_setting(m, :policy_maxit))
            excess_hi, μ = compute_excess(xswts, KF, bp, bg)

            if excess_hi > 0
                βhi = βhi_temp
            end
        elseif excess_lo >= 0
            βhi = βlo_temp
        end
    end

    while abs(excess) > tol && counter < maxit # clearing markets
        β = (βlo + βhi) / 2.0
        c, bp, Win, KF, reject = policy_hetdsgegovdebt(nx, ns, β, R, ω, H, η, T, γ, zhi, zlo,
                                                       xgrid, sgrid, xswts, Win_guess, f,
                                                       damp = get_setting(m, :policy_damp), maxit = get_setting(m, :policy_maxit))
        excess, μ = compute_excess(xswts, KF, bp, bg)
        # bisection
        if excess > 0
            βhi = β
        elseif excess < 0
            βlo = β
        end
        counter += 1
        if get_setting(m, :use_last_βstar) && !isnan(m[:βstar].value)
            break
        end
    end

  #=  @show β
    This all was for testing when we were trying to compare the old steadystate to the new steadystate
    GR.inline("png")
    p = Plots.plot(xgrid, μ[1:300], label = "lowskill")
    Plots.plot!(p, xgrid, μ[301:600], label = "highskill")
    Plots.savefig(p, "agrid_vs_D_together_old.png")
    p = Plots.plot(xgrid, μ[1:300], label = "lowskill")
    Plots.savefig(p, "agrid_vs_D_lowskill_old.png")
    p = Plots.plot(xgrid, μ[301:600], label = "highskill")
    Plots.savefig(p, "agrid_vs_D_highskill_old.png")

    Plots.plot(μ)
    Plots.savefig("mu_old.png")

    p = Plots.plot(xgrid, c[1:300], label = "lowskill")
    Plots.plot!(p, xgrid, c[301:600], label = "highskill")
    Plots.savefig(p, "agrid_vs_cpol_old.png")

    p = Plots.plot(xgrid, xgrid - c[1:300], label = "lowskill")
    Plots.plot!(p, xgrid, xgrid - c[301:600], label = "highskill")
    Plots.savefig(p, "agrid_minus_cpol_old.png")

    p = Plots.plot()
    for mult = 1:6
        for i = 1:mult*100
            Plots.plot!(p, KF[:, i])
        end
        Plots.savefig(p,"KF_old_mult=$mult.png")
    end
    for i = 1:600
        p = Plots.plot(KF[:, i])
        Plots.savefig(p, "KF_dir/KF_old_$i.png")
    end=#

    # If policy function does not converge, we signal to likelihood that should reject
    m <= Setting(:auto_reject, reject)
    m[:lstar]  = Win
    m[:cstar]  = c
    m[:μstar]  = μ
    m[:βstar]  = β
    m[:β_save] = β
    nothing
end

# These are for computing moments to compare against targets (the next few functions)
function ave_mpc(m::AbstractArray, c::AbstractArray, agrid::AbstractArray,
                 aswts::AbstractArray, na::Int, ns::Int)
	mpc = 0.
	for iss=1:ns
		i = na*(iss-1)+1
		mpc += aswts[i]*m[i]*(c[i+1] - c[i])/(agrid[2] - agrid[1])
		for ia=2:na
			i = na*(iss-1)+ia
			mpc += aswts[i]*m[i]*(c[i] - c[i-1])/(agrid[ia] - agrid[ia-1])
		end
	end
	return mpc
end

function frac_zero(m::AbstractArray, c::AbstractArray, agrid::AbstractArray,
                   aswts::AbstractArray, ns::Int)
	return sum(aswts .* m .* (c .== repeat(agrid, ns)))
end

function ssample(us::Matrix{S}, P::Matrix{S}, πss::AbstractArray,
                 ni::Int) where {S<:AbstractFloat}
	shist = ones(Int,ni,8)
	for i=1:ni
		if us[i,1] > πss[1]
			shist[i,1] = 2
		end
		for t=2:8
			if us[i,t] > P[shist[i,t-1],1]
				shist[i,t] = 2
			end
		end
	end
	return shist
end

function ln_annual_inc(zhist::Matrix{S}, us::Matrix{S}, zlo::S, P::Matrix{S},
                       πss::AbstractArray, sgrid::AbstractArray,
                       ni::Int) where {S<:AbstractFloat}
  	s_inds = ssample(us,P,πss,ni)
	linc1 = zeros(ni)
	linc2 = zeros(ni)
	for i=1:ni
		inc1 = 0.
		for t=1:4
			zshock = 1. + (1. - zlo)*(zhist[i,t]-1.)
			sshock = (sgrid[1] + (sgrid[2] - sgrid[1])*(s_inds[i,t] - 1.))
			inc1 += zshock*sshock
		end
		inc2 = 0.
		for t=5:8
			zshock = 1. + (1. - zlo)*(zhist[i,t]-1.)
			sshock = (sgrid[1] + (sgrid[2] - sgrid[1])*(s_inds[i,t] - 1.))
			inc2 += zshock*sshock
		end
		linc1[i] += log(inc1)
     	linc2[i] += log(inc2)
	end
	return linc1, linc2
end

function skill_moments(sH_over_sL::Real, zlo::Real, pLH::S, pHL::S, us::Matrix{S},
                       zs::Matrix{S}, ni::Int = 10000) where {S<:AbstractFloat}
	πL    = pHL / (pLH + pHL)
	πss   = [πL; 1.0-πL]
	P     = [[1.0-pLH pLH]; [pHL 1.0-pHL]]
	slo   = 1.0 / (πL + (1-πL) * sH_over_sL)
	shi   = sH_over_sL * slo
	sgrid = [slo; shi]
	linc1, linc2 = ln_annual_inc(zs, us, zlo, P, πss, sgrid, ni)
	return var(linc1), var(linc2 - linc1)
end

"""
```
function loss(x::Vector{S}, target::Vector{S}) where {S<:AbstractFloat}
```

Computes for given vector and target vector, the sum of absolute loss.
"""
loss(x::Vector{S}, target::Vector{S}) where {S<:AbstractFloat} = sum(abs.(x-target))

"""
```
function best_fit(pLH::S, pHL::S, target::Vector{S}, lower::Vector{S}, upper::Vector{S},
                  us::Matrix{S}, zs::Matrix{S}, max_iter::Int = 20,
                  initial_guess::Vector{S} = [6.3, 0.03]) where {S<:AbstractFloat}
```

Uses Nelder-Mead (gradient descent algorithm) to optimize for income moments.
"""
function best_fit(pLH::S, pHL::S, target::Vector{S}, lower::Vector{S}, upper::Vector{S},
                  us::Matrix{S}, zs::Matrix{S}, max_iter::Int = 20,
                  initial_guess::Vector{S} = [6.3, 0.03]) where {S<:AbstractFloat}

    skill_moments_f(x) = loss(collect(skill_moments(x[1], x[2], pLH, pHL, us, zs)), target)

    res = optimize(skill_moments_f, lower, upper, initial_guess, Fminbox(NelderMead()),
                   Optim.Options(f_calls_limit = max_iter))

    sH_over_sL_argmin, zlo_argmin = Optim.minimizer(res)
    min_varlinc, min_vardlinc = skill_moments(sH_over_sL_argmin, zlo_argmin, pLH,
                                              pHL, us, zs)
    return sH_over_sL_argmin, zlo_argmin, min_varlinc, min_vardlinc
end

function zsample(uz::Matrix{S}, zgrid::AbstractArray, zcdf::AbstractArray,
                 ni::Int, nz::Int) where {S<:AbstractFloat}
    zave = 0.5*zgrid[1:nz-1]+0.5*zgrid[2:nz]
	zs = zeros(ni,8)
	for i=1:ni
		for t=1:8
			for iz=1:nz-1
				if zcdf[iz] < uz[i,t] <= zcdf[iz+1]
					zs[i,t] = zave[iz]
				end
			end
		end
	end
	return zs
end

# This function computes excess savings. The KF that pops out of steadystate is nx*ns by nx*ns but we want to reduce that to a vector of length nx*ns which is the stationary distribution (this is what the eigenstuff is doing). Stationary because no matter where you are now, you want to have the same probability of going to all the different gridpoints tomorrow.
@inline function compute_excess(xswts::Vector{S}, KF::Matrix{S}, bp::Vector{S},
                                bg::S; print_warning::Bool = false,
                                tol::S = 2e-1) where {S<:Float64}
    LPMKF = xswts[1] * KF

    # Find eigenvalue closest to 1
    (D,V) = (eigen(LPMKF)...,)
    max_D = argmax(abs.(D))
    D = D[max_D]

    if abs(D - 1) > tol && print_warning
        @warn "Your eigenvalue is too far from 1, something is wrong."
    end

    # Pick eigenvector associated w/ largest eigenvalue and moving it back to values
    μ = real(V[:,max_D])
    # @show μ
    μ = μ ./ dot(xswts, μ) # Scale of eigenvectors not determinate: rescale to integrate to 1
    excess = dot(xswts, (μ .* bp)) - bg # Compute excess supply of savings, which is a fn of w
    return excess, μ
end

# policy function (consumption and ell = 1/consumption are the "policies" being decied on here)
function policy_hetdsgegovdebt(nx::Int, ns::Int, β::S, R::S, ω::S, H::S, η::S,
                               T::S, γ::S, zhi::S, zlo::S, xgrid::Vector{S},
                               sgrid::Vector{S}, xswts::Vector{S}, Win::Vector{S},
                               f::Matrix{S}, dist::S = 1., tol::S = 1e-4;
                               maxit::Int64 = 500, damp::S = 0.5) where {S<:AbstractFloat}
    n    = nx*ns
    c    = zeros(n)                  # consumption
    bp   = Vector{Float64}(undef, n) # savings
    Wout = Vector{Float64}(undef, length(Win))
    counter = 1
    reject = false
    qfunction(x::Float64) = mollifier_hetdsgegovdebt(x, zhi, zlo)

    while dist>tol && counter<maxit
        # compute c(w) given guess for Win = β*R*E[u'(c_{t+1})]
        for iss in 1:ns
            for ia in 1:nx
                c[nx*(iss-1)+ia] = min(1/Win[nx*(iss-1)+ia],xgrid[ia]+η)
            end
        end

        bp = R*(exp(-γ))*(repeat(xgrid, ns) - c)  # compute bp(w) given guess for Win
        # Wout and Win are the ell function (in the paper!). The point is we're trying to get the policy function (ell) to converge (which is the business with Win and Wout because you need to compare the difference between Win and Wout to see if it converged)
        Wout = parameterized_expectations_hetdsgegovdebt(nx, ns, β, R, ω, H, T, γ,
                                                         qfunction, xgrid,
                                                         sgrid, xswts, c, bp, f)
        dist = maximum(abs.(Wout - Win))
        Win  = damp*Wout + (1.0-damp) * Win
        counter += 1
    end
    if counter == maxit
        @warn "Euler iteration did not converge"
        reject = true
        return c, bp, Wout, zeros(n,n), reject
    end
    tr = kolmogorov_fwd_hetdsgegovdebt(nx, ns, ω, H, T, R, γ, qfunction, xgrid, sgrid, bp, f)
    return c, bp, Wout, tr, reject
end

# This function is the formula for ell_t in the middle top of page 7 of the paper
@inline function parameterized_expectations_hetdsgegovdebt(nx::Int, ns::Int, β::S,
                                                           R::S, ω::S, H::S, T::S, γ::S,
                                                           qfunc::Function, xgrid::Vector{S},
                                                           sgrid::Vector{S}, xswts::Vector{S},
                                                           c::Vector{S}, bp::Vector{S},
                                                           f::Matrix{S}) where {S<:AbstractFloat}
    l_out = zeros(nx*ns)
    for iss=1:ns     # loop off s (only 2 values)
        for ia=1:nx # loop over a
            sumn = 0.0
            for isp=1:ns # loop over s' (only 2 values)
                for iap=1:nx # loop over a'
                    # Name in code                                                 : notation in paper
                    #xswts                                                         : weights for the integral approximation
                    #c[nx*(isp-1)+iap]                                             : c_{t+1}(a', s')
                    # qfunc                                                        : g(⋅)
                    # qfunc((xgrid[iap] - bp[nx*(iss-1)+ia] - T) / (ω*H*sgrid[isp]): (a'-(1+r)exp(-γ)(a-c_t(a,s) - T))/ w*s'*H
                    # f[iss, isp]/ sgrid[isp]                                      : p(s'|s)/s'
                    # β*R*exp(-γ))/ω*H  (line 429)                                 β*(1+r)*exp(-z)/ω_{t+1}*H{t+1}
                    # When you're trying to match code to paper, be creative in the ordering of stuff as the variables in the paper are commutatively multipled in different orders
                    sumn += (xswts[nx*(isp-1)+iap]/c[nx*(isp-1)+iap]) *
                        qfunc((xgrid[iap] - bp[nx*(iss-1)+ia] - T) /
                              (ω*H*sgrid[isp])) * f[iss,isp] ./ sgrid[isp]
                end
            end
            l_out[nx*(iss-1)+ia] = (β*R*(exp(-γ))/ω*H)*sumn
        end
    end
    return l_out
end

# This function is the exprression for m_{t+1} in the middle of page 7 of the paper. See the key in lines 417-424 of this file to help you in matching stuff up
@inline function kolmogorov_fwd_hetdsgegovdebt(nx::Int, ns::Int, ω::S, H::S, T::S, R::S,
                                               γ::S, qfunc::Function, xgrid::Vector{S},
                                               sgrid::Vector{S}, bp::Vector{S},
                                               f::Matrix{S}) where {S<:AbstractFloat}
    tr = zeros(nx*ns,nx*ns)
    for iss=1:ns
        for ia=1:nx
            for isp=1:ns
                for iap=1:nx
                    tr[nx*(isp-1)+iap, nx*(iss-1)+ia] = qfunc((xgrid[iap] -
                                                               bp[nx*(iss-1)+ia] - T) /
                                                              (ω*H*sgrid[isp])) * f[iss,isp] ./
                                                              (ω * H * sgrid[isp])
                end
            end
        end
    end
    return tr
end

# also known as qfunction of g(⋅)
@inline function mollifier_hetdsgegovdebt(z::S, ehi::S, elo::S) where {S<:AbstractFloat}
    In = 0.443993816237631
    if z<ehi && z>elo
        temp = -1.0 + 2.0 * (z - elo) / (ehi - elo)
        return (2.0 / (ehi - elo)) * exp(-1.0 / (1.0 - temp^2)) / In
    end
    return 0.0
end
# also known as qfunction of g(⋅) but the derivative of it
@inline function dmollifier_hetdsgegovdebt(x::S, ehi::S, elo::S) where {S<:AbstractFloat}
    In = 0.443993816237631
    if x<ehi && x>elo
        temp = (-1.0 + 2.0*(x-elo)/(ehi-elo))
        out  = -(2*temp./((1 - temp.^2).^2)).*(2/(ehi-elo)) .*
            mollifier_hetdsgegovdebt(x, ehi, elo)
    else
        out = 0.0
    end
    return out
end

function calibrate_pLH_pHL(m::HetDSGEGovDebt)
    target_mpc, target_pc0 = get_setting(m, :targets)
    σt_mpc, σt_pc0         = get_setting(m, :target_σt)

    function minimize_penalty(m, pLHpHL::Vector{Float64})
        m[:pLH] = pLHpHL[1]
        m[:pHL] = pLHpHL[2]
        steadystate!(m)
        if get_setting(m, :auto_reject)
            m <= Setting(:auto_reject, false)
            @show pLHpHL, Inf
            return Inf
        end
        @show pLHpHL, m[:mpc].value, m[:pc0].value
        println(0.5*((log(m[:mpc])-log(target_mpc))^2/σt_mpc +
                    (log(m[:pc0].value)-log(target_pc0))^2/σt_pc0^2))
        return 0.5*((log(m[:mpc])-log(target_mpc))^2/σt_mpc +
                    (log(m[:pc0].value)-log(target_pc0))^2/σt_pc0^2)
    end
    minimize_me(x) = minimize_penalty(m, x)
    res = optimize(minimize_me, [0.005, 0.005], [0.095, 0.095], [0.01125, 0.03],
                   Fminbox(Optim.NelderMead()), Optim.Options(f_calls_limit=300))
    println(res)
    return res
end


text:
module NoMeta

"""
This docstring contains a doctest that needs a DocTestSetup

```jldoctest
julia> baz(20)
40
```
"""
function foo end

"""
This docstring contains a doctest that needs a DocTestSetup, but it's provided:

```@meta
DocTestSetup = quote
    qux(x) = 3x
end
```

```jldoctest
julia> qux(10)
30
```

```@meta
DocTestSetup = nothing
```
"""
function bar end

end


text:
module ForwardDiff2

include("dualnumber.jl")
include("tag.jl")
include("dualarray.jl")
include("custom_dispatch.jl")
include("dual_context.jl")
include("api.jl")

# Experimental
#include("experiment/aosoa.jl")

end # module


text:
"""
    RandomSampler{T<:AbstractRNG} <: Sampler

Sample a value for each parameter uniformly at random from the candidate vectors. Log-uniform sampling available by providing a log-spaced candidate vector.

Optionally, pass an `AbstractRNG` to initialize.
"""
struct RandomSampler{T<:AbstractRNG} <: Sampler
    rng_channel::RemoteChannel{Channel{T}}
    lock::ReentrantLock
end

RandomSampler() = RandomSampler(MersenneTwister(rand(1:1000)))

function RandomSampler(rng::T) where {T <: AbstractRNG}
    # We use a *single* RNG across threads and processors to ensure
    # that e.g. different processors don't have different copies of the
    # RNG.
    channel = RemoteChannel(()->Channel{T}(1), 1)
    put!(channel, rng)
    return RandomSampler(channel, ReentrantLock())
end

function (s::RandomSampler)(ho, iter)
    # Lock first, since RemoteChannel's may not be threadsafe:
    # <https://github.com/JuliaLang/julia/issues/37706>
    return lock(s.lock) do
        rng = take!(s.rng_channel)
        result = [list[rand(rng, 1:length(list))] for list in ho.candidates]
        put!(s.rng_channel, rng)
        return result
    end
end

# Latin Hypercube Sampler ======================================================

"""
Sample from a latin hypercube
"""
Base.@kwdef mutable struct LHSampler <: Sampler
    samples = zeros(0,0)
    iters = -1
end
function init!(s::LHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    ndims = length(ho.candidates)
    all(length(c) == ho.iterations for c in ho.candidates) ||
        throw(ArgumentError("Latin hypercube sampling requires all candidate vectors to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, [length(c) for c in ho.candidates]))))) with $(ho.iterations) iterations"))
        s.iters == -1 && (s.iters = (1000*100*2)÷ho.iterations÷ndims)
    X, fit = LHCoptim(ho.iterations,ndims,s.iters)
    s.samples = copy(X')
end

function (s::LHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

"""
    CLHSampler(dims=[Continuous(), Categorical(2), ...])
Sample from a categorical/continuous latin hypercube. All continuous variables must have the same length of the candidate vectors.
"""
Base.@kwdef mutable struct CLHSampler <: Sampler
    samples = zeros(0,0)
    dims = []
end
function init!(s::CLHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    all(zip(s.dims, ho.candidates)) do (d,c)
        d isa Categorical || length(c) == ho.iterations
    end || throw(ArgumentError("Latin hypercube sampling requires all candidate vectors for Continuous variables to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates)))))"))
    ndims = length(ho.candidates)
    initialSample = randomLHC(ho.iterations,s.dims)
    X,_ = LHCoptim!(initialSample, 500, dims=s.dims)
    s.samples = copy(X')
end

function (s::CLHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

# Hyperband ====================================================================

"""

"""
Base.@kwdef mutable struct Hyperband <: Sampler
    R
    η = 3
    minimum = (Inf,)
    inner::Sampler = RandomSampler()
end
Hyperband(R) = Hyperband(R=R)

function optimize(ho::Hyperoptimizer{<:Hyperband})
    try
        hyperband(ho)
    catch e
        if e isa InterruptException
            @info "Aborting hyperband"
        else
            rethrow(e)
        end
    end
    ho
end

function hyperband(ho::Hyperoptimizer{Hyperband}; threads=false)
    hb = ho.sampler
    R, η = hb.R, hb.η
    hb.minimum = (Inf,)
    smax = floor(Int, log(η,R))
    B = (smax + 1)*R # B is budget

    p = Progress(smax+1, 1, "Hyperband")
    for s in smax:-1:0
        n = ceil(Int, (B/R)*((η^s)/(s+1)))
        r = R / (η^s)
        minᵢ = successive_halving(ho, n, r, s; threads)
        if minᵢ[1] < hb.minimum[1]
            hb.minimum = minᵢ
        end
        ProgressMeter.next!(p; showvalues = [("bracket (of $(smax+1))",smax-s+1), ("minimum", hb.minimum[1]), ("with resources", hb.minimum[2]), ("minizer", hb.minimum[3])])
    end
    return hb.minimum
end

function successive_halving(ho, n, r=1, s=round(Int, log(hb.η, n)); threads=false)
    hb = ho.sampler
    costfun = ho.objective
    η = hb.η
    minimum = Inf
    T = [ hb.inner(ho, i) for i=1:n ]
    mapfun = threads ? ThreadPools.tmap : map
    for i in 0:s
        nᵢ = n/(η^i) # Flooring takes place below
        rᵢ = r*(η^i)
        if i == 0
            LTend = mapfun(t->costfun(rᵢ, t...), T)
        else
            LTend = mapfun(t->costfun(rᵢ, t), T)
        end
        L, T = first.(LTend), last.(LTend)
        append!(ho.history, T)
        append!(ho.results, L)
        if hb.inner isa BOHB
            update_observations(ho, rᵢ, T, L)
        end
        perm = sortperm(L)
        besti = perm[1]
        if L[besti] < minimum[1]
            minimum = (L[besti], rᵢ, T[besti])
        end
        T = T[perm[1:floor(Int,nᵢ/η)]]
        Base.CoreLogging.@logmsg -1 "successive_halving" progress=i/s
    end
    return minimum
end

function hyperband(f, candidates; R, η=3, inner = RandomSampler(), threads=false)
    sampler = Hyperband(; R, η, inner)
    objective(i, pars...) = f(i, [pars...]) # The objective needs two methods, accepting a vector and a list of args
    objective(i, state) = f(i, state)
    ho = Hyperoptimizer(;
        iterations = 1,
        params = [Symbol("$i") for i in eachindex(candidates)],
        candidates,
        history = [],
        results = Real[],
        sampler,
        objective,
    )
    if inner isa Union{LHSampler,CLHSampler}
        ho.iterations = length(candidates[1])
        init!(inner, ho)
    end
    hyperband(ho; threads)
    ho
end

# BOHB ====================================================================
# Acknowledgement: Code structure refers to official implementation of BOHB in ['HpBandSter'](https://github.com/automl/HpBandSter)
#
# Copyright of HpBandSter:
#
#     Copyright (c) 2017-2018, ML4AAD
#     All rights reserved.

# struct to record BOHB Observation
mutable struct ObservationsRecord
    dim::Int
    observation::Union{Vector, Tuple}
    loss::Real
end
function ObservationsRecord(observation, loss)
    ObservationsRecord(length(observation), observation, loss)
end

"""
    BOHB samplers
All variable names refer symbols in the [paper](`https://arxiv.org/pdf/1807.01774v1.pdf`)
- `ρ`: Fraction of random samples
- `q`: Fraction of best observations to build l and g
- `N_s`: Sample batch number
- `N_min`: Minimum number of points to build a model
- `bw_factor`: Bandwidth factor
- `D`: Evaluated observations
- `max_valid_budget`: Maximum budget i that |D_{i}| is big enough to fit a model
- `N_b`: |D_{max_valid_budget}|
- `KDE_good`: KDE consists of "good observations", see BOHB paper
- `KDE_bad`: KDE consists of "bad observations", see BOHB paper
"""
Base.@kwdef mutable struct BOHB <: Sampler
    dims::Union{Vector{DimensionType}, Nothing}=nothing
    # hyperparameters for BOHB
    N_min::Union{Int, Nothing} = nothing
    ρ::AbstractFloat = 1/3
    q::AbstractFloat = 0.15
    N_s::Int = 64
    bw_factor::Real = 3
    "minimum bandwidth: this parameter doesn't occur in the paper but used in the official implementation"
    min_bandwidth::Real = 1e-3
    "Random sampler used for random sampling in BOHB algorithm"
    random_sampler::RandomSampler = RandomSampler()
    # Context data
    ## Current observations, stored in a Dict, in which key is budget, value is a observation array to fit KDEs
    ## key of D: A real number represents budget
    ## value of D: An vector of ObservationsRecord, all the records of corresponding budget
    D::Dict{Real, Vector{ObservationsRecord}} = Dict{Real, Vector{ObservationsRecord}}()
    "Current maximum budget that |D_{b}| > N_{min}+2, means it is valid for fit KDEs"
    max_valid_budget::Union{Number, Nothing} = nothing
    "|D| of max_valid_budget"
    N_b::Union{Int, Nothing} = nothing
    ## Good and bad kernel density estimator
    KDE_good::Union{MultiKDE.KDEMulti, Nothing} = nothing
    KDE_bad::Union{MultiKDE.KDEMulti, Nothing} = nothing
end

# object call of BOHB sampler
function (s::BOHB)(ho, iter)
    # with probability ρ, return random sampled observations.
    # If max_valid_budget is nothing, which means currently we don't have enough sample for TPE, random sample as well.
    if rand() < s.ρ || s.max_valid_budget === nothing
        return s.random_sampler(ho, iter)
    end
    potential_samples = [sample_potential_hyperparam(s.KDE_good, s.min_bandwidth, s.bw_factor) for _ in 1:s.N_s]
    scores = [score(sample, s.KDE_good, s.KDE_bad) for sample in potential_samples]
    _, best_idx = findmax(scores)
    [potential_samples[best_idx]]
end

# Sample score l(x)/g(x), refers to line 6 of Algorithm2 in paper
function score(sample::Vector, KDE_good::MultiKDE.KDEMulti, KDE_bad::MultiKDE.KDEMulti)
    pdf(KDE_good, sample) / pdf(KDE_bad, sample)
end

# Update budget observations in BOHB
function update_observations(ho::Hyperoptimizer{Hyperband}, rᵢ, observations, losses)
    bohb = ho.sampler.inner
    if !haskey(bohb.D, rᵢ)
        bohb.D[rᵢ] = []
    end
    for (c, l) in zip(observations, losses)
        push!(bohb.D[rᵢ], ObservationsRecord(c, l))
    end
    D_length = length(bohb.D[rᵢ])
    if bohb.N_min === nothing
        bohb.N_min = length(ho.candidates)+1
    end
    if D_length > bohb.N_min+2 && (bohb.max_valid_budget===nothing || rᵢ >= bohb.max_valid_budget)
        bohb.max_valid_budget, bohb.N_b = rᵢ, D_length
        update_KDEs(ho)
    end
end

function update_KDEs(ho::Hyperoptimizer{Hyperband})
    bohb = ho.sampler.inner
    records = bohb.D[bohb.max_valid_budget]
    # fit KDEs according to Eqs. (2) and (3) in paper
    N_bl = max(bohb.N_min, floor(Int, bohb.q*bohb.N_b))
    N_bg = max(bohb.N_min, bohb.N_b-N_bl)
    sort_idx = sortperm(records, by=d->d.loss)
    idx_N_bl = sort_idx[begin:N_bl]
    idx_N_bg = reverse(sort_idx)[N_bg:end]
    bohb.KDE_good = KDEMulti(bohb.dims, records[idx_N_bl], bohb.min_bandwidth, ho.candidates)
    bohb.KDE_bad = KDEMulti(bohb.dims, records[idx_N_bg], bohb.min_bandwidth, ho.candidates)
end

# sample from KDEMulti
function sample_potential_hyperparam(kde::MultiKDE.KDEMulti, min_bandwidth, bw_factor)
    idx = rand(1:size(kde.mat_observations)[2])
    param = [kde.observations[i][idx] for i in 1:length(kde.observations)]
    sample = Vector()
    for (_i, _param, dim_type, _kde) in zip(1:length(kde.dims), param, kde.dims, kde.KDEs)
        bw = max(_kde.bandwidth, min_bandwidth)
        local ele
        if dim_type isa MultiKDE.ContinuousDim
            bw = bw*bw_factor
            ele = rand(truncated(Normal(_param, bw), -_param/bw, (1-_param)/bw))
        elseif dim_type isa Union{MultiKDE.CategoricalDim, MultiKDE.UnorderedCategoricalDim}
            ele = rand() < (1-bw) ? _param : rand(1:dim_type.levels)
        else
            error(string("Dim type ", string(dim_type), " not supported. "))
        end
        if kde.mapped[_i]
            ele = kde.index_to_unordered[_kde][ele]
        end
        push!(sample, ele)
    end
    sample
end

# Constructor extensions and adapters for MultiKDE.jl
const DIMENSION_TYPE = Dict(Categorical=>MultiKDE.CategoricalDim, Continuous=>MultiKDE.ContinuousDim, UnorderedCategorical=>MultiKDE.UnorderedCategoricalDim)

function MultiKDE.KDEMulti(dim_types::Vector{DimensionType}, records::Vector{ObservationsRecord}, min_bandwidth::Real, candidates::Tuple)
    # Get KDEMulti with min_bandwidth
    dim = records[1].dim
    observations = Vector{Vector}()
    for record in records
        @assert record.dim == dim "All observations need to be same dimension. "
        _observations = record.observation
        if _observations isa Tuple
            _observations = [_obs for _obs in _observations]
        end
        push!(observations, _observations)
    end
    multi_kde = KDEMulti(dim_types, observations, candidates)
    for i in 1:length(multi_kde.KDEs)
        multi_kde.KDEs[i].bandwidth = max(multi_kde.KDEs[i].bandwidth, min_bandwidth)
    end
    multi_kde
end

function MultiKDE.KDEMulti(dims::Vector{DimensionType}, observations::Vector, candidates::Tuple)
    dims = Vector{MultiKDE.DimensionType}([DIMENSION_TYPE[typeof(dim)] === MultiKDE.ContinuousDim ? DIMENSION_TYPE[typeof(dim)]() :
                                            DIMENSION_TYPE[typeof(dim)](dim.levels) for dim in dims])
    MultiKDE.KDEMulti(dims, observations, candidates)
end

function Base.getproperty(dim_type::Union{MultiKDE.CategoricalDim, MultiKDE.UnorderedCategoricalDim}, v::Symbol)
    if v === :levels
        return getfield(dim_type, :level)
    end
    getfield(dim_type, v)
end


text:
struct MessagePassing{A,B,C} <: Mill.AggregationFunction
    instance_model::A
    aggregation::B
    vertex_model::C
    n::Int
end
Flux.@treelike(MessagePassing)
(m::MessagePassing)(x, bags) = message_passing(x, bags, m.instance_model, m.aggregation, m.vertex_model, m.n)
(m::MessagePassing)(ds::BagNode) = ArrayNode(message_passing(ds.data.data, ds.bags, m.instance_model, m.aggregation, m.vertex_model, m.n))
(m::MessagePassing)(x::ArrayNode, bags) = ArrayNode(message_passing(x.data, bags, m.instance_model, m.aggregation, m.vertex_model, m.n))
function message_passing(xx, bags, instance_model, aggregation, vertex_model, n::Int)
    onepass(x, bags) = vertex_model(aggregation(instance_model(x), bags))
    x = onepass(xx, bags)
    for i in 1:n-1
        x = onepass(x, bags)
    end
    x
end
message_passing(ds::BagNode, instance_model, aggregation, vertex_model, n::Int) = message_passing(ds.data.data, ds.bags, instance_model, aggregation, vertex_model, n)

# #test case
# bags = ScatteredBags([setdiff(neighbors(g, i), i) for i in vertices(g)])
# ds = BagNode(vertex_properties, bags)
# k = 4
# rounds = 3
# instance_model = Dense(k, k, relu)
# aggregation = SegmentedMeanMax(k)
# vertex_model = Dense(2k, k)
# message_passing(ds, instance_model, aggregation, vertex_model, rounds) # put this to test
# a = Mill.Aggregation((MessagePassing(instance_model, aggregation, vertex_model, rounds),))
# m = BagModel(BagModel(identity, a, identity), SegmentedMeanMax(k), Chain(Dense(2k, k, relu), Dense(k,2)))
# dss = BagNode(ds, [1:nobs(ds)])
# m(dss)
# m(cat(dss,dss)).data


text:
module MaxwellSALT

using Reexport
@reexport using MaxwellFDFD
using SALTBase; export gen_gain_2lv
using AbbreviatedTypes

include("mxsolver/mxsolver.jl")
include("gainobj.jl")
include("model.jl")

end # module MaxwellSALT


text:
struct Composition{F<:Function}
    code::Dict{Symbol,String}
    f::F
    symbols::Vector{Vector{Symbol}}
end

function Composition(f::F, symbols) where {F<:Function}
    code = Dict{Symbol,String}()
    return Composition{F}(code, f, symbols)
end

function code(c::Composition, lang=:maths; name="composition")
    return get!(c.code, lang, generate(c, name, Val(lang)))
end
composition(c::Composition) = c.f
symbols(c::Composition) = c.symbols

"""
    compose(icn, weights=nothing)
Return a function composed by some of the operations of a given ICN. Can be applied to any vector of variables. If `weights` are given, will assign to `icn`.
"""
function compose(icn::ICN, weigths::BitVector=BitVector())
    !isempty(weigths) && weights!(icn, weigths)
    composition, symbols = _compose(icn)
    return Composition(composition, symbols)
end

function generate(c::Composition, name, ::Val{:maths})
    aux = map(s -> reduce_symbols(s, ", ", length(s) > 1), symbols(c))
    def = reduce_symbols(aux, " ∘ ", false)
    return "$name = $def"
end

function generate(c::Composition, name, ::Val{:Julia})
    symbs = symbols(c)
    @assert length(symbs) == 4 "Length of the decomposition ≠ 4"
    tr_length = length(symbs[1])

    CN = "CompositionalNetworks."
    tr = reduce_symbols(symbs[1], ", "; prefix=CN * "tr_")
    ar = reduce_symbols(symbs[2], ", ", false; prefix=CN * "ar_")
    ag = reduce_symbols(symbs[3], ", ", false; prefix=CN * "ag_")
    co = reduce_symbols(symbs[4], ", ", false; prefix=CN * "co_")

    documentation = """\"\"\"
        $name(x; X = zeros(length(x), $tr_length), param=nothing, dom_size)

    Composition `$name` generated by CompositionalNetworks.jl.
    ```
    $(code(c; name))
    ```
    \"\"\"
    """

    output = """
    function $name(x; X = zeros(length(x), $tr_length), param=nothing, dom_size)
        $(CN)tr_in(Tuple($tr), X, x, param)
        X[:, 1] .= 1:length(x) .|> (i -> $ar(@view X[i, 1:$tr_length]))
        return $ag(@view X[:, 1]) |> (y -> $co(y; param, dom_size, nvars=length(x)))
    end
    """
    return documentation * format_text(output, BlueStyle(); pipe_to_function_call=false)
end

function composition_to_file!(c::Composition, path, name, language=:Julia)
    output = code(c, language; name)
    write(path, output)
    return nothing
end


text:
module GridsTests

using Test
using Gridap
using Gridap.CellValuesGallery
using Gridap.Grids: GridFromData

c2v = [[1,2,4,5],[2,3,5,6]]
v2x = Point{2,Float64}[(0.,0.),(0.,1.),(0.,2.),(1.,0.),(1.,1.),(1.,2.)]

c2v = CellValueFromArray(c2v)
v2x = CellValueFromArray(v2x)

l = length(c2v)
order = 1
t = (HEX_AXIS,HEX_AXIS)
c2t = ConstantCellValue(t,l)
c2o = ConstantCellValue(order,l)

grid = GridFromData(v2x,c2v,c2t,c2o)

test_grid(grid,6,2)

polys = CellPolytopes(grid)

@test length(polys) == 2
@test isa(polys,ConstantCellValue{Polytope{2}})

@test string(grid) == "GridFromData object"

s = "GridFromData object:\n pointdim: 2\n celldim: 2\n npoints: 6\n ncells: 2"
@test sprint(show,"text/plain",grid) == s

end # module


text:
###convolution layers forprop

function layerForProp!(
    cLayer::Conv1D,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...
)

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    img2col = getindex(kwargs, :img2col) ? kwargs[:img2col] : false
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    n_Hi, ci, m = paddedSize(cLayer, Ai)
    s_H = cLayer.s
    f_H = cLayer.f
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1

    if NNlib
        NNConv!(cLayer, Ai)
    elseif img2col
        ## in case input different size than previous time
        if (n_H * c, n_Hi * ci) != size(cLayer.K)
            cLayer.K = unroll(cLayer, (n_Hi, ci, m))
        end
        img2colConvolve!(cLayer, Ai)
    else
        cLayer.Z = zeros(eltype(Ai), n_H, cLayer.channels, m)
        convolve!(cLayer, Ai)
    end

    cLayer.outputS = size(cLayer.A)

    Ai = nothing
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing

end #function layerForProp!(cLayer::Conv1D)

function layerForProp!(
    cLayer::Conv2D,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...
)
    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    img2col = getindex(kwargs, :img2col) ? kwargs[:img2col] : false
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    n_Hi, n_Wi, ci, m = paddedSize(cLayer, Ai)
    c = cLayer.channels
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1

    if NNlib
        NNConv!(cLayer, Ai)
    elseif img2col
        ## in case input different size than previous time
        if (n_W * n_H * c, n_Hi * n_Wi * ci) != size(cLayer.K)
            cLayer.K = unroll(cLayer, (n_Hi, n_Wi, ci, m))
        end #if (n_W*n_H* c, n_Hi*n_Wi*ci) != size(cLayer.K)
        img2colConvolve!(cLayer, Ai)
    else
        cLayer.Z = zeros(eltype(Ai), n_H, n_W, cLayer.channels, m)
        convolve!(cLayer, Ai)
    end #if img2colConvolve

    cLayer.outputS = size(cLayer.A)

    Ai = nothing
    cLayer.forwCount += 1

    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::Conv2D)

function layerForProp!(
    cLayer::Conv3D,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...
)
    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    img2col = getindex(kwargs, :img2col) ? kwargs[:img2col] : false
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    n_Hi, n_Wi, n_Di, ci, m = paddedSize(cLayer, Ai)
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f
    c = cLayer.channels

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    if NNlib
        NNConv!(cLayer, Ai)
    elseif img2col
        ## in case input different size than previous time
        if (n_W * n_H * n_D * c, n_Hi * n_Wi * n_Di * ci) != size(cLayer.K)
            cLayer.K = unroll(cLayer, (n_Hi, n_Wi, n_Di, ci, m))
        end #if (n_W*n_H* c, n_Hi*n_Wi*ci) != size(cLayer.K)
        img2colConvolve!(cLayer, Ai)
    else
        cLayer.Z = zeros(eltype(Ai), n_H, n_W, n_D, cLayer.channels, m)
        convolve!(cLayer, Ai)
    end #if img2colConvolve

    cLayer.outputS = size(cLayer.A)

    cLayer.forwCount += 1

    Ai = nothing

    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::Conv3D)

### Pooling Layers

#import only the needed parts not to have conflict
import NNlib.maxpool, NNlib.meanpool, NNlib.maxpool!, NNlib.meanpool!, NNlib.PoolDims

function layerForProp!(
    cLayer::OneD,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...
) where {OneD<:Union{MaxPool1D,AveragePool1D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    n_Hi, ci, m = paddedSize(cLayer, Ai)
    s_H = cLayer.s
    f_H = cLayer.f
    c = cLayer.channels

    n_H = (n_Hi - f_H) ÷ s_H + 1

    cLayer.A = zeros(eltype(Ai), n_H, c, m)

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            maxpool!(cLayer.A, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            meanpool!(cLayer.A, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer
    elseif f_H == s_H #to use the built-in reshape and maximum and mean
        fastPooling!(cLayer, Ai)
    else
        pooling!(cLayer, Ai)
    end #if NNlibConv

    cLayer.outputS = size(cLayer.A)
    Ai = nothing
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing
end #unction layerForProp!(cLayer::OneD) where {OneD <: Union{MaxPool1D, AveragePool1D}}

function layerForProp!(
    cLayer::TwoD,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
) where {TwoD<:Union{MaxPool2D,AveragePool2D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    n_Hi, n_Wi, ci, m = paddedSize(cLayer, Ai)
    s_H, s_W = S = cLayer.s
    f_H, f_W = F = cLayer.f
    c = cLayer.channels
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1

    cLayer.A =
        zeros(eltype(Ai), n_H, n_W, c, m)

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            maxpool!(cLayer.A, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            meanpool!(cLayer.A, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer
    elseif F == S #to use the built-in reshape, maximum and mean
        fastPooling!(cLayer, Ai)
    else
        pooling!(cLayer, Ai)
    end #if NNlibConv

    cLayer.outputS = size(cLayer.A)

    Ai = nothing
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing

end #function layerForProp!(cLayer::TwoD) where {TwoD <: Union{MaxPool2D, AveragePool2D}}

function layerForProp!(
    cLayer::ThreeD,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
) where {ThreeD<:Union{MaxPool3D,AveragePool3D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = size(Ai)
    # Ai = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    n_Hi, n_Wi, n_Di, ci, m = paddedSize(cLayer, Ai)
    s_H, s_W, s_D = S = cLayer.s
    f_H, f_W, f_D = F = cLayer.f
    c = cLayer.channels

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    cLayer.A =
        zeros(eltype(Ai), n_H, n_W, n_D, c, m)

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            maxpool!(cLayer.A, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            meanpool!(cLayer.A, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer
    elseif F == S #to use the built-in reshape and maximum and mean
        fastPooling!(cLayer, Ai)
    else
        pooling!(cLayer, Ai)
    end #if NNlibConv

    cLayer.outputS = size(cLayer.A)

    Ai = nothing
    cLayer.forwCount += 1

    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::ThreeD) where {ThreeD <: Union{MaxPool3D, AveragePool3D}}


text:
# Autogenerated wrapper script for qwtw_jll for x86_64-w64-mingw32-cxx11
export qwproc, qwtw

using Libglvnd_jll
using boost_jll
using Qt_jll
using qwt_jll
using CompilerSupportLibraries_jll
using FreeType2_jll
using OpenSSL_jll
using marble_jll
using MathGL_jll
JLLWrappers.@generate_wrapper_header("qwtw")
JLLWrappers.@declare_executable_product(qwproc)
JLLWrappers.@declare_library_product(qwtw, "libqwtw.dll")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, boost_jll, Qt_jll, qwt_jll, CompilerSupportLibraries_jll, FreeType2_jll, OpenSSL_jll, marble_jll, MathGL_jll)
    JLLWrappers.@init_executable_product(
        qwproc,
        "bin\\qwproc.exe",
    )

    JLLWrappers.@init_library_product(
        qwtw,
        "bin\\libqwtw.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const SRect = Rect{2, Cshort}
const SPoint = Point{2, Cshort}

struct FaceLocationRaw
    confidence::Cshort
    boundingbox::SRect
    points::NTuple{5, SPoint}
    padding::NTuple{127, Cshort}
end

"""
    confidence: value between 0..1
    boundingbox: axis aligned boundingbox of the face
    points: 5 points, 2xeye, nose, 2xmouth
"""
struct FaceLocation
    confidence::Float64
    boundingbox::Rect{2, Float32}
    points::Vector{Point2f0}
end

convert_point(rows, p) = Point2f0(p[1], rows - p[2])
getpoints(rows, face) = convert_point.(rows, collect(face.points))

function convert_rect(rows, face)
    r = face.boundingbox
    o = origin(r)
    wh = widths(r)
    return Rect(Vec(convert_point(rows, o)) - Vec2f0(0, wh[2]), Vec2f0(wh))
end

function convert_location(rows::Int, location::FaceLocationRaw)
    return FaceLocation(location.confidence/100.0,
                        convert_rect(rows, location),
                        getpoints(rows, location))
end

function detect_faces(img)
    buffer = fill(Cuchar(0), 0x20000)
    img_conv = BGR{N0f8}.(permutedims(img))
    # buffer needs to be preserved outside ccall since i think
    # it's used for the return value
    GC.@preserve buffer begin
        result = ccall((:_Z14facedetect_cnnPhS_iii, libfacedetection), Ptr{Cint}, (Ptr{Cuchar}, Ptr{Cuchar}, Cint, Cint, Cint),
                        buffer, img_conv, size(img, 2), size(img, 1), size(img, 2)*3)
        result == C_NULL && return nothing
        n_faces = unsafe_load(result)
        faces = unsafe_wrap(Array, Ptr{FaceLocationRaw}(Ptr{Cshort}(result) + 4), n_faces)
        return convert_location.(size(img, 1), faces)
    end
end


text:
struct CheapestRute
    rute::Vector{Int};
    rute_total_cost;
    rute_intervals_cost::Vector;
end

import Base.show;
function Base.show(io::IO, cr::CheapestRute)
    println("CheapestRute");
    println("rute: ", cr.rute);
    println("rute_total_cost: ", cr.rute_total_cost);
    println("rute_intervals_cost: ", cr.rute_intervals_cost);
end


text:
module ImageUtils

using ImageMagick
using Reexport
#using FileIO
#import FileIO.save
using ColorVectorSpace
using AxisArrays
const axes = Base.axes
using Random
using Statistics
using LinearAlgebra
using Interpolations
using ImageCore
using Unitful
using FFTW

@reexport using Colors
@reexport using FixedPointNumbers
@reexport using Images
@reexport using ImageMetadata
@reexport using ImageAxes

include("Arrays.jl")
include("Coloring.jl")
include("Slicing.jl")
include("Transformation.jl")
include("PerfusionMaps.jl")
include("FileExport.jl")
include("Permutation.jl")
include("AnalyzeNifti.jl")

end # module


text:
using MyExample
using Test

@testset "Function tests" begin
    #2x+3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
    @test my_f(0,3) == 9
end

@testset "Derivative tests" begin
    @test derivative_of_my_f(2,1) == 2
end