text:
"""
    interp(spl::CubicSpline, xrange)

Returns tuple of vectors of spline values `s` at locations `xs`,
as defined by the grid points at `x+n*dλ` (smoothed curve).

# Examples:

for x=[0,10,30,50,70,90,100] and y=[30,130,150,150,170,220,320]

```jldoctest
julia> cs = spline3(x,y,:natural)
CubicSpline([-0.015813397129186602, 0.009126794258373205, -0.0006937799043062205, -0.0013516746411483252, 0.007350478468899521, -0.01305023923444976], [0.0, -0.4744019138755981, 0.07320574162679427, 0.03157894736842105, -0.049521531100478466, 0.3915071770334928], [11.58133971291866, 6.83732057416268, -1.1866028708133973, 0.9090909090909092, 0.5502392344497609, 7.389952153110048], [30.0, 130.0, 150.0, 150.0, 170.0, 220.0], [0.0, 10.0, 30.0, 50.0, 70.0, 90.0, 100.0])

julia> interp(cs, 30.0:0.2:320.0)
([30.0, 30.2, 30.4, 30.6, 30.8, 31.0, 31.2, 31.4, 31.6, 31.8  …  318.2, 318.4, 318.6, 318.8, 319.0, 319.2, 319.4, 319.6, 319.8, 320.0], [150.00000000000003, 149.76560210526316, 149.53702736842106, 149.31424248803827, 149.09721416267942, 148.88590909090908, 148.68029397129186, 148.48033550239234, 148.28600038277511, 148.09725531100477  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function interpolate(spl::CubicSpline, xrange)
    xs=collect(xrange)
    s=zeros(length(xs))         # spline spectral vector
    j=1                         # count through knots
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while (xs[i] > spl.x[j+1])
                j+=1
            end
            delta = xs[i]-spl.x[j]
            delta2 = delta * delta
            s[i] = spl.d[j] * delta2 * delta + spl.c[j] * delta2 + spl.b[j] * delta + spl.a[j]
        end
    end
    return xs,s
end

"""
    interpolate(spl::HermiteSpline, xrange)

Create interpolated vectors based on a Hermite spline

# Examples:

for x=[0,10,30,50,70,90,100] and y=[30,130,150,150,170,220,320]

```jldoctest
julia> hs = hermitespline(x,y,:kruger)
HermiteSpline([14.909090909090908, 0.18181818181818182, 0.0, 0.0, 0.5714285714285714, 0.16, 14.92], [30.0, 130.0, 150.0, 150.0, 170.0, 220.0, 320.0], [0.0, 10.0, 30.0, 50.0, 70.0, 90.0, 100.0])

julia> interp(hs, 30.0:0.2:320.0)
([30.0, 30.2, 30.4, 30.6, 30.8, 31.0, 31.2, 31.4, 31.6, 31.8  …  318.2, 318.4, 318.6, 318.8, 319.0, 319.2, 319.4, 319.6, 319.8, 320.0], [150.0, 150.0, 150.00000000000003, 149.99999999999997, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function interpolate(spl::HermiteSpline, xrange)
    xs = collect(xrange)
    s = Vector{Float64}(undef, length(xs))
    j = 1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while (xs[i] > spl.x[j+1])
                j+=1
            end
            #scaling and translation factor
            xb = (xs[i] - spl.x[j]) / (spl.x[j+1] - spl.x[j])
            xb2 = xb * xb
            #basis functions for function values
            Φ0 = 2.0 * xb2 * xb - 3.0 * xb2 + 1.0
            Φ1 = xb2 * xb - 2.0 * xb2 + xb
            #basis functions for derivatives
            Ψ0 = -2.0 * xb2 * xb + 3.0 * xb2
            Ψ1 = xb2 * xb - xb2
            Δx = spl.x[j+1] - spl.x[j]

            s[i] = spl.a[j] * Φ0 + Δx * spl.b[j] * Φ1 + spl.a[j+1] * Ψ0 + Δx * spl.b[j+1] * Ψ1
        end
    end
    return xs, s
end

"""
    interpolate(spl::QuadraticSpline, xrange)

Create interpolated vectors based on a quadratic spline

# Examples:

for x=[0,10,30,50,70,90,100] and y=[30,130,150,150,170,220,320]

```jldoctest
julia> qs = spline2(x,y)
QuadraticSpline([1.0, -0.95, 0.9, -0.85, 0.925, -1.1, 0.0], [0.0, 20.0, -18.0, 18.0, -16.0, 21.0, -1.0], [30.0, 130.0, 150.0, 150.0, 170.0, 220.0, 0.0], [0.0, 10.0, 30.0, 50.0, 70.0, 90.0, 100.0])

julia> interp(qs, 30.0:0.2:320.0)
([30.0, 30.2, 30.4, 30.6, 30.8, 31.0, 31.2, 31.4, 31.6, 31.8  …  318.2, 318.4, 318.6, 318.8, 319.0, 319.2, 319.4, 319.6, 319.8, 320.0], [150.0, 146.436, 142.94400000000002, 139.52399999999997, 136.176, 132.9, 129.69600000000003, 126.56400000000002, 123.50399999999998, 120.51599999999999  …  2.2006465129835e-310, 2.2034094153224e-310, 2.2061740509752e-310, 2.20894041994186e-310, 2.2117085222224e-310, 2.21447835781683e-310, 2.2172499267252e-310, 2.2200232289474e-310, 2.22279826448355e-310, 2.22557503333353e-310])
```
"""
function interpolate(spl::QuadraticSpline, xrange)
    xs = collect(xrange)
    s = zeros(length(xs))         # spline spectral vector
    j = 1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            delta = xs[i] - spl.x[j]
            s[i] = spl.c[j] * delta * delta + spl.b[j] * delta + spl.a[j]
        end
    end
    return xs,s
end

"""
    interpolate(spl::LinearSpline, xrange)

Create interpolated vectors based on a linear spline

# Examples:

for x=[0,10,30,50,70,90,100] and y=[30,130,150,150,170,220,320]

```jldoctest
julia> ls = spline1(x,y)
LinearSpline([10.0, 1.0, 0.0, 1.0, 2.5, 10.0], [30.0, 130.0, 150.0, 150.0, 170.0, 220.0], [0.0, 10.0, 30.0, 50.0, 70.0, 90.0, 100.0])

julia> interp(ls, 30.0:0.2:320.0)
([30.0, 30.2, 30.4, 30.6, 30.8, 31.0, 31.2, 31.4, 31.6, 31.8  …  318.2, 318.4, 318.6, 318.8, 319.0, 319.2, 319.4, 319.6, 319.8, 320.0], [150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function interpolate(spl::LinearSpline, xrange)
    xs = collect(xrange)
    s = zeros(length(xs))         # spline spectral vector
    j = 1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            s[i] = spl.b[j] * (xs[i] - spl.x[j]) + spl.a[j]
        end
    end
    return xs,s
end

"""
    interpolate(spl::NearestNeighborSpline, xrange)

Create interpolated vectors based on a nearest neighbor spline.

# Examples:

for x=[0,10,30,50,70,90,100] and y=[30,130,150,150,170,220,320]

```jldoctest
julia> nns = spline0(x,y)
NearestNeighborSpline([30.0, 130.0, 150.0, 150.0, 170.0, 220.0, 320.0], [0.0, 10.0, 30.0, 50.0, 70.0, 90.0, 100.0])

julia> interp(nns, 30.0:0.2:320.0)
([30.0, 30.2, 30.4, 30.6, 30.8, 31.0, 31.2, 31.4, 31.6, 31.8  …  318.2, 318.4, 318.6, 318.8, 319.0, 319.2, 319.4, 319.6, 319.8, 320.0], [150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function interpolate(spl::NearestNeighborSpline, xrange)
    xs = collect(xrange)
    s = zeros(length(xs))
    j = 1
    @inbounds for i in 1:length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            if xs[i] == spl.x[end]
                s[i] = spl.a[end]
            elseif 2 * xs[i] < spl.x[j] + spl.x[j+1]
                s[i]=spl.a[j]
            else
                s[i]=spl.a[j+1]
            end
        end
    end
    return xs, s
end

# alias
const interp = interpolate


text:
# Creating default envs
function create_env(; seed=37, env=:simple, kwargs...)
    if env === :single
        dt = 1.0
        instance_cost = 1.0
        jobtypes = JobParams[
            JobParams(
                path = collect(1:3),
                time = ones(3),
                deadline = 3.5,
                value = 3 * (instance_cost + 1) + 10,
                arrival = FlippingArrival(dt/10, 0:3),
            ),
        ]
        ServiceMeshEnv(;
            kwargs...,
            microservices = 3,
            jobtypes = jobtypes,
            min_scale = 1,
            max_scale = 10,
            max_queue = 5,
            close_time = 0.0,
            boot_time = 1.0,
            dt = dt,
            instance_cost = instance_cost,
            seed = seed,
        )
    elseif env === :double
        dt = 1.0
        instance_cost = 1.0
        jobtypes = JobParams[
            JobParams(
                path = [1, 2, 4],
                time = Float64[1, 1, 1],
                deadline = 3.5,
                value = 3 * (instance_cost + 1) + 10,
                arrival = FlippingArrival(dt/10, 0:4),
            ),
            JobParams(
                path = [1, 3, 4],
                time = Float64[1, 1, 1],
                deadline = 3.5,
                value = 3 * (instance_cost + 1) + 10,
                arrival = FlippingArrival(dt/10, 0:2),
            ),
        ]
        ServiceMeshEnv(;
            kwargs...,
            microservices = 4,
            jobtypes = jobtypes,
            min_scale = 1,
            max_scale = 10,
            max_queue = 5,
            close_time = 0.0,
            boot_time = 1.0,
            dt = dt,
            instance_cost = instance_cost,
            seed = seed,
        )
    else
        throw(ArgumentError("no env for that tag."))
    end
end

# Wrap env in common RL wrappers
function rl_wrappers(env; frames=1, reward_scaling=0.001, state_smoothing=0.0, kwargs...)
    A = action_space(env)
    alow = [x isa UnitRange ? x.start : x.left for x in A]
    ahigh = [x isa UnitRange ? x.stop : x.right for x in A]
    S = state_space(env)
    slow = [x.left for x in S]
    shigh = [x.right for x in S]

    wrapped_env = env

    wrapped_env = ActionTransformedEnv(
        StateTransformedEnv(
            RewardTransformedEnv(wrapped_env, x -> x * reward_scaling);
            state_mapping = x -> 2 .* (x .- slow) ./ (shigh .- slow) .- 1,
            state_space_mapping = x -> Space([-1f0..1f0 for _ in x])
        );
        action_mapping = x -> round.(Int, clamp.(alow .+ (x .+ 1) .* 0.5 .* (ahigh .- alow), alow, ahigh)),
        action_space_mapping = x -> Space([-1f0..1f0 for _ in x])
    )

    if frames != 1
        wrapped_env = StateStackedEnv(wrapped_env; n=frames)
    elseif state_smoothing != 0
        wrapped_env = SmoothedStateEnv(wrapped_env; smoothing=state_smoothing)
    end

    return wrapped_env
end

# Wrap env in common wrappers
function default_wrappers(env; closing_cost=0.0, kwargs...)
    if closing_cost != 0
        env = RewardOverriddenEnv(env,
            function(env)
                num_closing = sum(ms->length(ms.closing_nodes), env.microservices)
                reward(env) - closing_cost * num_closing
            end
        )
    end
    return env
end

# StateStacked wrapper
mutable struct StateStackedEnv{E<:AbstractEnv, T} <: AbstractEnvWrapper
    env::E
    state::Matrix{T}
    ptr::Int
end

StateStackedEnv(env; n=1) =
    StateStackedEnv(env, repeat(state(env), 1, n), 1)

function (env::StateStackedEnv)(a)
    env[](a)
    env.ptr %= size(env.state, 2)
    env.ptr += 1
    env.state[:, env.ptr] .= state(env[])
    nothing
end

RLBase.state(env::StateStackedEnv) = vcat(
    reshape(env.state[:, (env.ptr+1):end], :),
    reshape(env.state[:, 1:env.ptr], :)
)

RLBase.state_space(env::StateStackedEnv, args...; kwargs...) =
    Space(vcat((state_space(env.env, args...; kwargs...).s for i in 1:size(env.state, 2))...))

# RunningState wrapper
mutable struct SmoothedStateEnv{E<:AbstractEnv, T} <: AbstractEnvWrapper
    env::E
    running_state::Vector{T}
    smoothing::Float64
end

SmoothedStateEnv(env; smoothing) =
    SmoothedStateEnv(env, state(env), smoothing)

function (env::SmoothedStateEnv)(a)
    env[](a)
    env.running_state .= (1 - env.smoothing) * env.running_state + env.smoothing * state(env[])
    nothing
end

RLBase.state(env::SmoothedStateEnv) = vcat(
    env.running_state,
    state(env[])
)

RLBase.state_space(env::SmoothedStateEnv, args...; kwargs...) =
    Space(vcat((state_space(env.env, args...; kwargs...).s for i in 1:2)...))

ext: jl

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
"""
Grammar-based Monte Carlo Tree Search algorithm
Modifications: max instead of mean. Number of actions varying with state. Deterministic assumed.
Subtree explorations are tracked.
"""
module GBMCTS

using POMDPs

export
    MCTSSolver,
    MCTSPolicy,
    action,
    simulate,
    rollout,
    # SPW
    StateNode

typealias Reward Float64

include("policies.jl")
include("vanilla.jl")

end # module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
const CRLF = "\r\n"

const Byte = UInt8
const Bytes = Vector{UInt8}

"zero_index(0) == zero_index(1), not sure if it's OK"
zero_index(x) = ifelse(x > 0, x - 1, x)

"use with caution"
function chomp!(s::Bytes)
    # @assert s[end-1] == '\r' && s[end] == '\n'
    ccall(:jl_array_del_end, Nothing, (Any, UInt), s, 2)
    return s
end

car(x) = x[1]
cdr(x) = x[2:end]
cadr(x) = x[2]

<<(x::IO, y) = (print(x, y); x)
<<(x::IO, y::Byte) = (write(x, y); x)
<<(x::IO, y::Bytes) = (write(x, y); x)
<<(x::IO, f::Function) = (f(x); x)

>>(x::IO, y) = read(x, y)
>>(x::IO, f::Function) = f(x)

macro some(x)
    quote
        s = $(esc(x))
        s == nothing && return nothing
        s
    end
end


text:
"""
    Constraints(mode::String)
    Constraints(;
        linear_acceleration::AbstractArray{<:Integer,1} = [-IPVACX,IPVACY,IPVACZ],
        angular_acceleration::AbstractArray{<:Integer,1} = -[IPRACX,IPRACY,IPRACZ],
        velocity::AbstractArray{<:Integer,1} = -[IPVIAS,IPBETA,IPALFA],
        rotation_rate::AbstractArray{<:Integer,1} = -[IPROTX,IPROTY,IPROTZ],
        position::AbstractArray{<:Integer,1} = -[IPPOSX, IPPOSY, IPPOSZ],
        phi::Integer = -IPBANK,
        theta::Integer = -IPELEV,
        psi::Integer = -IPHEAD,
        flap_defl_ctrl_var::AbstractArray{<:Integer,1} =
            -(IPFLP1:(IPFLP1-1+NFLPX)),
        eng_pwr_ctrl_var::AbstractArray{<:Integer,1} =
            -(IPENG1:(IPENG1-1+NENGX))
        err_int_Vinf::Integer = IPVIAS,
        err_int_beta::Integer = IPBETA,
        err_int_alpha::Integer = IPALFA,
        err_int_phi::Integer = IPROTX,
        err_int_theta::Integer = IPROTY,
        err_int_psi::Integer = IPROTZ,
        err_int_ROTxyz::AbstractArray{<:Integer,1} = [0,0,0],
        err_int_VACxyz::AbstractArray{<:Integer,1} = [0,0,0])

    Defines operating point constraints as done within ASWING.
    Implemented modes: ["default", "anchored", "free", "static"]
"""
struct Constraints
    linear_acceleration::Array{Int,1}
    angular_acceleration::Array{Int,1}
    velocity::Array{Int,1}
    rotation_rate::Array{Int,1}
    position::Array{Int,1}
    phi::Int
    theta::Int
    psi::Int
    flap_defl_ctrl_var::Array{Int,1}
    eng_pwr_ctrl_var::Array{Int,1}
    err_int_Vinf::Int
    err_int_beta::Int
    err_int_alpha::Int
    err_int_phi::Int
    err_int_theta::Int
    err_int_psi::Int
    err_int_ROTxyz::Array{Int,1}
    err_int_VACxyz::Array{Int,1}
end

function Constraints(;
    linear_acceleration::AbstractArray{<:Integer,1} = -[IPVACX,IPVACY,IPVACZ],
    angular_acceleration::AbstractArray{<:Integer,1} = -[IPRACX,IPRACY,IPRACZ],
    velocity::AbstractArray{<:Integer,1} = -[IPVIAS,IPBETA,IPALFA],
    rotation_rate::AbstractArray{<:Integer,1} = -[IPROTX,IPROTY,IPROTZ],
    position::AbstractArray{<:Integer,1} = -[IPPOSX,IPPOSY,IPPOSZ],
    phi::Integer = -IPBANK,
    theta::Integer = -IPELEV,
    psi::Integer = -IPHEAD,
    flap_defl_ctrl_var::AbstractArray{<:Integer,1} =
        -(IPFLP1:(IPFLP1-1+NFLPX)),
    eng_pwr_ctrl_var::AbstractArray{<:Integer,1} =
        -(IPENG1:(IPENG1-1+NENGX)),
    err_int_Vinf::Integer = IPVIAS,
    err_int_beta::Integer = IPBETA,
    err_int_alpha::Integer = IPALFA,
    err_int_phi::Integer = IPROTX,
    err_int_theta::Integer = IPROTY,
    err_int_psi::Integer = IPROTZ,
    err_int_ROTxyz::AbstractArray{<:Integer,1} = [0,0,0],
    err_int_VACxyz::AbstractArray{<:Integer,1} = [0,0,0]) where Int<:Integer
    Constraints(linear_acceleration, angular_acceleration, velocity,
        rotation_rate, position, phi, theta, psi, flap_defl_ctrl_var,
        eng_pwr_ctrl_var, err_int_Vinf, err_int_beta, err_int_alpha,
        err_int_phi, err_int_theta, err_int_psi, err_int_ROTxyz, err_int_VACxyz)
end

function Constraints(mode::String) where Int<:Integer
    if lowercase(mode) == "default"
        linear_acceleration = [-IPVACX,IPVACY,IPVACZ]
        angular_acceleration = -[IPRACX,IPRACY,IPRACZ]
        velocity = -[IPVIAS,IPBETA,IPALFA]
        rotation_rate = -[IPROTX,IPROTY,IPROTZ]
        position = -[IPPOSX,IPPOSY,IPPOSZ]
        phi = -IPBANK
        theta = -IPELEV
        psi = -IPHEAD
        flap_defl_ctrl_var = -(IPFLP1:(IPFLP1-1+NFLPX))
        eng_pwr_ctrl_var = -(IPENG1:(IPENG1-1+NENGX))
    elseif lowercase(mode) == "anchored"
        linear_acceleration = [IPVACX,IPVACY,IPVACZ]
        angular_acceleration = [IPRACX,IPRACY,IPRACZ]
        velocity = [IPVIAS,IPBETA,IPALFA]
        rotation_rate = [IPROTX,IPROTY,IPROTZ]
        position = -[IPPOSX,IPPOSY,IPPOSZ]
        phi = -IPBANK
        theta = -IPELEV
        psi = -IPHEAD
        flap_defl_ctrl_var = IPFLP1:(IPFLP1-1+NFLPX)
        eng_pwr_ctrl_var = IPENG1:(IPENG1-1+NENGX)
    elseif lowercase(mode) == "free"
        linear_acceleration = [IPFORX,IPFORY,IPFORZ]
        angular_acceleration = [IPMOMX,IPMOMY,IPMOMZ]
        velocity = [IPVIAS,IPBETA,IPALFA]
        rotation_rate = [IPROTX,IPROTY,IPROTZ]
        position = [IPPOSX,IPPOSY,IPPOSZ]
        phi = IPBANK
        theta = IPELEV
        psi = IPHEAD
        flap_defl_ctrl_var = IPFLP1:(IPFLP1-1+NFLPX)
        eng_pwr_ctrl_var = IPENG1:(IPENG1-1+NENGX)
    elseif lowercase(mode) == "static"
        linear_acceleration = [IPVACX,IPVACY,IPVACZ]
        angular_acceleration = [IPRACX,IPRACY,IPRACZ]
        velocity = [IPFORX,IPFORY,IPFORZ]
        rotation_rate = [IPMOMX,IPMOMY,IPMOMZ]
        position = -[IPPOSX,IPPOSY,IPPOSZ]
        phi = -IPBANK
        theta = -IPELEV
        psi = -IPHEAD
        flap_defl_ctrl_var = IPFLP1:(IPFLP1-1+NFLPX)
        eng_pwr_ctrl_var = IPENG1:(IPENG1-1+NENGX)
    else
        error("Undefined default constraint set")
    end

    err_int_Vinf = IPVIAS
    err_int_beta = IPBETA
    err_int_alpha = IPALFA
    err_int_phi = IPROTX
    err_int_theta = IPROTY
    err_int_psi = IPROTZ
    err_int_ROTxyz = [0,0,0]
    err_int_VACxyz = [0,0,0]

    Constraints(linear_acceleration, angular_acceleration, velocity,
        rotation_rate, position, phi, theta, psi, flap_defl_ctrl_var,
        eng_pwr_ctrl_var, err_int_Vinf, err_int_beta, err_int_alpha,
        err_int_phi, err_int_theta, err_int_psi, err_int_ROTxyz, err_int_VACxyz)
end

"""
    Parameters(;
        linear_acceleration::AbstractArray{<:Real,1} = zeros(Float64,3),
        angular_acceleration::AbstractArray{<:Real,1} = zeros(Float64,3),
        velocity::Real = 0.0,
        beta::Real = 0.0,
        alpha::Real = 0.0,
        rotation_rate::AbstractArray{<:Real,1} = zeros(Float64,3),
        position::AbstractArray{<:Real,1} = zeros(Float64,3),
        phi::Real = 0.0,
        theta::Real = 0.0,
        psi::Real = 0.0,
        flap_defl_ctrl_var::AbstractArray{<:Real,1} = zeros(Float64,NFLPX),
        eng_pwr_ctrl_var::AbstractArray{<:Real,1} = zeros(Float64,NENGX),
        sum_force::AbstractArray{<:Real,1} = zeros(Float64,3),
        sum_mom::AbstractArray{<:Real,1} = zeros(Float64,3),
        lift::Real = 0.0,
        climb_angle::Real = 0.0,
        radial_acceleration::Real = 0.0,
        usr1::Real = 0.0,
        usr2::Real = 0.0,
        least_squared::Real = 0.0)

    Defines operating point parameters as done within ASWING.
"""
struct Parameters
    linear_acceleration::Array{Float64,1}
    angular_acceleration::Array{Float64,1}
    velocity::Float64
    beta::Float64
    alpha::Float64
    rotation_rate::Array{Float64,1}
    position::Array{Float64,1}
    phi::Float64
    theta::Float64
    psi::Float64
    flap_defl_ctrl_var::Array{Float64,1}
    eng_pwr_ctrl_var::Array{Float64,1}
    sum_force::Array{Float64,1}
    sum_mom::Array{Float64,1}
    lift::Float64
    climb_angle::Float64
    radial_acceleration::Float64
    usr1::Float64
    usr2::Float64
    least_squared::Float64
    function Parameters(
        linear_acceleration::AbstractArray{<:Real,1},
        angular_acceleration::AbstractArray{<:Real,1},
        velocity::Real, beta::Real, alpha::Real,
        rotation_rate::AbstractArray{<:Real,1},
        position::AbstractArray{<:Real,1},
        phi::Real, theta::Real, psi::Real,
        flap_defl_ctrl_var::AbstractArray{<:Real,1},
        eng_pwr_ctrl_var::AbstractArray{<:Real,1},
        sum_force::AbstractArray{<:Real,1},
        sum_mom::AbstractArray{<:Real,1},
        lift::Real, climb_angle::Real, radial_acceleration::Real,
        usr1::Real, usr2::Real, least_squared::Real)

        if ((length(linear_acceleration) != 3) || (length(angular_acceleration) != 3) ||
            (length(rotation_rate) != 3) ||
            (length(position) != 3) || (length(flap_defl_ctrl_var) > NFLPX) ||
            (length(eng_pwr_ctrl_var) > NENGX) || (length(sum_force) != 3) ||
            (length(sum_mom) != 3))
            error("One or more input arrays are not the correct dimensions")
        end

        new(linear_acceleration, angular_acceleration, velocity, beta, alpha,
            rotation_rate, position, phi, theta, psi, flap_defl_ctrl_var,
            eng_pwr_ctrl_var, sum_force, sum_mom,lift, climb_angle,
            radial_acceleration, usr1, usr2, least_squared)
    end
end

function Parameters(;
    linear_acceleration::AbstractArray{<:Real,1} = zeros(Float64,3),
    angular_acceleration::AbstractArray{<:Real,1} = zeros(Float64,3),
    velocity::Real = 0.0,
    beta::Real = 0.0,
    alpha::Real = 0.0,
    rotation_rate::AbstractArray{<:Real,1} = zeros(Float64,3),
    position::AbstractArray{<:Real,1} = zeros(Float64,3),
    phi::Real = 0.0,
    theta::Real = 0.0,
    psi::Real = 0.0,
    flap_defl_ctrl_var::AbstractArray{<:Real,1} = zeros(Float64,NFLPX),
    eng_pwr_ctrl_var::AbstractArray{<:Real,1} = zeros(Float64,NENGX),
    sum_force::AbstractArray{<:Real,1} = zeros(Float64,3),
    sum_mom::AbstractArray{<:Real,1} = zeros(Float64,3),
    lift::Real = 0.0,
    climb_angle::Real = 0.0,
    radial_acceleration::Real = 0.0,
    usr1::Real = 0.0,
    usr2::Real = 0.0,
    least_squared::Real = 0.0)

    if ((length(linear_acceleration) != 3) || (length(angular_acceleration) != 3) ||
        (length(rotation_rate) != 3) ||
        (length(position) != 3) || (length(flap_defl_ctrl_var) > NFLPX) ||
        (length(eng_pwr_ctrl_var) > NENGX) || (length(sum_force) != 3) ||
        (length(sum_mom) != 3))
        error("One or more input arrays are not the correct dimensions")
    end

    Parameters(linear_acceleration, angular_acceleration, velocity, beta, alpha,
        rotation_rate, position, phi, theta, psi, flap_defl_ctrl_var,
        eng_pwr_ctrl_var, sum_force, sum_mom,lift, climb_angle,
        radial_acceleration, usr1, usr2, least_squared)
end

"""
    OperatingPoint(;
        constraints::Constraints=Constraints(),
        mach_from_airspeed::Bool = false,
        machpg::Real = 0.0,
        ground_image::Real = 0,
        parameters::Parameters = Parameters(),
        wflap::AbstractArray{<:Real,1} = ones(Float64,NFLPX),
        wpeng::AbstractArray{<:Real,1} = ones(Float64,NENGX))

    contains operating point inputs
"""
struct OperatingPoint
    constraints::Constraints # Solution constraints
    mach_from_airspeed::Bool # Set Prandtl-Glauert Mach based on airspeed?
    machpg::Float64 # Prandtl-Glauert Mach
    ground_image::Int # ground image flag, 0 no image,+1 solid ground, -1 free surface (anti-image)
    parameters::Parameters
    wflap::Array{Float64,1} # flap weights for least squared constraint
    wpeng::Array{Float64,1} # engine weights for least squared constraint
end

function OperatingPoint(;
    constraints::Constraints=Constraints(),
    mach_from_airspeed::Bool = false,
    machpg::Real = 0.0,
    ground_image::Real = 0,
    parameters::Parameters = Parameters(),
    wflap::AbstractArray{<:Real,1} = ones(Float64,NFLPX),
    wpeng::AbstractArray{<:Real,1} = ones(Float64,NENGX))

    OperatingPoint(constraints, mach_from_airspeed, machpg, ground_image,
        parameters, wflap, wpeng)
end

function OperatingPoint(;
    constraints::Constraints=Constraints(),
    mach_from_airspeed::Bool = false,
    machpg::Real = 0.0,
    ground_image::Real = 0,
    parameters::Parameters = Parameters(),
    wflap::AbstractArray{<:Real,1} = ones(Float64,NFLPX),
    wpeng::AbstractArray{<:Real,1} = ones(Float64,NENGX))

    OperatingPoint(constraints, mach_from_airspeed, machpg, ground_image,
        parameters, wflap, wpeng)
end


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

module GaussianSimulation

using GeoStatsBase
using Variography
using KrigingEstimators
using Distributions
using LinearAlgebra
using Statistics
using FFTW
using CpuId

import GeoStatsBase: preprocess, solvesingle

include("lu.jl")
include("fft.jl")
include("seq.jl")

export LUGS, FFTGS, SGS

end


text:
mutable struct ARPACKAlloc{T}

    n::Int
    nev::Int
    ncv::Int
    maxiter::Int
    bmat::String
    which::String
    mode::Int
    A::Function
    Amat::Symmetric{T,Matrix{T}}
    lworkl::Int
    TOL::Vector{T}
    v::Matrix{T}
    workd::Vector{T}
    workl::Vector{T}
    rwork::Vector{T}
    resid::Vector{T}
    info::Vector{BlasInt}
    iparam::Vector{BlasInt}
    ipntr::Vector{BlasInt}
    ido::Vector{BlasInt}
    zernm1::UnitRange{Int}
    howmny::String
    select::Vector{BlasInt}
    info_e::Vector{BlasInt}
    d::Vector{T}
    sigmar::Vector{T}
    converged::Bool
    arpackerror::Bool
    x::Vector{T}
    y::Vector{T}
    tol::T

    rng::Random.MersenneTwister

    function ARPACKAlloc{T}() where T
        new{T}()
    end
end

hasconverged(arc::ARPACKAlloc) = arc.converged

function unsafe_getvalues(arc::ARPACKAlloc)
    return arc.d
end

function unsafe_getvectors(arc::ARPACKAlloc)
    return arc.v
end

function ARPACKAlloc(T::DataType, n::Int64)::ARPACKAlloc
    arc = ARPACKAlloc{T}()
    @timeit "init_arc" _init_arc!(arc::ARPACKAlloc, Symmetric(Matrix{T}(I, n, n), :U), 1, n)
    return arc
end

function _init_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64) where T

    # Dimension of the eigenproblem
    arc.n = n
    # Number of eigenvalues of OP to be computed
    arc.nev = nev
    # Standard eigenvalue problem
    arc.bmat = "I"
    # Compute the NEV largest (algebraic) eigenvalues
    arc.which = "LA"
    # Flag for A*x = lambda*x
    arc.mode = 1

    # How many Lanczos vectors are generated
    # arc.ncv = arc.nev + 1
    arc.ncv = max(20, 2 * arc.nev + 1)
    # Stopping criterion
    arc.tol = 1e-6
    arc.TOL = arc.tol * ones(T, 1)
    # Maximum number of iterations
    arc.maxiter = Int(1e+3)

    # Build linear operator
    arc.Amat = A
    arc.x = Vector{T}(undef, arc.n)
    arc.y = Vector{T}(undef, arc.n)

    # Lanczos basis vectors (output)
    arc.v = Matrix{T}(undef, arc.n, arc.ncv)

    # If info != 0, RESID contains the initial residual vector
    arc.info = ones(BlasInt, 1)
    arc.info_e = ones(BlasInt, 1)
    # Resid contains the initial residual vector
    arc.rng = Random.MersenneTwister(1234)
    arc.resid = rand(arc.rng, arc.n)

    # Workspace
    arc.workd = Vector{T}(undef, 3 * arc.n)
    arc.lworkl = arc.ncv * (arc.ncv + 8)
    arc.workl = Vector{T}(undef, arc.lworkl)
    arc.rwork = Vector{T}()
    arc.iparam = ones(BlasInt, 11)
    arc.iparam[1] = BlasInt(1)
    arc.iparam[3] = BlasInt(arc.maxiter)
    arc.iparam[4] = BlasInt(1)
    arc.iparam[7] = BlasInt(arc.mode)
    arc.ipntr = zeros(BlasInt, 11)
    # IDO must be zero on the first call to dsaupd
    arc.ido = zeros(BlasInt, 1)

    # Parameters for _EUPD! routine
    arc.zernm1 = 0:(arc.n-1)
    arc.howmny = "A"
    arc.select = Vector{BlasInt}(undef, arc.ncv)
    arc.d = Vector{T}(undef, arc.nev)
    arc.sigmar = zeros(T, 1)

    # Flags created for ProxSDP use
    arc.converged = false
    arc.arpackerror = false

    return nothing
end

function _update_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, iter::Int64) where T

    if arc.nev < nev
        # Number of eigenvalues of OP to be computed
        arc.nev = nev
        # How many Lanczos vectors are generated
        # arc.ncv = arc.nev + 1
        arc.ncv = max(20, 2 * arc.nev + 1)
        # Lanczos basis vectors (output)
        arc.v = Matrix{T}(undef, arc.n, arc.ncv)
        # Workspace
        arc.lworkl = arc.ncv * (arc.ncv + 8)
        arc.workl = Vector{T}(undef, arc.lworkl)
        # Parameters for _EUPD! routine
        arc.select = Vector{BlasInt}(undef, arc.ncv)
        arc.d = Vector{T}(undef, arc.nev)
    else
        # Number of eigenvalues of OP to be computed
        arc.nev = nev
        # How many Lanczos vectors are generated
        # arc.ncv = arc.nev + 1
        arc.ncv = max(20, 2 * arc.nev + 1)
        # Lanczos basis vectors (output)
        # arc.v = Matrix{T}(undef, arc.n, arc.ncv)
        # Workspace
        arc.lworkl = arc.ncv * (arc.ncv + 8)
        # arc.workl = Vector{T}(undef, arc.lworkl)
        # Parameters for _EUPD! routine
        # arc.select = Vector{BlasInt}(undef, arc.ncv)
        # arc.d = Vector{T}(undef, arc.nev)
    end

    # Build linear operator
    arc.Amat = A

    # If info != 0, RESID contains the initial residual vector
    arc.info[] = BlasInt(1)#ones(BlasInt, 1)
    arc.info_e[] = BlasInt(1)#ones(BlasInt, 1)

    # Iparam
    # arc.iparam = ones(BlasInt, 11)
    fill!(arc.iparam, BlasInt(1))
    arc.iparam[1] = BlasInt(1)
    arc.iparam[3] = BlasInt(arc.maxiter)
    arc.iparam[4] = BlasInt(1)
    arc.iparam[7] = BlasInt(arc.mode)
    # arc.ipntr = zeros(BlasInt, 11)
    fill!(arc.ipntr, BlasInt(0))
    # IDO must be zero on the first call to dsaupd
    # arc.ido = zeros(BlasInt, 1)
    fill!(arc.ido, BlasInt(0))

    # Stopping criterion
    arc.tol = max(arc.tol / (iter / 10.), 1e-10)
    # arc.TOL = arc.tol * ones(T, 1)
    arc.TOL[] = arc.tol * one(T)

    return nothing
end

function _saupd!(arc::ARPACKAlloc{T})::Nothing where T

    while true
        Arpack.saupd(arc.ido, arc.bmat, arc.n, arc.which, arc.nev, Ref(arc.TOL[1]), arc.resid, arc.ncv, arc.v, arc.n,
        arc.iparam, arc.ipntr, arc.workd, arc.workl, arc.lworkl, arc.info)

        if arc.ido[] == 1
            @inbounds @simd for i in 1:arc.n
                arc.x[i] = arc.workd[i-1+arc.ipntr[1]]
            end
            mul!(arc.y, arc.Amat, arc.x)
            @inbounds @simd for i in 1:arc.n
                 arc.workd[i-1+arc.ipntr[2]] = arc.y[i]
            end

        elseif arc.ido[] == 99
            # In this case, don't call _EUPD! (according to https://help.scilab.org/docs/5.3.3/en_US/dseupd.html)
            break
        else
            arc.converged = false
            arc.arpackerror = true
            break
        end
    end

    # Check if _AUPD has converged properly
    if !(0 <= arc.info[] <= 1)
        arc.converged = false
        arc.arpackerror = true
    else
        arc.converged = true
    end

    return nothing
end

function _seupd!(arc::ARPACKAlloc{T})::Nothing where T

    # Check if _AUPD has converged properly
    if !arc.converged
        arc.converged = false
        return nothing
    end

    Arpack.seupd(true, arc.howmny, arc.select, arc.d, arc.v, arc.n, arc.sigmar,
    arc.bmat, arc.n, arc.which, arc.nev, Ref(arc.TOL[1]), arc.resid, arc.ncv, arc.v, arc.n,
    arc.iparam, arc.ipntr, arc.workd, arc.workl, arc.lworkl, arc.info_e)

    # Check if seupd has converged properly
    if arc.info_e[] != 0
        arc.converged = false
        arc.arpackerror = true
        return nothing
    end

    # Check number of converged eigenvalues (maybe some can still be used?)
    nconv = arc.iparam[5]
    if nconv < arc.nev
        arc.converged = false
        arc.arpackerror = false
        return nothing
    end

    arc.converged = true
    arc.arpackerror = false

    return nothing
end

function eig!(arc::ARPACKAlloc, A::Symmetric{T1,Matrix{T1}}, nev::Integer, iter::Int64, warm_start_eig::Bool)::Nothing where T1

    # Initialize parameters and do memory allocation
    @timeit "update_arc" _update_arc!(arc, A, nev, iter)::Nothing

    # Top level reverse communication interface to solve real double precision symmetric problems.
    @timeit "saupd" _saupd!(arc)::Nothing

    # Post processing routine (eigenvalues and eigenvectors purification)
    @timeit "seupd" _seupd!(arc)::Nothing

    return nothing
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(quadruped)
v1 = zeros(quadruped.nq)

# ## Time
h = 0.01
T = 100

# ## Simulator
s = Simulator(quadruped, T, h=h)

# ## Simulate
simulate!(s, q1, v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)


text:
# composite type (like object)
struct User
    name::String
    email::String
    age::Int64
    option # -> Any
end

println(fieldnames(User))

# construct
u = User(
    "John Doe",
    "john@example.com",
    22,
    "Human"
)
println(u)
# error, fields are immutable
# u.email = "john2@example.com"

mutable struct MutableUser
    name::String
    email::String
    age::Int64
    option # -> Any
end

u2 = MutableUser(
    "John Doe",
    "john@example.com",
    22,
    "Human"
)
println(u2)
u2.email = "john2@example.com"

# inner constructor method
struct UserWithValidation
    name::String
    age::Int64
    UserWithValidation(name, age) = age < 0 ? error("Age out of range") : new(name, age)
end
UserWithValidation("John", 18)
# error
# UserWithValidation("John", -1)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
# # Decaying barotropic QG turbulence over topography
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_decaying_topography.ipynb).
#
# An example of decaying barotropic quasi-geostrophic turbulence over topography.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Plots, Printf, Random, Statistics"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Plots, Printf, Random

using Statistics: mean

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n²
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 2000           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
nothing # hide

# Define the topographic potential vorticity, ``\eta = f_0 h(x, y)/H``. The topography here is
# an elliptical mount at ``(x, y) = (1, 1)``, and an elliptical depression at ``(x, y) = (-1, -1)``.
σx, σy = 0.4, 0.8
topographicPV(x, y) = 3exp(-(x-1)^2/(2σx^2) -(y-1)^2/(2σy^2)) - 2exp(-(x+1)^2/(2σx^2) -(y+1)^2/(2σy^2))
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing a viscosity
# coefficient `ν` leads to the module's default value: `ν=0`. In this example numerical instability
# due to accumulation of enstrophy in high wavenumbers is taken care with the `FilteredTimestepper`
# we picked. Thus, we choose not to do any dealiasing by providing `aliased_fraction=0`.
#
# The topophic PV is prescribed via keyword argument `eta`.
prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, eta=topographicPV,
                                  dt=dt, stepper=stepper, aliased_fraction=0)
nothing # hide

# and define some shortcuts
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# and let's plot the topographic PV. Note that when plotting, we decorate the variable to be
# plotted with `Array()` to make sure it is brought back on the CPU when the variable lives
# on the GPU.
contourf(grid.x, grid.y, Array(params.eta'),
          aspectratio = 1,
            linewidth = 0,
               levels = 10,
                    c = :balance,
                 clim = (-3, 3),
                xlims = (-grid.Lx/2, grid.Lx/2),
                ylims = (-grid.Ly/2, grid.Ly/2),
               xticks = -3:3,
               yticks = -3:3,
               xlabel = "x",
               ylabel = "y",
                title = "topographic PV η=f₀h/H")

# ## Setting initial conditions

# Our initial condition consist of a flow that has power only at wavenumbers with
# ``6 < \frac{L}{2\pi} \sqrt{k_x^2 + k_y^2} < 12`` and initial energy ``E_0``.
# `ArrayType()` function returns the array type appropriate for the device, i.e., `Array` for
# `dev = CPU()` and `CuArray` for `dev = GPU()`.

E₀ = 0.04 # energy of initial condition

K = @. sqrt(grid.Krsq)                             # a 2D array with the total wavenumber

Random.seed!(1234)
qih = ArrayType(dev)(randn(Complex{eltype(grid)}, size(sol)))
@. qih = ifelse(K < 6  * 2π/L, 0, qih)
@. qih = ifelse(K > 12 * 2π/L, 0, qih)
qih *= sqrt(E₀ / SingleLayerQG.energy(qih, vars, params, grid))  # normalize qi to have energy E₀
qi = irfft(qih, grid.nx)

SingleLayerQG.set_q!(prob, qi)
nothing # hide

# Let's plot the initial vorticity field:

p1 = heatmap(x, y, Array(vars.q'),
         aspectratio = 1,
              c = :balance,
           clim = (-8, 8),
          xlims = (-grid.Lx/2, grid.Lx/2),
          ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity ∂v/∂x-∂u/∂y",
     framestyle = :box)

p2 = contourf(x, y, Array(vars.ψ'),
        aspectratio = 1,
             c = :viridis,
        levels = range(-0.25, stop=0.25, length=11),
          clim = (-0.25, 0.25),
         xlims = (-grid.Lx/2, grid.Lx/2),
         ylims = (-grid.Ly/2, grid.Ly/2),
        xticks = -3:3,
        yticks = -3:3,
        xlabel = "x",
        ylabel = "y",
         title = "initial streamfunction ψ",
    framestyle = :box)

layout = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout=layout, size = (800, 360))

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps)
Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files.
filepath = "."
filename = joinpath(filepath, "decayingbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
out = Output(prob, filename, (:sol, get_sol))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction and
# their corresponding zonal mean structure.

function plot_output(prob)
  q = prob.vars.q
  ψ = prob.vars.ψ
  η = prob.params.eta

  pq = heatmap(x, y, Array(q'),
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-6, 6),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ∂v/∂x-∂u/∂y",
        framestyle = :box)

  contour!(pq, x, y, Array(η'),
          levels=0.5:0.5:3,
          lw=2, c=:black, ls=:solid, alpha=0.7)

  contour!(pq, x, y, Array(η'),
          levels=-2:0.5:-0.5,
          lw=2, c=:black, ls=:dash, alpha=0.7)

  pψ = contourf(x, y, Array(ψ'),
       aspectratio = 1,
            legend = false,
                 c = :viridis,
            levels = range(-0.75, stop=0.75, length=31),
              clim = (-0.75, 0.75),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  l = @layout Plots.grid(1, 2)
  p = plot(pq, pψ, layout = l, size = (800, 360))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)

  if j % (1000 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = Array(vars.q)
  p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t)
  p[2][1][:z] = Array(vars.ψ)

  stepforward!(prob, diags, nsubs)
  SingleLayerQG.updatevars!(prob)
end

mp4(anim, "singlelayerqg_decaying_topography.mp4", fps=12)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

include(joinpath(@__DIR__, "config.jl"))
include(joinpath(@__DIR__, "lstm.jl"))

using MXNet

vis_n_layer = 1
vis_seq_len = 2
vis_n_class = 128

lstm  = LSTM(vis_n_layer, vis_seq_len, DIM_HIDDEN, DIM_EMBED, vis_n_class, name=NAME, output_states=true)

open("visualize.dot", "w") do io
  println(io, mx.to_graphviz(lstm))
end
run(pipeline(`dot -Tsvg visualize.dot`, stdout="visualize.svg"))


text:
struct Particle{N,T} <: GeometryPrimitive{N,T}
    position::Point{N,T}
    velocity::Vec{N,T}
end


text:
"""
    gapfill!(x t, fs; m::Bool=true)
    y = gapfill(x, t, fs; m::Bool=true)

Fill gaps in `x`, sampled at `fs`, with gap indices given by `t[:,1]` and
gap lengths in μs given by `t[:,2]`.

Specify `m=false` to fill with NaNs; else, fill with the mean of non-NaN
values in `x`.
"""
function gapfill!(x::Array{T,1}, t::Array{Int64,2}, fs::Float64; m::Bool=true) where T<: Real
  (fs == 0.0 || isempty(x)) && (return x)
  nt = size(t,1)
  ng = nt - (t[nt,2] == 0 ? 2 : 1)
  (ng == 0) && (return x)
  mx = m ? mean(skipmissing(x)) : T(NaN)
  u = max(Int64(20), round(Int64, 0.2*fs))
  Δ = round(Int64, 1.0/(fs*μs))
  nx = length(x) + div(sum(t[2:end,2]), Δ)
  resize!(x, nx)

  w = t_win(t, Δ)
  broadcast!(-, w, w, w[1,1])
  broadcast!(div, w, w, Δ)
  broadcast!(+, w, w, 1)
  nw = size(w,1)
  for i = nw:-1:1
    N = w[i,2]-w[i,1]+1
    j = t[i,1]
    copyto!(x, w[i,1], x, j, N)
    if i > 1
      fill_s = w[i-1,2]+1
      fill_e = w[i,1]-1
      x[fill_s:fill_e] .= mx
    end
  end
  return nothing
end

function gapfill(x::Array{T,1}, t::Array{Int64,2}, fs::Float64; m::Bool=true) where T<: Real
  y = deepcopy(x)
  gapfill!(y, t, fs, m=m)
  return y
end

# replace NaNs with the mean
function nanfill!(x::Array{T,1}) where T<: Real
  J = findall(isnan.(x))
  if !isempty(J)
    if length(J) == length(x)
      fill!(x, zero(T))
    else
      x[J] .= T(mean(findall(isnan.(x).==false)))
    end
  end
  return nothing
end


text:
using SignedDistanceFunctions
using StaticArrays
using GeometryBasics
using Rotations
using CoordinateTransformations
using Test

@testset "SignedDistanceFunctions.jl" begin
    T = Float64
    spheresdf1 = SphereSignedDistanceFunction(0.5)
    spheresdf2 = SphereSignedDistanceFunction(1.0)
    unionsdf = UnionSignedDistanceFunction((spheresdf1, spheresdf2))
    intersectsdf = IntersectSignedDistanceFunction((spheresdf1, spheresdf2))
    negatesdf1 = NegatedSignedDistanceFunction(spheresdf1)
    negatesdf2 = NegatedSignedDistanceFunction(spheresdf2)

    p1 = zero(Point{3, T})
    p2 = Point{3, T}(1.3, 0., 0.)
    p3 = Point{3, T}(1.2, 0., 0.)
    @test spheresdf1(p1) ≈ -0.5
    @test negatesdf1(p1) ≈ 0.5
    @test spheresdf2(p1) ≈ -1.0
    @test negatesdf2(p1) ≈ 1.0
    @test unionsdf(p1) ≈ -1.0
    @test intersectsdf(p1) ≈ -0.5

    @test spheresdf1(p2) ≈ 0.8
    @test spheresdf2(p2) ≈ 0.3
    @test intersectsdf(p2) ≈ 0.8
    @test unionsdf(p2) ≈ 0.3
    @test normal(spheresdf1, p2) ≈ Vec{3, T}(1., 0., 0.)
    @test normal(negatesdf1, p2) ≈ Vec{3, T}(-1., 0., 0.)
    @test normal(spheresdf1, p3) ≈ Vec{3, T}(1., 0., 0.)
    @test normal(negatesdf1, p3) ≈ Vec{3, T}(-1., 0., 0.)

    bvsdf = BoundingVolumeSignedDistanceFunction(0.25, spheresdf2, spheresdf1)

    @test bvsdf(p1) ≈ -0.5
    @test bvsdf(p2) ≈ 0.3
    @test bvsdf(p3) ≈ 0.7

    # translation transformation tests
    translation1 = Translation(1., 0., 0.)
    translated_sphere1 = transform_sdf(translation1, spheresdf1)
    translated_sphere2 = transform_sdf(translation1, spheresdf2)
    @test translated_sphere1(p1) ≈ 0.5
    @test translated_sphere1(p2) ≈ -0.2
    @test translated_sphere2(p1) ≈ 0.0
    @test translated_sphere2(p2) ≈ -0.7

    translation2 = Translation(-0.5, 0., 0.)
    translated2x_sphere1 = transform_sdf(translation2, translated_sphere1)
    translated2x_sphere2 = transform_sdf(translation2, translated_sphere2)
    @test translated2x_sphere1(p1) ≈ 0.0
    @test translated2x_sphere1(p2) ≈ 0.3
    @test translated2x_sphere2(p1) ≈ -0.5
    @test translated2x_sphere2(p2) ≈ -0.2

    rotmat = LinearMap(RotZ(pi/2))
    affinemat = compose(rotmat, translation1)
    invaffinemat = inv(affinemat)

    p4 = Point{3, T}(1., 0., 0.)
    p5 = Point{3, T}(0., 1., 0.)
    p6 = Point{3, T}(0., -1., 0.)

    @test isapprox(invaffinemat(p5), Point{3, T}(0., 0., 0.); atol=1e-5)
    @test isapprox(invaffinemat(p1), Point{3, T}(-1., 0., 0.); atol=1e-5)
    @test isapprox(invaffinemat(p2), Point{3, T}(-1., -1.3, 0.); atol=1e-5)

    affine_sphere1 = transform_sdf(affinemat, spheresdf1)
    @test normal(affine_sphere1, p1) ≈ Vec{3, T}(0., 1., 0.)

end


text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

# using pkgs
using ForwardDiff
using Test
using TFT

# pkg constants
const _FD = ForwardDiff

# tolerances
atol = 1e-6

@testset "TFT.jl" begin

    include("prob/tft.jl")
    include("prob/ftft.jl")

end


text:
export blockBiCGSTB

function blockBiCGSTB{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,2}; kwargs...)
	TYPE = promote_type(T1,T2);
	Ax = zeros(TYPE,size(b));
	return blockBiCGSTB(x -> A_mul_B!(one(TYPE),A,x,zero(TYPE),Ax),b;kwargs...);
end

blockBiCGSTB(A,b; kwargs...) =  blockBiCGSTB(x -> A*x,b; kwargs...)

"""
x,flag,err,iter,resvec = blockBiCGSTB(A,b,tol=1e-6,maxIter=100,M1=identity,M2=identity,x=[],out=0)

Block version of BiConjugate Gradient Stabilized Method applied to the linear system Ax=b with multiple right hand sides.

Based on the paper:
El Guennouni, A., K. Jbilou, and H. Sadok. "A block version of BiCGSTAB for linear systems with multiple right-hand sides."
Electronic Transactions on Numerical Analysis 16.129-142 (2003): 2.

Input:

  A       - computing A*x
  b       - right hand side vector
  tol     - error tolerance
  maxIter - maximum number of iterations
  M1,M2   - preconditioners, either matrices or functions computing M1\\x or M2\\x
		  - If memory is reused in M1 and M2, make sure each has its own copy of return vector
  x       - starting guess
  out     - flag for output (-1: no output, 0 : only errors, 1 : final status, 2: relres at each iteration)

Output:

  x       - solution
  flag    - exit flag (  	0 : desired tolerance achieved,
  						-1 : maxIter reached without converging
  						-2 : ||rho||_F equals to zero
  						-3 : norm(s)/bnrm2 < tol
  						-4 : omega < 1e-16
  						-9 : right hand side was zero)
  err     - error, i.e., vecnorm(A*x-b)/vecnorm(b)
  iter    - number of iterations
  resvec  - error at each iteration
"""

function blockBiCGSTB{T}(A::Function, b::Array{T}; tol::Real=1e-6, maxIter::Int=100, M1=identity, M2=identity,x::Array=[],out::Int=0)

	n   = size(b,1);
	m   = size(b,2);

	N = m*n;

	if vecnorm(b)==0;
		return zeros(eltype(b),n,m),-9;
	end
	M1f =  isa(M1,Function) ? M1 : x -> M1\x
	M2f =  isa(M2,Function) ? M2 : x -> M2\x

	constOne = one(T);
	constZero = zero(T);

	if length(x)==0
		x = zeros(eltype(b),n,m)
		r = copy(b)
		r_tld = b; # We are not changing r_tld in this function, so we assign it to point to b that is also held constant.
	elseif vecnorm(x)==0.0
		r = copy(b)
		r_tld = b; # We are not changing r_tld in this function, so we assign it to point to b that is also held constant.
	else
		r = b - A(x);
		r_tld = copy(r); # We are not changing r_tld in this function, but we do change r;
	end
	v 	= zeros(T,n,m);
    p   = copy(r);

	resvec = zeros(maxIter+1)
	bnrm2 = vecnorm( b )

	resid   = vecnorm( r ) / bnrm2;
	resvec[1] = resid

	alpha = 1.0
	omega = 1.0

	iter = 1
	flag = -1
	rho1 = 0.0

	if out==2
		println("=== block bicgstb ===")
		println(@sprintf("%4s\t%7s","iter","relres"))
	end
	RtV = [];
	t = [];
	for iter = 1:maxIter
		rho = BLAS.gemm('C','N', constOne, r_tld, r); 				# equivalent to rho   = r_tld'*r;

		if ( vecnorm(rho) < 1e-13 ); flag = -2; break; end

		if ( iter > 1 )
			RtT = BLAS.gemm('C','N', constOne, r_tld, t); # RtT  = r_tld'*t;
			beta = - RtV \ RtT;
			# the next lines do: p = r + (p-v*omega)*beta; We run over v here...
			BLAS.scal!(N,-omega,v,1); 									# v = -v*omega
			BLAS.axpy!(N,constOne,p,1,v,1); 							# v += p;
			BLAS.gemm!('N','N', constOne, v, beta,constZero,p); 		# p = v*beta;
			BLAS.axpy!(N,constOne,r,1,p,1); 							# p = p + r
		end

		p_hat = M1f(p);      # compute M1\p
		p_hat = M2f(p_hat) ; # compute M2\phat

		t     = A(p_hat);    # compute A*phat

		v[:]  = t;

		RtV = BLAS.gemm('C','N', constOne, r_tld, v); 					# RtV = r_tld'*v;
		alpha = RtV\rho;
		BLAS.gemm!('N','N', constOne, p_hat, alpha,constOne,x); 		# x = x + p_hat*alpha # After this line p_hat is done.
		BLAS.gemm!('N','N', -constOne, v, alpha,constOne,r);  			# r = r - v*alpha;

		resid = vecnorm( r ) / bnrm2
		if ( resid < tol )
			iter -=1
			flag  = -3; break
		end

		s_hat = M1f(r)      # compute M1\r Note that it might be that M1f and M2f return the same memory vector.
		s_hat = M2f(s_hat)  # compute M2\shat

		t     = A(s_hat)    # compute A*shat

		omega = vecdot(t,r) / vecdot(t,t) ;

		BLAS.axpy!(N,omega,s_hat,1,x,1); 	# x = x + s_hat*omega
		BLAS.axpy!(N,-omega,t,1,r,1); 		# r = r - omega * t

		resid = vecnorm( r ) / bnrm2
		resvec[iter+1] = resid
		if out==2
			println(@sprintf("%3d\t%1.2e",iter,resvec[iter+1]))
		end
		if ( resid <= tol )
			flag = 0; break
		end
		if norm(omega) < 1e-16; flag = -4; break; end

	end
	if out>=0
		if flag==-1
			println(@sprintf("blockBiCGSTB iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
																								maxIter,resvec[iter],tol))
		elseif flag==-2
			println(@sprintf("blockBiCGSTB: rho equal to zero at iteration %d. Returned residual has norm %1.2e.", iter,resvec[iter+1]))
		elseif flag==-4
			println(@sprintf("blockBiCGSTB: omega < 1e-16"))
		elseif out>=1
			if flag == -3
				println(@sprintf("blockBiCGSTB achieved desired tolerance at iteration %d.5. Residual norm is %1.2e.",iter,resid))
			else
				println(@sprintf("blockBiCGSTB achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,resvec[iter+1]))
			end
		end
	end
	return x, flag,resvec[iter+1],iter,resvec[1:iter+1]
end


text:
################################################################################
# Solver Methods
################################################################################

function newton_solve!(prob::GameProblem{KN,n,m,T,SVd,SVx}) where {KN,n,m,T,SVd,SVx}
	model = prob.model
	core = prob.core
	game_con = prob.game_con
	opts = prob.opts

	# Set initial trajectory
	Random.seed!(opts.seed)
	init_traj!(prob.pdtraj; x0=prob.x0, f=opts.f_init, amplitude=opts.amplitude_init, s=opts.shift)
	init_traj!(prob.pdtraj_trial; x0=prob.x0, f=opts.f_init, amplitude=opts.amplitude_init, s=opts.shift)
	init_traj!(prob.Δpdtraj; x0=prob.x0, f=zeros, amplitude=0.0)

	rollout!(RK3, prob.model, prob.pdtraj.pr)
	rollout!(RK3, prob.model, prob.pdtraj_trial.pr)
	# Set the initial penalties
	prob.pen.ρ = SVector{1,T}([opts.ρ_0])
	prob.pen.ρ_trial = SVector{1,T}([opts.ρ_trial])

	# Reset Statistics and constraints
	reset!(prob.stats)
	opts.dual_reset ? reset!(game_con) : nothing
	# Iterative solve
	out = 0
	t_elap = 0.0
	Δ = 0.0
    for k = 1:opts.outer_iter
		out = k
		# plot!(model, prob.pdtraj.pr)
		# Initialize regularization and failed line search count.
		set!(opts.reg, opts.reg_0)
		LS_count = 0
		opts.inner_print ? display_solver_header() : nothing
		#XXX opts.inner_print ? display_condition_header() : nothing
        for l = 1:opts.inner_iter
			set!(opts.reg, opts.reg_0*l^4)
			t_elap = @elapsed begin
				LS_count, control_flow, Δ = inner_iteration(prob, LS_count, t_elap, Δ, k, l)
			end
			LS_count >= 1 || control_flow == :break ? break : nothing
		end

		# Visualize trajectory
		#XXX live_vis && (opts.outer_iter-k)%1==0 ? visualize!(vis, model, pdtraj.q) : nothing
		# Kick out before updating the penalty and duals
		if k == opts.outer_iter || (
			prob.stats.dyn_vio[end].max < opts.ϵ_dyn &&
			prob.stats.con_vio[end].max < opts.ϵ_con &&
			prob.stats.sta_vio[end].max < opts.ϵ_sta &&
			prob.stats.opt_vio[end].max < opts.ϵ_opt)
			break
		end
		# Dual Ascent
		evaluate!(game_con, prob.pdtraj.pr)
		dual_update!(game_con)
		# Increasing Schedule
		prob.pen.ρ = min.(prob.pen.ρ * opts.ρ_increase, opts.ρ_max)
		penalty_update!(game_con)
    end
	record!(prob.stats, prob, model, game_con, prob.pdtraj, t_elap, Δ, out)
    return nothing
end

function inner_iteration(prob::GameProblem, LS_count::Int, t_elap::T, Δ::T, k::Int, l::Int) where {T}
	core = prob.core
	opts = prob.opts
	# plot!(prob.model, prob.pdtraj.pr)

	# Residual
	residual!(prob, prob.pdtraj)
	opts.regularize ? regularize_residual!(core, opts, prob.pdtraj, prob.pdtraj) : nothing # should do nothing since we regularize around pdtraj
	record!(prob.stats, prob, prob.model, prob.game_con, prob.pdtraj, t_elap, Δ, k)
	res_norm = norm(core.res, 1)/length(core.res)

	# Reset Δ_traj to zero
	Δ = 0.0
	if prob.stats.opt_vio[end].max < opts.ϵ_opt
		return LS_count, :break, Δ
	end
	# Residual Jacobian
	residual_jacobian!(prob)
	regularize_residual_jacobian!(prob)

	Δtraj = - \(lu(core.jac), core.res)
	set_traj!(core, prob.Δpdtraj, Δtraj)

	# Line Search
	α, j = line_search(prob, res_norm)
	failed_ls = j == opts.ls_iter
	failed_ls ? LS_count += 1 : LS_count = 0
	update_traj!(prob.pdtraj, prob.pdtraj, α, prob.Δpdtraj)
	Δ = Δ_step(prob.Δpdtraj, α)
	if Δ < opts.Δ_min
		return LS_count, :break, Δ
	end

	opts.inner_print ? display_solver_data(k, l, j, Δ, res_norm, opts.reg) : nothing
	#XXX opts.inner_print ? display_condition_data(res) : nothing
	return LS_count, :continue, Δ
end

function line_search(prob::GameProblem, res_norm::T) where {T}
	core = prob.core
	opts = prob.opts

	j = 1
	α = 1.0
	while j < opts.ls_iter
		update_traj!(prob.pdtraj_trial, prob.pdtraj, α, prob.Δpdtraj)
		residual!(prob, prob.pdtraj_trial)
		opts.regularize ? regularize_residual!(core, opts, prob.pdtraj_trial, prob.pdtraj) : nothing# should add regularization term to prevent deviation from pdtraj
		res_norm_trial = norm(core.res, 1)/length(core.res)
		if res_norm_trial <= (1.0-α*opts.β)*res_norm
			LS_count = 0
			break
		else
			α *= opts.α_decrease # step size decrease
			j += 1
		end
	end
	return α, j
end

################################################################################
# Iterative Best Response Solver Methods
################################################################################

function ibr_newton_solve!(prob::GameProblem{KN,n,m,T,SVd,SVx};
		ibr_opts::IBROptions=IBROptions()) where {KN,n,m,T,SVd,SVx}
	opts = prob.opts
	p = prob.probsize.p

	# Reset Statistics
	reset!(prob.stats)
	# Set initial trajectory
	Random.seed!(opts.seed)
	init_traj!(prob.pdtraj; x0=prob.x0, f=opts.f_init, amplitude=opts.amplitude_init, s=opts.shift)
	init_traj!(prob.pdtraj_trial; x0=prob.x0, f=opts.f_init, amplitude=opts.amplitude_init, s=opts.shift)
	init_traj!(prob.Δpdtraj; x0=prob.x0, f=zeros, amplitude=0.0)

	rollout!(RK3, prob.model, prob.pdtraj.pr)
	rollout!(RK3, prob.model, prob.pdtraj_trial.pr)

	# Records all the players that changed their strategies at the last iteration:
	# when all the players stuck to their strategies we can safely exit the IBR iteration.
	Δ_change = trues(p)
	for q = 1:ibr_opts.ibr_iter # need a better stopping criterion, based on convergence
		for id = 1:p
			i = ibr_opts.ordering[id]
			ibr_newton_solve!(prob, i)
			Δ_change[i] = !(ibr_opts.Δ_min > maximum(prob.stats.Δ_traj))
			ibr_opts.live_plotting && plot!(prob.model, prob.pdtraj.pr)
			ibr_opts.live_plotting && plot!(prob_ibr.stats)
		end
		residual!(prob, prob.pdtraj)
		res = norm(prob.core.res, 1)/length(prob.core.res)
		# If we see no change for one full IBR iteration, then we break out of the loop.
		all(.!Δ_change) && break
	end
	return nothing
end

function ibr_newton_solve!(prob::GameProblem{KN,n,m,T,SVd,SVx}, i::Int) where {KN,n,m,T,SVd,SVx}
	model = prob.model
	core = prob.core
	game_con = prob.game_con
	opts = prob.opts

	# Set the initial penalties
	prob.pen.ρ = SVector{1,T}([opts.ρ_0])
	prob.pen.ρ_trial = SVector{1,T}([opts.ρ_trial])

	# Reset Statistics and constraints
	# reset!(prob.stats)
	if opts.dual_reset
		reset!(game_con)
		reset_duals!(prob.pdtraj)
		reset_duals!(prob.pdtraj_trial)
	end

	# Iterative solve
	out = 0
	t_elap = 0.0
	Δ = 0.0
    for k = 1:opts.outer_iter
		out = k
		# plot!(model, prob.pdtraj.pr)
		# Initialize regularization and failed line search count.
		set!(opts.reg, opts.reg_0)
		LS_count = 0
		opts.inner_print ? display_solver_header() : nothing
		#XXX opts.inner_print ? display_condition_header() : nothing
        for l = 1:opts.inner_iter
			set!(opts.reg, opts.reg_0*l^4)
			t_elap = @elapsed begin
				LS_count, control_flow, Δ = ibr_inner_iteration(prob, LS_count, t_elap, Δ, k, l, i)
			end
			LS_count >= 1 || control_flow == :break ? break : nothing
		end

		# Visualize trajectory
		#XXX live_vis && (opts.outer_iter-k)%1==0 ? visualize!(vis, model, pdtraj.q) : nothing
		# Kick out before updating the penalty and duals
		if k == opts.outer_iter || (
			prob.stats.dyn_vio[end].max < opts.ϵ_dyn &&
			prob.stats.con_vio[end].max < opts.ϵ_con &&
			prob.stats.sta_vio[end].max < opts.ϵ_sta &&
			prob.stats.opt_vio[end].max < opts.ϵ_opt)
			break
		end
		# Dual Ascent
		evaluate!(game_con, prob.pdtraj.pr)
		dual_update!(game_con)
		# Increasing Schedule
		prob.pen.ρ = min.(prob.pen.ρ * opts.ρ_increase, opts.ρ_max)
		penalty_update!(game_con)
    end
	record!(prob.stats, prob, model, game_con, prob.pdtraj, t_elap, Δ, out, i)
    return nothing
end

function ibr_inner_iteration(prob::GameProblem, LS_count::Int, t_elap, Δ::T, k::Int, l::Int, i::Int) where {T}
	core = prob.core
	opts = prob.opts
	verti_mask = vertical_mask(prob.core, i)
	horiz_mask = horizontal_mask(prob.core, i)
	# plot!(prob.model, prob.pdtraj.pr)

	# Residual
	ibr_residual!(prob, prob.pdtraj, i)
	opts.regularize ? regularize_ibr_residual!(core, opts, prob.pdtraj, prob.pdtraj, i) : nothing # should do nothing since we regularize around pdtraj
	record!(prob.stats, prob, prob.model, prob.game_con, prob.pdtraj, t_elap, Δ, k, i)
	res_norm = norm(core.res[verti_mask], 1)/length(core.res[verti_mask])

	# Reset Δ_traj to zero
	Δ = 0.0
	if prob.stats.opt_vio[end].max < opts.ϵ_opt
		return LS_count, :break, Δ
	end
	# Residual Jacobian
	ibr_residual_jacobian!(prob, i)
	regularize_ibr_residual_jacobian!(prob, i)

	Δtraj = zeros(prob.probsize.S)
	Δtraj[horiz_mask] = - \(lu(core.jac[verti_mask, horiz_mask]), core.res[verti_mask])
	set_traj!(core, prob.Δpdtraj, Δtraj)

	# Line Search
	α, j = ibr_line_search(prob, res_norm, i)
	failed_ls = j == opts.ls_iter
	failed_ls ? LS_count += 1 : LS_count = 0
	update_traj!(prob.pdtraj, prob.pdtraj, α, prob.Δpdtraj)
	Δ = Δ_step(prob.Δpdtraj, α)
	if Δ < opts.Δ_min
		return LS_count, :break, Δ
	end

	opts.inner_print ? display_solver_data(k, l, j, Δ, res_norm, opts.reg) : nothing
	#XXX opts.inner_print ? display_condition_data(res) : nothing
	return LS_count, :continue, Δ
end

function ibr_line_search(prob::GameProblem, res_norm::T, i::Int) where {T}
	core = prob.core
	opts = prob.opts
	verti_mask = vertical_mask(prob.core, i)

	j = 1
	α = 1.0
	while j < opts.ls_iter
		update_traj!(prob.pdtraj_trial, prob.pdtraj, α, prob.Δpdtraj)
		ibr_residual!(prob, prob.pdtraj_trial, i)
		opts.regularize ? regularize_ibr_residual!(core, opts, prob.pdtraj_trial, prob.pdtraj, i) : nothing# should add regularization term to prevent deviation from pdtraj
		res_norm_trial = norm(core.res[verti_mask], 1)/length(core.res[verti_mask])
		if res_norm_trial <= (1.0-α*opts.β)*res_norm
			LS_count = 0
			break
		else
			α *= opts.α_decrease # step size decrease
			j += 1
		end
	end
	return α, j
end


text:
using Test
using Newick
using TranscodingStreams

using Base.CoreLogging:Debug

@testset "Newick" begin

@testset "no nodes are named" begin

    #         12345678
    newick = "(,,(,));"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "nodefinish 2"),
        (:debug, "nodefinish 3"),
        (:debug, "cladestart 4"),
        (:debug, "nodefinish 5"),
        (:debug, "nodefinish 6"),
        (:debug, "cladefinish 6"),
        (:debug, "nodefinish 7"),
        (:debug, "cladefinish 7"),
        (:debug, "finish 8"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "no nodes are named" newick root

    Newick.print_tree(root)

end

@testset "leaf nodes are named" begin

    #         000000000111
    #         123456789012
    newick = "(A,B,(C,D));"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "name 3"),
        (:debug, "nodefinish 3"),
        (:debug, "name 5"),
        (:debug, "nodefinish 5"),
        (:debug, "cladestart 6"),
        (:debug, "name 8"),
        (:debug, "nodefinish 8"),
        (:debug, "name 10"),
        (:debug, "nodefinish 10"),
        (:debug, "cladefinish 10"),
        (:debug, "nodefinish 11"),
        (:debug, "cladefinish 11"),
        (:debug, "finish 12"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "leaf nodes are named" newick root

    Newick.print_tree(root)

end

@testset "all nodes are named" begin
    #         00000000011111
    #         12345678901234
    newick = "(A,B,(C,D)E)F;"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "name 3"), #A
        (:debug, "nodefinish 3"),
        (:debug, "name 5"), #B
        (:debug, "nodefinish 5"),
        (:debug, "cladestart 6"),
        (:debug, "name 8"), #C
        (:debug, "nodefinish 8"),
        (:debug, "name 10"), #D
        (:debug, "nodefinish 10"),
        (:debug, "cladefinish 10"),
        (:debug, "name 12"), #E
        (:debug, "nodefinish 12"),
        (:debug, "cladefinish 12"),
        (:debug, "name 14"), #F
        (:debug, "finish 14"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "all nodes are named" newick root

    Newick.print_tree(root)

end

@testset "all but root node have a distance to parent" begin

    #         0000000001111111111222222222
    #         1234567890123456789012345678
    newick = "(:0.1,:0.2,(:0.3,:0.4):0.5);"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "distance 6"),
        (:debug, "nodefinish 6"),
        (:debug, "distance 11"),
        (:debug, "nodefinish 11"),
        (:debug, "cladestart 12"),
        (:debug, "distance 17"),
        (:debug, "nodefinish 17"),
        (:debug, "distance 22"),
        (:debug, "nodefinish 22"),
        (:debug, "cladefinish 22"),
        (:debug, "distance 27"),
        (:debug, "nodefinish 27"),
        (:debug, "cladefinish 27"),
        (:debug, "finish 28"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "all but root node have a distance to parent" newick root

    Newick.print_tree(root)

end

@testset "all have a distance to parent" begin

    #         00000000011111111112222222222333
    #         12345678901234567890123456789012
    newick = "(:0.1,:0.2,(:0.3,:0.4):0.5):0.0;"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "distance 6"),
        (:debug, "nodefinish 6"),
        (:debug, "distance 11"),
        (:debug, "nodefinish 11"),
        (:debug, "cladestart 12"),
        (:debug, "distance 17"),
        (:debug, "nodefinish 17"),
        (:debug, "distance 22"),
        (:debug, "nodefinish 22"),
        (:debug, "cladefinish 22"),
        (:debug, "distance 27"),
        (:debug, "nodefinish 27"),
        (:debug, "cladefinish 27"),
        (:debug, "distance 32"),
        # (:debug, "nodefinish 32"), #TODO: not sure about this one. But, it does yeild the correct result.
        (:debug, "finish 32"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "all have a distance to parent" newick root

    Newick.print_tree(root)

end

@testset "distances and leaf names (popular)" begin

    #         00000000011111111112222222222333
    #         12345678901234567890123456789012
    newick = "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "name 3"),
        (:debug, "distance 7"),
        (:debug, "nodefinish 7"),
        (:debug, "name 9"),
        (:debug, "distance 13"),
        (:debug, "nodefinish 13"),
        (:debug, "cladestart 14"),
        (:debug, "name 16"),
        (:debug, "distance 20"),
        (:debug, "nodefinish 20"),
        (:debug, "name 22"),
        (:debug, "distance 26"),
        (:debug, "nodefinish 26"),
        (:debug, "cladefinish 26"),
        (:debug, "distance 31"),
        (:debug, "nodefinish 31"),
        (:debug, "cladefinish 31"),
        (:debug, "finish 32"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "distances and leaf names (popular)" newick root

    Newick.print_tree(root)

end

@testset "distances and all names" begin

    #         0000000001111111111222222222233333
    #         1234567890123456789012345678901234
    newick = "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "name 3"),
        (:debug, "distance 7"),
        (:debug, "nodefinish 7"),
        (:debug, "name 9"),
        (:debug, "distance 13"),
        (:debug, "nodefinish 13"),
        (:debug, "cladestart 14"),
        (:debug, "name 16"),
        (:debug, "distance 20"),
        (:debug, "nodefinish 20"),
        (:debug, "name 22"),
        (:debug, "distance 26"),
        (:debug, "nodefinish 26"),
        (:debug, "cladefinish 26"),
        (:debug, "name 28"),
        (:debug, "distance 32"),
        (:debug, "nodefinish 32"),
        (:debug, "cladefinish 32"),
        (:debug, "name 34"),
        (:debug, "finish 34"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "distances and all names" newick root

    Newick.print_tree(root)

end

@testset "a tree rooted on a leaf node (rare)" begin

    #         00000000011111111112222222222333333
    #         12345678901234567890123456789012345
    newick = "((B:0.2,(C:0.3,D:0.4)E:0.5)A:0.1)F;"

    stream = TranscodingStreams.NoopStream(IOBuffer(newick))

    root = Newick.Record()

    @test_logs(
        (:debug, "cladestart 1"),
        (:debug, "cladestart 2"),
        (:debug, "name 4"),
        (:debug, "distance 8"),
        (:debug, "nodefinish 8"),
        (:debug, "cladestart 9"),
        (:debug, "name 11"),
        (:debug, "distance 15"),
        (:debug, "nodefinish 15"),
        (:debug, "name 17"),
        (:debug, "distance 21"),
        (:debug, "nodefinish 21"),
        (:debug, "cladefinish 21"),
        (:debug, "name 23"),
        (:debug, "distance 27"),
        (:debug, "nodefinish 27"),
        (:debug, "cladefinish 27"),
        (:debug, "name 29"),
        (:debug, "distance 33"),
        (:debug, "nodefinish 33"),
        (:debug, "cladefinish 33"),
        (:debug, "name 35"),
        (:debug, "finish 35"),
        min_level=Debug,
        Newick.read!(stream, root)
    )

    @info "a tree rooted on a leaf node (rare)" newick root

    Newick.print_tree(root)

end

end # testset "Newick"


text:
module Unpack

export
    unpack,
    CompileOptions,
    CompiledField,
    @compile

struct CompiledField
    field::Symbol
    # aliases::Vector{Symbol}
    # remap::Tuple{Bool, Symbol}
    # case::Symbol
    # default::Tuple{Bool, Any}
    # skip::Bool
    # key_resolver::Function
    # callback::Function
end

struct CompileOptions
    parent::Symbol
    fields::Vector{CompiledField}
end

unpacker_name(dt) = Symbol(:_UNPACK_, dt)
unpacker_type(dt) = eval(unpacker_name(dt))

unpack() = nothing
unpack(dt::DataType, ustore) = unpack(unpacker_type(dt)(), ustore)

function field_processor(field::CompiledField)
    fnm = QuoteNode(:($(field.field)))
    e = :(us -> us[])
    append!(e.args[2].args[2].args, [fnm])
    e
end

function struct_processor(upk::Symbol, opts::CompileOptions)::Expr
    fn = :(unpack(v::$upk, ustore) = $(opts.parent)())

    for opt in opts.fields
        λ = field_processor(opt)
        append!(fn.args[2].args[2].args, [:($λ(ustore))])
    end

    fn
end

macro compile(ex)
    struct_name = unpacker_name(ex.args[1])
    constructor = :(struct $struct_name end)
    eval(constructor)
    e = struct_processor(struct_name, eval(ex.args[2]))
    eval(e)
end

end # module


text:
using Flux, FluxTools

# Explicit model definition
m0 = Model((Dense(3, 5), Dense(5, 4), Dense(3, 1))) do x, (l1, l2, l3)
    x = l1(x)
    x = l2(x)
    return l3(x[1:3, :])
end

m0(rand(3, 12))

@show m0.layers

# Model definition via macro. As many layers as there are `$` signs.
m1 = @model x begin
    y = $(Dense(2, 3))(x)
    z = y .* 1.2
    return $(Dense(3, 4))(abs.(z))
end

m1(rand(2, 10))

@show m1.layers

# Nested example

m2 = @model x begin
    y = $m1(x)
    return softmax($(Dense(4, 2))(y))
end

m2(rand(2, 10))

@show length(params(m2)) == length(params(m1)) + length(params(Dense(4, 2)))


text:
###########################
#Main function that iterates through sampling new parameters
###########################
function ptMCMC(prob::DiffEqBase.DEProblem,alg,priors,parBounds,lossFunc,iterations = 1000)
#This is a Metropolis-Hastings MCMC algorithm designed to optimize ODE parameters to fit a given dataset.

    ###########################
    #Basic setup for the run
    ###########################
    parNum = length(priors) #Number of parameters
    chainNum = 3 #number of chains
    βmin, βmax = [0.1, 1.0] #Bounds for chain temperatures
    β = chainNum==1 ? 1 : range(βmin, βmax, length=chainNum) #assign each chain a temperature

    #Initialize objects to store chains and probabilities
    parCurr = [rand.(priors) for i = 1:chainNum] #Starting parameters
    chain = zeros(iterations,parNum,chainNum) #Initialize chain
    logProb = zeros(iterations,chainNum)  #log-posterior chain

    #Storage for when the proposals are accepted
    acceptRate = zeros(iterations,chainNum) #Acceptance ratio
    acceptNum = zeros(Int64,chainNum) #Number of acceptances

    #Generate the log likelihood function (need to include std dev of data)
    SumSqErr = build_loss_objective(prob,alg,lossFunc)
    LogLikeLossFunc(pars) = -0.5*SumSqErr(pars)
    logProbCurr = LogLikeLossFunc.(parCurr)

    #Create proposal distributions for each parameter to general new guesses
    σ = 0.7 # "Tuning parameter", i.e. how big are your jumps

    #Swap Information
    swapRatio = zeros(iterations) #Acceptance ratio of swaps
    swapNum = 0 #Number of swap acceptances

    #Input the first sample for the chain and log probability
    chain[1,:,:] = hcat(parCurr...) #hcat: array of vectors → 2D array
    logProb[1,:] = logProbCurr

    ###########################
    #Main Loop: Go through number of samples desired
    ###########################
    @showprogress 1 "Computing MCMC Samples..." for i = 2:iterations
        ###########################
        #Apply MH_MCMC step to each chain
            #(similar to lapply in R or map in python)
        ###########################
        #for pmap I'll have to be careful about rand...
        chainResult = MH_MCMC.(parCurr,logProbCurr,acceptNum,β,σ,[parBounds],LogLikeLossFunc)

        #Store the results from the MCMC step
        for (chnIndex, chn) in enumerate(chainResult)
            #Update the current parameters
            parCurr[chnIndex] = chn[1]
            logProbCurr[chnIndex] = chn[2]
            acceptNum[chnIndex] += chn[3]

            #Update the chains storage containers
            chain[i,:,chnIndex] = parCurr[chnIndex]
            logProb[i,chnIndex] = logProbCurr[chnIndex]
            acceptRate[i,chnIndex] = acceptNum[chnIndex] / (i-1)
        end

        ###########################
        #Let's swap some chains (going from highest to lowest temp chain)
        ###########################
    if chainNum > 1
        for chnIndex = chainNum:-1:2
            #Calculate the error and temperature change between chains
            ΔE = logProbCurr[chnIndex] - logProbCurr[chnIndex-1]
            Δβ = β[chnIndex] - β[chnIndex-1]

            #Check if the swap is g0od enough
            if ΔE*Δβ > rand(LogNormal(0,1))
                #swapIndex = chnIndex:-1:chnIndex-1
                swapIndex = reverse(chnIndex-1:chnIndex) #idk if more readable
                #Swap the parameters
                parCurr[reverse(swapIndex)] = parCurr[swapIndex]
                #Swap the error/energy
                logProbCurr[reverse(swapIndex)] = logProbCurr[swapIndex]

                #Count the swap
                swapNum += 1
            end
        end
        #Save swap percent
        swapRatio[i] = swapNum / (i-1)
    end

end #iterations

#return the final results of the chains
return chain, logProb, acceptRate

end #parSample2

###########################
#MH algorithm for one chain
###########################
function MH_MCMC(parCurr,logProbCurr,acceptNum,β,σ,parBounds,LogLikeLossFunc)
    #Propose a new parameter set
    parPropose = NewParameter!(parCurr,parBounds,σ)
    #Calculate change in error
    logProbPropose = LogLikeLossFunc(parPropose)
    ΔE = logProbPropose - logProbCurr
    #Calculate the acceptance ratio correction
    ΔC = cdfCalculator(parPropose,parBounds,σ) - cdfCalculator(parCurr,parBounds,σ)

    #Calculate the acceptance ratio
    acceptRatio = min(1,exp(β*ΔE + ΔC))

    if rand() < acceptRatio #If the acceptance passes...
        accept = 1
        #Update the current parameters with the proposed parameters
        return parPropose, logProbPropose, accept
    else
        accept = 0
        #Rejected, return the same values inputted
        return parCurr, logProbCurr, accept
    end

end #MH_MCMC

###########################
#Given the current parameter set, generate a new parameter set
###########################
function NewParameter!(parameters,parBounds,σ=0.3)
    #Add random normal number to current Parameter set
    #parameters += rand(Normal(0,σ^2), length(parameters))
    parameters = [rand(TruncatedNormal(p,σ^2,lb,ub)) for (p,(lb,ub)) in zip(parameters,parBounds)]
end

###########################
#Given a parameter set, calculate the CDF
#see: https://darrenjw.wordpress.com/tag/truncate/
###########################
function cdfCalculator(pars,bounds,σ=0.3)
    p = logcdf.([TruncatedNormal(0,σ^2,lb,ub) for (lb,ub) in bounds], pars)
    return sum(p)
end

max_stars_repo_path: Julia Monte Carlo Fitting/main.jl

max_issues_repo_path: Julia Monte Carlo Fitting/main.jl

max_forks_repo_path: Julia Monte Carlo Fitting/main.jl


text:
"""
$(SIGNATURES)

Convenience function to assemble the html of a page out of its parts:
- head
- content and pgfoot (which will be wrapped to form the body)
- foot.
"""
function build_page(head, content, pgfoot, foot)
    # (legacy support) if div_content is offered explicitly, it takes
    # precedence, otherwise use defaults
    dc = globvar("div_content")
    if isempty(dc)
        content_tag   = globvar("content_tag")
        content_class = globvar("content_class")
        content_id    = globvar("content_id")
    else
        content_tag   = "div"
        content_class = dc
        content_id    = ""
    end
    body = content * pgfoot
    # wrap the body in appropriate tags
    # <$tag class=$cclass>$content</$tag>
    # if tag is empty, do not wrap it (can be useful for some templates)
    if !isempty(content_tag)
        body = html_content(content_tag, body;
                            class=content_class, id=content_id)
    end
    return head * body * foot
end

"""
$(SIGNATURES)

Take a generated HTML string and apply prerendering and link fixing as
appropriate.
"""
function postprocess_page(pg)
    # Prerender if required (using JS tools)
    # order in which things are done matter a bit.
    if FD_ENV[:PRERENDER]
        # Maths (KATEX)
        if locvar(:hasmath)::Bool
            pg = js_prerender_katex(pg)
        end
        # Code (HIGHLIGHT.JS)
        if locvar(:hascode)::Bool && FD_CAN_HIGHLIGHT
            pg = js_prerender_highlight(pg)
            # remove script
            pg = replace(pg, r"<script.*?(?:highlight\.pack\.js|initHighlightingOnLoad).*?<\/script>"=>"")
        end
        if locvar(:hasmath)::Bool
            # remove katex scripts
            pg = replace(pg, r"<script.*?(?:katex\.min\.js|auto-render\.min\.js|renderMathInElement).*?<\/script>" => "")
        end
    end
    # append pre-path to links if required (see optimize)
    if FD_ENV[:FINAL_PASS]
        pg = fix_links(pg)
    end
    return pg
end

const PAGINATED = Set{String}()

"""
$(SIGNATURES)

Write a html page at the appropriate location and with the appropriate
structure. This is usually called specifying the scaffolding but can be done
without in which case the scaffolding is read from `layout`.
"""
function write_page(output_path::AS, content::AS;
                    head::T=nothing, pgfoot::T=nothing, foot::T=nothing,
                    )::String where T <: Union{Nothing,AS}
    slug = locvar(:slug)
    if !isempty(slug)
        output_path = form_custom_output_path(slug)
    end
    # NOTE
    #   - output_path is assumed to exist // see form_output_path
    #   - head/pgfoot/foot === nothing --> read (see franklin.jl)
    layout = path(:layout)
    if isnothing(head)
        head = read(joinpath(layout, "head.html"), String)
    end
    if isnothing(pgfoot)
        pgfoot = read(joinpath(layout, "page_foot.html"), String)
    end
    if isnothing(foot)
        foot = read(joinpath(layout, "foot.html"), String)
    end
    # convert the pieces
    head, ctt, pgfoot, foot = map(convert_html, (head, content, pgfoot, foot))

    # remove dirs from past pagination attempts (so we limit risk of stale
    # pagination folder with stale files)
    outdir = dirname(output_path)
    ispaginated = outdir in PAGINATED
    if ispaginated
        for elem in readdir(outdir)
            if all(isnumeric, elem) && first(elem) != '0'
                dpath = joinpath(outdir, elem)
                isdir(dpath) && rm(dpath, recursive=true)
            end
        end
    end

    # the previous convert call possibly resolved a {{paginate}} which will
    # have stored a :paginate_itr var, so we must branch on that
    if locvar(:paginate_itr) !== nothing
        union!(PAGINATED, (outdir,))
        name    = locvar(:paginate_itr)
        iter    = locvar(name)
        npp     = locvar(:paginate_npp)
        niter   = length(iter)
        n_pages = ceil(Int, niter / npp)
        for pgi = 1:n_pages
            # form the content multiple times
            sta_i = (pgi - 1) * npp + 1
            end_i = min(sta_i + npp - 1, niter)
            rge_i = sta_i:end_i
            ins_i = prod(String(e) for e in iter[rge_i])
            ctt_i = replace(ctt, PAGINATE => ins_i)
            # assemble, optimize and write
            pg = build_page(head, ctt_i, pgfoot, foot)
            pg = postprocess_page(pg)
            pgi == 1 && write(output_path, pg)
            dst = mkpath(joinpath(outdir, "$pgi"))
            write(joinpath(dst, "index.html"), pg)
        end
    else
        # maybe it was paginated but isn't anymore
        ispaginated && setdiff(PAGINATED, outdir)
        # convert any `{{...}}` that may be left and form the full page string
        pg = build_page(head, ctt, pgfoot, foot)
        pg = postprocess_page(pg)
        # write the html file where appropriate
        write(output_path, pg)
    end
    return pg
end

"""
$(SIGNATURES)

Take a path to an input markdown file (via `root` and `file`), then construct
the appropriate HTML page (inserting `head`, `pgfoot` and `foot`) and finally
write it at the appropriate place.
"""
function convert_and_write(root::String, file::String, head::String,
                           pgfoot::String, foot::String, output_path::String
                           )::Nothing
    # 1. read the markdown into string, convert it and extract definitions
    # 2. eval the definitions and update the variable dictionary, also retrieve
    # document variables (time of creation, time of last modif) and add those
    # to the dictionary.
    fpath = joinpath(root, file)
    # The curpath is the relative path starting after /src/ so for instance:
    # f1/blah/page1.md or index.md etc... this is useful in the code evaluation
    # and management of paths
    set_cur_rpath(fpath)
    # conversion
    content = convert_md(read(fpath, String))

    # adding document variables to the dictionary
    # note that some won't change and so it's not necessary to do this every
    # time but it takes negligible time to do this so ¯\_(ツ)_/¯
    # (and it's less annoying than keeping tabs on which file has
    # already been treated etc).
    s = stat(fpath)
    set_var!(LOCAL_VARS, "fd_ctime", fd_date(unix2datetime(s.ctime)))
    mtime = unix2datetime(s.mtime)
    set_var!(LOCAL_VARS, "fd_mtime_raw", Date(mtime))
    set_var!(LOCAL_VARS, "fd_mtime", fd_date(mtime))

    # Check if should add item
    #   should we generate ? otherwise no
    #   are we in the full pass ? otherwise no
    #   is there a `rss` or `rss_description` ? otherwise no
    cond_add = globvar(:generate_rss) &&   # should we generate?
                    FD_ENV[:FULL_PASS] &&  # are we in the full pass?
                    !all(e -> isempty(locvar(e)), ("rss", "rss_description"))
    # otherwise yes
    cond_add && add_rss_item()

    # Same for the sitemap
    cond_add = globvar(:generate_sitemap) && FD_ENV[:FULL_PASS]
    cond_add && add_sitemap_item()

    # Same for disallow robots locally
    cond_add = locvar(:robots_disallow_this_page)
    cond_add && add_disallow_item()

    pg = write_page(output_path, content; head=head, pgfoot=pgfoot, foot=foot)

    # 6. possible post-processing via the "on-write" function.
    FD_ENV[:ON_WRITE](pg, LOCAL_VARS)
    return nothing
end


text:
# Testing part 1
using Test, ReferenceTests

# dir to the part 1 script
include("../scripts-part1/Diffusion_3D.jl")

# Unit tests
Hτ = rand(5, 5, 5);
D_dx, D_dy, D_dz = 1.0, 1.0, 1.0;
@testset "Test heat flux macros" begin
    @test @qx(2,1,3) ≈ -(Hτ[3,2,4] - Hτ[2,2,4])
    @test @qx(4,3,1) ≈ -(Hτ[5,4,2] - Hτ[4,4,2])
    @test @qy(1,1,1) ≈ -(Hτ[2,2,2] - Hτ[2,1,2])
    @test @qy(3,2,3) ≈ -(Hτ[4,3,4] - Hτ[4,2,4])
    @test @qz(1,4,2) ≈ -(Hτ[2,5,3] - Hτ[2,5,2])
    @test @qz(2,3,4) ≈ -(Hτ[3,4,5] - Hτ[3,4,4])
end

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2; atol = 1e-5) for (v1,v2) in zip(values(d1), values(d2))])
inds = Int.(ceil.(LinRange(1, length(xc_g), 12)))
d = Dict(:X=> xc_g[inds], :H=>H_g[inds, inds, 15])

@testset "Reference Test" begin
    @test_reference "reftest-files/test_part_1.bson" d by=comp
end


text:
include("problems.jl")

function testProblem(::Type{so}, ::Type{l}, p::NamedTuple, ::Type{T} = Float64)::Bool where {so <: AbstractSolver, l <: AbstractLP, T}
    (sol, vsol,  A, b, c, isoptimal, isbounded, isfeasible) = (p.sol, p.vsol,  p.A, p.b, p.c, p.isoptimal, p.isbounded, p.isfeasible)
    linearform = l{T}(A, b, c)
    changeSolver!(so)
    solve!(linearform)
    goodstatus = false
    if isoptimal
        goodstatus = (linearform.status == Optimal)
    elseif !isbounded
        goodstatus = (linearform.status == Unbounded)
        return true
    elseif !isfeasible
        goodstatus = (linearform.status == Infeasible)
        return true
    end

    goodsol = isapprox(minimum([norm(xstar(linearform) - s) for s in sol]), 0, atol = 1e-6)
    goodsolvalue = isapprox(vstar(linearform), vsol, atol = 1e-6)
    #@show goodsol
    #@show goodsolvalue
    #@show goodstatus
    #@show linearform.status
    return goodsol && goodsolvalue && goodstatus
end

solverlist = [BaseEnumeration, StandardSimplexe]
representationlist = [LinearProblem]

for so in solverlist
    for rep in representationlist
        for p in ProblemTest.PL
            @test testProblem(so, rep, p)
        end
    end
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
using CSV
using DataFrames
using PGFPlots
using Statistics

"""
Plotting signals of 2 algorithms for visual comparison (similar to Fig. A.5)

NOTE:
-----
A lot of manual setting of things... signalnames, algonames, legendentry,
scalesurprise, ymax, xmax, etc...
"""

function compare_MB_signals(;
            results_fit_path =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/fits"),
            signalname1 = "SPE",
            # signalname1 = "gamma_surprise",
            # signalname1 = "Surprise",
            # signalname1 = "RPE",
            # signalname2 = "RPE",
            signalname2 = "Surprise",
            # signalname2 = "DeltaPolicy",
            # signalname2 = "SPE",
            # algoname1 = "forward_learner",
            algoname1 = "hybrid_actorcritic_FWD",
            # algoname1 = "actor_critic_TD_modulatedActor_1degree",
            # algoname1 = "sarsa_lambda",
            algoname2 = "actor_critic_TD_modulatedActor_1degree",
            fitType = "pop"
            )
    save_path1 = makehomesavepath("M_signalsInTime")
    mkdir(save_path1)
    save_path2 = makehomesavepath("M_signalsVsEachOther")
    mkdir(save_path2)
    save_path3 = makehomesavepath("M_signalsDifferenceInTime")
    mkdir(save_path3)

    for i in 1:21
        df_signal1, df_signal2 = getSubjSignalofTwoAlgos(subjID = i,
                results_fit_path = results_fit_path,
                signalname1 = signalname1,
                signalname2 = signalname2,
                algoname1 = algoname1,
                algoname2 = algoname2,
                fitType = fitType)
        df_surprise = getSubjSurpriseEvents(subjID = i)

        plotSubjSignalsInTime(df_signal1, df_signal2, df_surprise,
                         subjID = i,
                         scalesurprise = 3, #1.5, #3, #1.3, #1,#3, #1, #10, #1,
                         signalname1 = signalname1,
                         signalname2 = signalname2,
                         save_path = save_path1)
        plotSubjSignalsVsEachOther(df_signal1, df_signal2,
                        subjID = i,
                        signalname1 = signalname1,
                        signalname2 = signalname2,
                        save_path = save_path2)
        # plotSubjSignalDifferenceInTime(df_signal1, df_signal2, df_surprise,
        #                 subjID = i,
        #                 signalname1 = signalname1,
        #                 signalname2 = signalname2,
        #                 scalesurprise = 0.7,
        #                 save_path = save_path3)
    end
end

function getSubjSignalofTwoAlgos(;subjID = 3,
        results_fit_path =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/fits"),
        signalname1 = "SPE",
        signalname2 = "Surprise",
        algoname1 = "hybrid_actorcritic_FWD",
        algoname2 = "actor_critic_TD_modulatedActor_1degree",
        fitType = "pop",
        )
    df_signal1 = getSubjSignal(subjID = subjID,
                    results_fit_path = results_fit_path,
                    signalname = signalname1,
                    fitType = fitType,
                    algoname = algoname1)
    df_signal2 = getSubjSignal(subjID = subjID,
                    results_fit_path = results_fit_path,
                    signalname = signalname2,
                    fitType = fitType,
                    algoname = algoname2)

    df_signal1, df_signal2
end

function getSubjSignal(;subjID = 3,
                results_fit_path =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/fits"),
                signalname = "SPE",
                algoname = "actor_critic_TD_modulatedActor_1degree",
                fitType = "pop",
                )
    df_signal = CSV.read(joinpath(results_fit_path, algoname, "maxLL_fit",
                    "subj_"*string(subjID), signalname*"_"*fitType*".csv"),
                    header=false)
    rename!(df_signal, Dict(:Column1 => Symbol(algoname)))
end
function getSubjSurpriseEvents(;subjID = 3,
                results_fit_path =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/all_events"),
                eventname = "is_surprise_trial")
    df_event = getSubjEvent(;subjID = subjID,
                    results_fit_path = results_fit_path,
                    eventname = eventname)
    rename!(df_event, Dict(:Column1 => Symbol("is_surprise_trial")))
end
function getSubjEvent(;subjID = 3,
                results_fit_path =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/all_events"),
                eventname = "is_surprise_trial")
    df_event = CSV.read(joinpath(results_fit_path,
                    "subj_"*string(subjID), eventname*".csv"),
                    header=false)
end

function plotSubjSignalsInTime(df_signal1, df_signal2;
                    subjID = 3,
                    signalname1 = "SPE",
                    signalname2 = "Surprise",
                    save_path = "./",
                    )
        legendentry1 = replace(signalname1, "_"=>"") * ": " * replace(string(names(df_signal1)[1]), "_"=>"") # Remove "_" cause latex problem
        legendentry2 = replace(signalname2, "_"=>"") * ": " * replace(string(names(df_signal2)[1]), "_"=>"") # Remove "_" cause latex problem

        ax = Axis( #1:size(df_signal1,1)
            [Plots.Scatter(1:size(df_signal1,1), df_signal1[:,1],
                    legendentry = legendentry1,
                    markSize=1),
                    # , style = "violet!95!black"),
            Plots.Scatter(1:size(df_signal2,1), df_signal2[:,1],
                    legendentry = legendentry2,
                    markSize=1)],
           title = replace(signalname1, "_"=>"") *" and "*replace(signalname2, "_"=>""),
           # legendStyle = "{at={(1.05,1.0)},anchor=north west}")
           legendStyle = "{at={(0.3,-0.1)},anchor=north west}",style = "grid=both")
           # ymin = -0.1, ymax = 0.1)
        save(joinpath(save_path, "subj_"*string(subjID)*"_"*signalname1*"_"*signalname2*".pdf"), ax)
end
function plotSubjSignalsInTime(df_signal1, df_signal2, df_surprise;
                    subjID = 3,
                    scalesurprise = 1, #10
                    signalname1 = "SPE",
                    signalname2 = "Surprise",
                    save_path = "./",
                    )

        legendentry1 = L"SPE$_{\mathrm{FW}}$"

        legendentry2 = L"$\textbf{S}_{\mathrm{BF}}$"

        signal1 = df_signal1[:,1]
        signal2 = df_signal2[:,1] #.* 100

        ax = Axis( #1:size(df_signal1,1)
            [Plots.Scatter(1:size(df_signal1,1), signal1,
                    # mark = *,
                    legendentry = legendentry1,
                    markSize=1,
                    style = "orange, mark options={fill=orange}")
                    # scatterClasses=sc1)
            Plots.Scatter(1:size(df_signal2,1), signal2, #df_signal2[:,1],
                    legendentry = legendentry2,
                    # legendentry = legendentry2 * "Norm",
                    markSize=1,
                    style = "violet!95!black, mark options={fill=violet!95!black}")
                    # This is correct, cause it shifts surprise trials one step earlier (which correct given the way things are calculated
                    # and are encoded in all_data, all_stats = RL_Fit.load_SARSPEI(data_path))
            Plots.Linear(1:size(df_surprise,1)-1, scalesurprise .* df_surprise[2:end,1], style=:ycomb, mark = :none)
                    ],
           legendStyle = "{at={(0.01,0.95)},anchor=north west}",
           # ymin = -0.1, ymax = 1.,
           ymax = 3.1,
           # ymax = 1.6,
           xlabel= "time steps",
           style="ymajorgrids=true")
        save(joinpath(save_path, "subj_"*string(subjID)*"_"*signalname1*"_"*signalname2*".pdf"), ax)
end

function plotSubjSignalsVsEachOther(df_signal1, df_signal2;
                    subjID = 3,
                    signalname1 = "SPE",
                    signalname2 = "Surprise",
                    save_path = "./",
                    )

        legendentry1 = L"SPE$_{\mathrm{FW}}$ - Hybrid Actor-critic"
        legendentry2 = L"$\textbf{S}_{\mathrm{BF}}$ - Surprise Actor-critic"

        signal1 = df_signal1[:,1]
        signal2 = df_signal2[:,1]

        # # --- Theory curve
        x = collect(0.:0.05:0.95)
        y = 1. ./ (7. .* (1. .- x) )

        ax = Axis(
            [Plots.Scatter(signal1, signal2,
                    # legendentry = legendentry1,
                    markSize=1,
                    style = "black, mark options={fill=black}"
                    )
                    # ],
            ##### Plot theory curve of Sbf vs SPE
            Plots.Linear(x, y,
                    # legendentry = legendentry1,
                    # markSize=1,
                    style = "red, thick, mark = :none"
                    )
                    ],
           xlabel=legendentry1,
           ylabel=legendentry2,
           xmin = 0.,
           ymin = 0.,
           ymax=3., # for Surprise
           xmax=1.03,
           style = "grid=both"
           )
           # ymin = -0.1, ymax = 0.1)
        save(joinpath(save_path, "subj_"*string(subjID)*"_"*signalname1*"_"*signalname2*".pdf"), ax)
end

function plotSubjSignalDifferenceInTime(df_signal1, df_signal2;
                    subjID = 3,
                    signalname1 = "SPE",
                    signalname2 = "Surprise",
                    save_path = "./",
                    )
        legendentry1 = replace(string(names(df_signal1)[1]), "_"=>"") # Remove "_" cause latex problem
        legendentry2 = replace(string(names(df_signal2)[1]), "_"=>"") # Remove "_" cause latex problem

        ax = Axis(
            [Plots.Scatter(1:length(df_signal1[:,1]) , df_signal1[:,1] .- df_signal2[:,1],
                    markSize=1)],
           title = "Difference " * replace(signalname1, "_"=>"") *" - "*replace(signalname2, "_"=>""),
           xlabel="Time",
           ylabel="Difference",
           style = "grid=both"
           )

        save(joinpath(save_path, "subj_"*string(subjID)*"_"*signalname1*"_"*signalname2*".pdf"), ax)
end
function plotSubjSignalDifferenceInTime(df_signal1, df_signal2, df_surprise;
                    subjID = 3,
                    signalname1 = "RPE",
                    signalname2 = "RPE",
                    scalesurprise = 1.,
                    save_path = "/Users/vasia/Desktop/",
                    )
        legendentry1 = replace(signalname1, "_"=>"") * ": " * replace(string(names(df_signal1)[1]), "_"=>"") # Remove "_" cause latex problem
        legendentry2 = replace(signalname2, "_"=>"") * ": " * replace(string(names(df_signal2)[1]), "_"=>"") # Remove "_" cause latex problem

        ax = Axis(
            [Plots.Scatter(1:length(df_signal1[:,1]) , df_signal1[:,1] .- df_signal2[:,1],
                    # legendentry = legendentry1,
                    markSize=1),
            # This is correct, cause it shifts surprise trials one step earlier (which correct given the way things are calculated
            # and are encoded in all_data, all_stats = RL_Fit.load_SARSPEI(data_path))
            Plots.Linear(1:size(df_surprise,1)-1, scalesurprise .* df_surprise[2:end,1], style=:ycomb, mark = :none)
            ],
           title = "Difference " * replace(signalname1, "_"=>"") *" - "*replace(signalname2, "_"=>""),
           xlabel="Time",
           ylabel="Difference",
           style = "grid=both"
           )

        save(joinpath(save_path, "subj_"*string(subjID)*"_"*signalname1*"_"*signalname2*".pdf"), ax)
end