text:
using MATLAB

calibrate(file, xyt) = last(splitext(file)) == ".csv" ? calibrate_csv(file, xyt) : calibrate_mat(file, xyt)

function calibrate_csv(csvfile, xyt::P) where {P <: AbstractPeriod}
    n = size(xyt.data,1)
    xy1 = [xyt.data[:,1:2] ones(n)]
    tform = readdlm(csvfile, ',', Float64)
    xy2 = xy1*tform
    xy3 = xy2[:,1:2]./xy2[:,3]
    P([xy3 xyt.data[:,3]])
end

function spawnmatlab(targetname, check, intrinsic, extrinsic)
    mat"""
    warning('off','all')
    [imagePoints, boardSize, imagesUsed] = detectCheckerboardPoints($intrinsic);
    $kept = 1:length($intrinsic);
    $kept = $kept(imagesUsed);
    extrinsicI = imread($extrinsic);
    sz = size(extrinsicI);
    worldPoints = generateCheckerboardPoints(boardSize, $check);
    %%
    params = estimateCameraParameters(imagePoints, worldPoints, 'ImageSize', sz, 'EstimateSkew', true, 'NumRadialDistortionCoefficients', 3, 'EstimateTangentialDistortion', true, 'WorldUnits', 'cm');
    n = size(imagePoints, 3);
    errors = zeros(n,1);
    for i = 1:n
        [R,t] = extrinsics(imagePoints(:,:,i), worldPoints, params);
        newWorldPoints = pointsToWorld(params, R, t, imagePoints(:,:,i));
        errors(i) = mean(vecnorm(worldPoints - newWorldPoints, 1, 2));
    end
    kill = errors > 1;
    while any(kill)
        imagePoints(:,:,kill) = [];
        $kept(kill) = [];
        params = estimateCameraParameters(imagePoints, worldPoints, 'ImageSize', sz, 'EstimateSkew', true, 'NumRadialDistortionCoefficients', 3, 'EstimateTangentialDistortion', true, 'WorldUnits', 'cm');
        n = size(imagePoints, 3);
        errors = zeros(n,1);
        for i = 1:n
            [R,t] = extrinsics(imagePoints(:,:,i), worldPoints, params);
            newWorldPoints = pointsToWorld(params, R, t, imagePoints(:,:,i));
            errors(i) = mean(vecnorm(worldPoints - newWorldPoints, 1, 2));
        end
        kill = errors > 1;
    end
    $mean_error = mean(errors);
    %%
    MinCornerMetric = 0.15;
    xy = detectCheckerboardPoints(extrinsicI, 'MinCornerMetric', MinCornerMetric);
    MinCornerMetric = 0.;
    while size(xy,1) ~= size(worldPoints, 1)
        MinCornerMetric = MinCornerMetric + 0.05;
        xy = detectCheckerboardPoints(extrinsicI, 'MinCornerMetric', MinCornerMetric);
    end
    [R,t] = extrinsics(xy, worldPoints, params);
    save($targetname, 'params', 'R', 't')
    """
    kill = setdiff(1:length(intrinsic), kept)
    for i in kill
        rm(intrinsic[i])
    end
    mean_error
end

function spawnmatlab(targetname, check, extrinsic)
    mat"""
    warning('off','all')
    [imagePoints, boardSize] = detectCheckerboardPoints($extrinsic);
    worldPoints = generateCheckerboardPoints(boardSize, $check);
    tform_ = fitgeotrans(imagePoints, worldPoints, 'projective');
    $tform = tform_.T;
    %%
    [x, y] =  transformPointsForward(tform_, imagePoints(:,1), imagePoints(:,2));
    $mean_error = mean(vecnorm(worldPoints - [x, y], 2, 2))
    """
    writedlm(targetname, tform, ',')
    mean_error
end

function calibrate_mat(matfile, xyt::P) where {P <: AbstractPeriod}
    xy = xyt.data[:,1:2]
    mat"""
    a = load($matfile);
    $xy = pointsToWorld(a.params, a.R, a.t, $xy);
    """
    P([xy xyt.data[:,3]])
end

function build_calibration(coffeesource, calibId, c::Calibration{Temporal{V1, P}, Temporal{V2, I}}) where {V1 <: AbstractTimeLine, V2 <: AbstractTimeLine, P <: Prolonged, I <: Instantaneous} # moving
    path = joinpath(coffeesource, "calibration_images", string(calibId))
    mkpath(path)
    extract(joinpath(path, "extrinsic.png"), c.extrinsic, coffeesource)
    extract(path, c.intrinsic, coffeesource)
    intrinsic = [joinpath(path, f) for f in readdir(path) if f ≠ "extrinsic.png"]
    extrinsic = joinpath(path, "extrinsic.png")
    mkpath(joinpath(coffeesource, "autocalib"))
    targetname = joinpath(coffeesource, "autocalib", "$calibId.mat")
    ϵ = spawnmatlab(targetname, c.board.checker_width_cm, intrinsic, extrinsic)
    (filename = targetname, error = ϵ)
end

function build_calibration(coffeesource, calibId, c::Calibration{Missing, Temporal{V2, I}}) where {V1 <: AbstractTimeLine, V2 <: AbstractTimeLine, I <: Instantaneous}
    path = joinpath(coffeesource, "calibration_images", string(calibId))
    mkpath(path)
    extract(joinpath(path, "extrinsic.png"), c.extrinsic, coffeesource)
    extrinsic = joinpath(path, "extrinsic.png")
    mkpath(joinpath(coffeesource, "autocalib"))
    targetname = joinpath(coffeesource, "autocalib", "$calibId.csv")
    ϵ = spawnmatlab(targetname, c.board.checker_width_cm, extrinsic)
    (filename = targetname, error = ϵ)
end

ext: jl


text:
#
# CART-related types and functions
#

struct Cart <: Wrapper
    pyo::PyObject

    function Cart(x::DataSet, y::Vector; threshold = nothing, threshold_type = ">",
                         include = nothing, exclude = nothing, kwargs...)

        pandasDF = pandas_dataframe(x; include = include, exclude = exclude)

        cart = rhodium.Cart(pandasDF, y; threshold = threshold, threshold_type = threshold_type,
                            kwargs...)
        new(cart)
    end
end

function show_tree(c::Cart; kwargs...)
    fig = c.pyo.show_tree(kwargs...)
    return fig
end

function print_tree(c::Cart; coi = nothing, all = true, kwargs...)
    c.pyo.print_tree(;coi = coi, all = all, kwargs...)
    nothing
end

function save(c::Cart, file, format = "png"; kwargs...)
    c.pyo[:save](file, format; kwargs...)
end

function save_pdf(c::Cart, file; feature_names = nothing, kwargs...)
    save(c, file, "pdf"; feature_names = feature_names, kwargs...)
end

function save_png(c::Cart, file; feature_names = nothing, kwargs...)
    save(c, file, "png"; feature_names = feature_names, kwargs...)
end


text:
"""
    Trajectory of a dipole in an electric field
    ===========================================
"""

function voltages_t(voltages::Vector{Float64}, t::Float64)
    """ electrode voltages at time t
    """
    return voltages
end

function voltages_t(voltages::Function, t::Float64)
    """ electrode voltages at time t
    """
    return voltages(t)
end

function kinetic_energy(v::Vector{Float64}, mass::Float64)
    """ kinetic energy of a moving particle
    """
    return 0.5 * mass * (v[1]^2.0 + v[2]^2.0 + v[3]^2.0)
end

function potential_energy(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, xg::Float64, yg::Float64, zg::Float64, dipole::Float64; get_famp::Bool=false)
    """ potential energy of a static dipole
    """
    famp = amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg)
    if get_famp
        return dipole * famp, famp
    else
        return dipole * famp
    end
end

function potential_energy(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, xg::Float64, yg::Float64, zg::Float64, dipole::Function; get_famp::Bool=false)
    """ potential energy of a field-dependent dipole
    """
    famp = amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg)
    dpm = dipole(famp)
    if get_famp
        return dpm * famp, famp
    else
        return dpm * famp
    end
end

function force(fa::FastAdjust3D, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, dipole_f::Function)
    """ gradient of the potential energy of a dipole in an inhomegeous electric field (3D)
    """
    xg, yg, zg = grid_r(fa, r)
    famp = amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg)
    dipole = dipole_f(famp)
    fx = -(potential_energy(fa, voltages, t, xg + 0.5, yg, zg, dipole) -
           potential_energy(fa, voltages, t, xg - 0.5, yg, zg, dipole)) / fa.dx
    fy = -(potential_energy(fa, voltages, t, xg, yg + 0.5, zg, dipole) -
           potential_energy(fa, voltages, t, xg, yg - 0.5, zg, dipole)) / fa.dy
    fz = -(potential_energy(fa, voltages, t, xg, yg, zg + 0.5, dipole) -
           potential_energy(fa, voltages, t, xg, yg, zg - 0.5, dipole)) / fa.dz
    return Vector([fx, fy, fz])
end

function force(fa::FastAdjust2D, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, dipole_f::Function)
    """ gradient of the potential energy of a dipole in an inhomegeous electric field (2D)
    """
    xg, yg, zg = grid_r(fa, r)
    famp = amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg)
    dipole = dipole_f(famp)
    fx = -(potential_energy(fa, voltages, t, xg + 0.5, yg, zg, dipole) -
           potential_energy(fa, voltages, t, xg - 0.5, yg, zg, dipole)) / fa.dx
    fy = -(potential_energy(fa, voltages, t, xg, yg + 0.5, zg, dipole) -
           potential_energy(fa, voltages, t, xg, yg - 0.5, zg, dipole)) / fa.dy
    fz = 0.0
    return Vector([fx, fy, fz])
end

function force(fa::FastAdjust3D, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, dipole::Float64)
    """ gradient of the potential energy of a dipole in an inhomegeous electric field (3D)
    """
    xg, yg, zg = grid_r(fa, r)
    fx = -(potential_energy(fa, voltages, t, xg + 0.5, yg, zg, dipole) -
           potential_energy(fa, voltages, t, xg - 0.5, yg, zg, dipole)) / fa.dx
    fy = -(potential_energy(fa, voltages, t, xg, yg + 0.5, zg, dipole) -
           potential_energy(fa, voltages, t, xg, yg - 0.5, zg, dipole)) / fa.dy
    fz = -(potential_energy(fa, voltages, t, xg, yg, zg + 0.5, dipole) -
           potential_energy(fa, voltages, t, xg, yg, zg - 0.5, dipole)) / fa.dz
    return Vector([fx, fy, fz])
end

function force(fa::FastAdjust2D, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, dipole::Float64)
    """ gradient of the potential energy of a dipole in an inhomegeous electric field (2D)
    """
    xg, yg, zg = grid_r(fa, r)
    fx = -(potential_energy(fa, voltages, t, xg + 0.5, yg, zg, dipole) -
           potential_energy(fa, voltages, t, xg - 0.5, yg, zg, dipole)) / fa.dx
    fy = -(potential_energy(fa, voltages, t, xg, yg + 0.5, zg, dipole) -
           potential_energy(fa, voltages, t, xg, yg - 0.5, zg, dipole)) / fa.dy
    fz = 0.0
    return Vector([fx, fy, fz])
end

function acceleration(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, dipole::Union{Float64, Function}, mass::Float64)
    """ accelleration of a dipole in an inhomegeous electric field
    """
    return  force(fa, voltages, t, r, dipole) ./ mass
end

function euler(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, v::Vector{Float64}, a::Vector{Float64}, dipole::Union{Float64, Function}, mass::Float64, dt::Float64)
    """ Euler method of integration
    """
    a = acceleration(fa, voltages, t, r, dipole, mass)
    rn = r .+ v .* dt
    vn = v .+ a .* dt
    return rn, vn, a
end

function leapfrog(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, v::Vector{Float64}, a::Vector{Float64}, dipole::Union{Float64, Function}, mass::Float64, dt::Float64)
    """ leapfrog method of integration
    """
    if any(isnan.(a))
        a = acceleration(fa, voltages, t, r, dipole, mass)
    end
    rn = r .+ v .* dt + a .* (0.5 * dt^2.0)
    an = acceleration(fa, voltages, t + dt, rn, dipole, mass)
    vn = v .+ (a .+ an) .* (0.5 * dt)
    return rn, vn, an
end

function rk4(fa::FastAdjust, voltages::Union{Function, Vector{Float64}}, t::Float64, r::Vector{Float64}, v::Vector{Float64}, a::Vector{Float64}, dipole::Union{Float64, Function}, mass::Float64, dt::Float64)
    """ fourth order Runge-Kutta method of integration
    """
    # initial position
    a = acceleration(fa, voltages, t, r, dipole, mass)
    # half-step
    r1 = r .+ (0.5 * dt) .* v
    v1 = v .+ (0.5 * dt) .* a
    a1 = acceleration(fa, voltages, t + 0.5 * dt, r1, dipole, mass)
    # half-step again
    r2 = r .+ (0.5 * dt) .* v1
    v2 = v .+ (0.5 * dt) .* a1
    a2 = acceleration(fa, voltages, t + 0.5 * dt, r2, dipole, mass)
    # full step
    r3 = r .+ dt .* v2
    v3 = v .+ dt .* a2
    a3 = acceleration(fa, voltages, t + dt, r3, dipole, mass)
    # next position
    r4 = r .+ (dt / 6.0) .* (v .+ 2.0 .* v1 .+ 2.0 .* v2 .+ v3)
    v4 = v .+ (dt / 6.0) .* (a .+ 2.0 .* a1 .+ 2.0 .* a2 .+ a3)
    return r4, v4, a3
end

function traj(fa::FastAdjust, voltages::Union{Function, Vector{Float64}},
              t0::Float64, r0::Vector{Float64}, v0::Vector{Float64},
              dipole::Union{Float64, Function}, mass::Float64, dt::Float64, method::Function;
              max_field::Float64=NaN, max_t::Float64=1.0, max_iterations::Int64=Int64(1e6), df::Bool=false)
    """ full trajectory
    """
    # initialise
    i = 1
    t, r, v = t0, r0, v0
    a = Vector([NaN, NaN, NaN])
    xg, yg, zg = grid_r(fa, r)
    KE = kinetic_energy(v, mass)
    PE = potential_energy(fa, voltages, t, xg, yg, zg, dipole)
    famp = amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg)
    result = [t r[1] r[2] r[3] KE PE famp;]
    # step-by-step trajectory
    while i < max_iterations
        try
            r, v, a = method(fa, voltages, t, r, v, a, dipole, mass, dt)
            t += dt
        catch
            break
        end
        if any(isnan, r) | any(isnan, v)
            break
        end
        # grid position
        xg, yg, zg = grid_r(fa, r)
        # energy
        KE = kinetic_energy(v, mass)
        PE, famp = potential_energy(fa, voltages, t, xg, yg, zg, dipole, get_famp=true)
        # checks
        if electrode_g(fa, xg, yg, zg)
            # hit an electrode or left the pa
            break
        elseif ~isnan(max_field) && famp > max_field
            # ionisation field
            break
        elseif ~isnan(max_t) && t > max_t
            # time's up!
            break
        end
        # record
        result = vcat(result, [t r[1] r[2] r[3] KE PE famp])
        # next step
        i += 1
    end
    # output
    if df
        result = DataFrame(result)
        rename!(result, [:t, :x, :y, :z, :KE, :PE, :famp]);
    end
    return result
end

function traj_last(fa::FastAdjust, voltages::Union{Function, Vector{Float64}},
                   t0::Float64, r0::Vector{Float64}, v0::Vector{Float64},
                   dipole::Union{Float64, Function}, mass::Float64, dt::Float64, method::Function;
                   max_field::Float64=NaN, max_t::Float64=NaN, max_iterations::Int64=Int64(1e6))
    """ last point of a trajectory
    """
    # initialise
    i = 1
    t, r, v = t0, r0, v0
    a = Vector([NaN, NaN, NaN])
    result = [t r[1] r[2] r[3] v[1] v[2] v[3];]
    # step-by-step trajectory
    while i < max_iterations
        try
            r, v, a = method(fa, voltages, t, r, v, a, dipole, mass, dt)
            t += dt
        catch
            break
        end
        if any(isnan, r) | any(isnan, v)
            break
        end
        # grid position
        xg, yg, zg = grid_r(fa, r)
        # checks
        if electrode_g(fa, xg, yg, zg)
            # hit an electrode or left the pa
            break
        elseif ~isnan(max_field) && amp_field_g(fa, voltages_t(voltages, t), xg, yg, zg) > max_field
            # ionisation field
            break
        elseif ~isnan(max_t) && t > max_t
            # time's up!
            break
        end
        # record
        result = [t r[1] r[2] r[3] v[1] v[2] v[3];]
        # next step
        i += 1
    end
    # output
    return result
end


text:
# The "default" Associative (a Dict)
struct Dic{I, T} <: Assoc{I, T}
    d::Dict{I, T}
end

Dic(pairs::Pair...) = Dic(Dict(pairs...))
Dic{I, T}(pairs::Pair{>:I, >:T}...) where {I, T} = Dic(Dict(pairs...))
Dic() = Dic{Any, Any}()
Dic{I, T}() where {I, T} = Dic{I, T}(Dict{I, T}())

indices(d::Dic) = keys(d.d) # TODO: should this return an `AbstractSet` or something?
getindex(d::Dic{I}, i::I) where {I} = d.d[i]
@propagate_inbounds function setindex!(d::Dic{I,T}, v::T, i::I) where {I, T}
    d.d[i] = v
    return v
end

# This is a hack to get `similar` to work at all (probably not necessary with some more
# work / a new inner constructor for `Dict`)
default(::Type{T}) where {T} = isbits(T) ? Ref{T}()[] : zero(T)
default(::Type{String}) = ""
default(::Type{Array{T, N}}) where {T, N} = Array{T, N}()

# The default `similar` returns a `Dic`
function similar(::Assoc, ::Type{T}, inds) where {T}
    d = Dict{eltype(inds), T}()
    for i in inds
        d[i] = default(T)
    end
    return Dic(d)
end

empty(::Assoc, ::Type{I}, ::Type{T}) where {T, I} = Dic{I,T}(Dict{I,T}())


text:
export TABound

mutable struct TABound{T}
    strict::Bool
    value::T
    infinite::Bool
    function TABound(s,v::T) where T
        b=new{T}(s,v,false)
        b.value=v
        b
    end
end
function Base.setproperty!(b::TABound{T}, name::Symbol, value) where T
    if name == :strict
        setfield!(b, :strict, value)
    elseif name == :value
        if isinf(value) || value==typemax(value) || value==typemin(value)
            if value > 0
                setfield!(b, :value, typemax(T))
            else
                setfield!(b, :value, typemin(T))
            end
            setfield!(b, :infinite, true)
        else
            setfield!(b, :value, value)
            setfield!(b, :infinite, false)
        end
    else
        throw(ArgumentError("Can't set property '$(name)' of '$(b)'!"))
    end
end

Base.isinf(b::TABound) = b.infinite
Base.iszero(b::TABound) = iszero(b.value)
Base.hash(b::TABound) = hash((hash(b.strict), hash(b.value)))
Base.:(==)(b1::TABound, b2::TABound) = ((b1.value == b2.value) && (b1.strict == b2.strict)) || (isinf(b1) && isinf(b2))
Base.isequal(b1::TABound, b2::TABound) = b1==b2
Base.isless(b1::TABound{T}, b2::TABound) where {T} = (!isinf(b1) && isinf(b2) && b2.value > zero(T)) || (b1.value < b2.value) || ((b1.value == b2.value) && b1.strict && !b2.strict )
(Base.zero(::Type{TABound{T}}) where T) = TABound(false, zero(T))
(Base.typemax(::Type{TABound{T}}) where T) = TABound(false, typemax(T))
(Base.typemin(::Type{TABound{T}}) where T) = TABound(false, typemin(T))
function Base.:+(b1::TABound{T}, b2::TABound) where T
    if isinf(b1) && isinf(b2) && sign(b1.value)!=sign(b2.value)
        throw(ArgumentError("Sum of bounds undefined! (∞ - ∞)"))
    elseif isinf(b1) && b1.value > 0 || isinf(b2) && b2.value > 0
        typemax(TABound{T})
    elseif isinf(b1) && b1.value < 0 || isinf(b2) && b2.value < 0
        typemin(TABound{T})
    else
        TABound(b1.strict || b2.strict, b1.value + b2.value)
    end
end
Base.:-(c::TABound) = TABound(c.strict, -c.value)
Base.:-(c1::TABound, c2::TABound) = c1 + (-c2)

Base.copy(b::TABound) = TABound(b.s,b.v)
Base.deepcopy(b::TABound) = TABound(deepcopy(b.s),deepcopy(b.v))
Base.iterate(x::TABound) = (x, nothing)
Base.iterate(::TABound, ::Any) = nothing
Base.length(x::TABound) = 1
Base.Broadcast.broadcastable(b::TABound) = Ref(b)
Base.adjoint(c::TABound) = c

function Base.show(io::IO, M::MIME"text/plain", b::TABound)
    s = if isinf(b) && b.value > 0
        "∞"
    elseif isinf(b) && b.value ≤ 0
        "-∞"
    elseif b.strict
        "< "*repr(b.value)
    else
        "≤ "*repr(b.value)
    end
    print(io, s)
end


text:
using Joseki, JSON, HTTP

### Create some endpoints

# This function takes two numbers x and y from the query string and returns x^y
# In this case they need to be identified by name and it should be called with
# something like 'http://localhost:8000/pow/?x=2&y=3'
function pow(req::HTTP.Request)
    j = HTTP.queryparams(HTTP.URI(req.target))
    has_all_required_keys(["x", "y"], j) || return error_responder(req, "You need to specify values for x and y!")
    # Try to parse the values as numbers.  If there's an error here the generic
    # error handler will deal with it.
    x = parse(Float32, j["x"])
    y = parse(Float32, j["y"])
    json_responder(req, x^y)
end

# This function takes two numbers n and k from a JSON-encoded request
# body and returns binomial(n, k)

# This function takes two numbers x and y from the query string and returns x^y
# In this case they need to be identified by name and it should be called with
# something like 'http://localhost:8000/pow/?x=2&y=3'
function returnmx(req::HTTP.Request)
    j = HTTP.queryparams(HTTP.URI(req.target))
    has_all_required_keys(["x", "y"], j) || return error_responder(req, "You need to specify values for x and y!")
    # Try to parse the values as numbers.  If there's an error here the generic
    # error handler will deal with it.
    x = parse(Int, j["x"])
    y = parse(Int, j["y"])
    A=rand(x,y)
    json_responder(req, A)
end

# This function takes two numbers n and k from a JSON-encoded request
# body and returns binomial(n, k)

function bin(req::HTTP.Request)
    j = try
        body_as_dict(req)
    catch err
        return error_responder(req, "I was expecting a json request body!")
    end
    has_all_required_keys(["n", "k"], j) || return error_responder(req, "You need to specify values for n and k!")
    json_responder(req, binomial(j["n"],j["k"]))
end

### Create and run the server

# Make a router and add routes for our endpoints.
endpoints = [
    (pow, "GET", "/pow"),
    (returnmx, "GET", "/returnmx"),
    (bin, "POST", "/bin"),
    (req -> req.response, "OPTIONS", "*")
]
r = Joseki.router(endpoints)

# If there is a PORT environment variable defined us it, otherwise use 8000
haskey(ENV, "PORT") ? port = ENV["PORT"] : port = 8000
# Fire up the server, binding to all ips
HTTP.serve(r, "0.0.0.0", port; verbose=false)

max_stars_repo_name: wadehenning/fda-app

max_issues_repo_name: wadehenning/fda-app

max_forks_repo_name: wadehenning/fda-app


text:
using Random

@testset "predict" begin
    Random.seed!(100)

    @model function linear_reg(x, y, σ = 0.1)
        β ~ Normal(0, 1)

        for i ∈ eachindex(y)
            y[i] ~ Normal(β * x[i], σ)
        end
    end

    @model function linear_reg_vec(x, y, σ = 0.1)
        β ~ Normal(0, 1)
        y ~ MvNormal(β .* x, σ)
    end

    f(x) = 2 * x + 0.1 * randn()

    Δ = 0.1
    xs_train = 0:Δ:10; ys_train = f.(xs_train);
    xs_test = [10 + Δ, 10 + 2 * Δ]; ys_test = f.(xs_test);

    # Infer
    m_lin_reg = linear_reg(xs_train, ys_train);
    chain_lin_reg = sample(m_lin_reg, NUTS(100, 0.65), 200);

    # Predict on two last indices
    m_lin_reg_test = linear_reg(xs_test, fill(missing, length(ys_test)));
    predictions = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    ys_pred = vec(mean(Array(group(predictions, :y)); dims = 1))

    @test sum(abs2, ys_test - ys_pred) ≤ 0.1

    # Ensure that `rng` is respected
    predictions1 = let rng = MersenneTwister(42)
        predict(rng, m_lin_reg_test, chain_lin_reg[1:2])
    end
    predictions2 = let rng = MersenneTwister(42)
        predict(rng, m_lin_reg_test, chain_lin_reg[1:2])
    end
    @test all(Array(predictions1) .== Array(predictions2))

    # Predict on two last indices for vectorized
    m_lin_reg_test = linear_reg_vec(xs_test, missing);
    predictions_vec = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)
    ys_pred_vec = vec(mean(Array(group(predictions_vec, :y)); dims = 1))

    @test sum(abs2, ys_test - ys_pred_vec) ≤ 0.1

    # Multiple chains
    chain_lin_reg = sample(m_lin_reg, NUTS(100, 0.65), MCMCThreads(), 200, 2);
    m_lin_reg_test = linear_reg(xs_test, fill(missing, length(ys_test)));
    predictions = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    @test size(chain_lin_reg, 3) == size(predictions, 3)

    for chain_idx in MCMCChains.chains(chain_lin_reg)
        ys_pred = vec(mean(Array(group(predictions[:, :, chain_idx], :y)); dims = 1))
        @test sum(abs2, ys_test - ys_pred) ≤ 0.1
    end

    # Predict on two last indices for vectorized
    m_lin_reg_test = linear_reg_vec(xs_test, missing);
    predictions_vec = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    for chain_idx in MCMCChains.chains(chain_lin_reg)
        ys_pred_vec = vec(mean(
            Array(group(predictions_vec[:, :, chain_idx], :y));
            dims = 1
        ))
        @test sum(abs2, ys_test - ys_pred_vec) ≤ 0.1
    end
end


text:
using StructJuMP

firststage = StructuredModel()
@variable(firststage, x[1:2])
@objective(firststage, :Min, (x[1]+x[2])^2)
@NLconstraint(firststage, x[1] * x[2] == 10)

for scen in 1:2
    bl = StructuredModel(parent=firststage, id=scen)
    @variable(bl, y)
    @constraint(bl, x[2]^2 + x[1]*y ≤ 5)
    @NLobjective(bl, :Min, (x[1]+x[2])*y)
end


text:
export Pathfinding

"""
    Pathfinding
Submodule containing functionality for path-finding based on the A* algorithm.
Currently available only for [`GridSpace`](@ref).

You can enable path-finding and set it's options by passing an instance of a
[`Pathfinding.Pathfinder`](@ref) struct to the `pathfinder`
parameter of the [`GridSpace`](@ref) constructor. During the simulation, call
[`Pathfinding.set_target!`](@ref) to set the target destination for an agent.
This triggers the algorithm to calculate a path from the agent's current position to the one
specified. You can alternatively use [`Pathfinding.set_best_target!`](@ref) to choose the
best target from a list. Once a target has been set, you can move an agent one step along
its precalculated path using the [`move_along_route!`](@ref) function.

Refer to the [Maze Solver](@ref) and [Mountain Runners](@ref) examples using path-finding
and see the available functions below as well.
"""
module Pathfinding

using Agents
using DataStructures
using LinearAlgebra

include("metrics.jl")
include("grid_pathfinder.jl")

export CostMetric,
    DirectDistance,
    MaxDistance,
    HeightMap,
    Pathfinder,
    delta_cost,
    set_target!,
    set_best_target!,
    heightmap,
    walkmap

end


text:
################################################################################
## ReverseKullbackLeibler
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate{T<:AbstractFloat}(div::RKL, a::AbstractVector{T}, b::AbstractVector{T})
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end

    r = zero(T)

    for i in eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        ui = ai/bi
        if ui > 0
            r += -bi*log(ui) + ai - bi
        else
            r = convert(T, Inf)
            break
        end
    end
    return r
end

function evaluate{T<:AbstractFloat}(div::RKL, a::AbstractVector{T})
    r = zero(T)
    for i in eachindex(a)
        @inbounds ai = a[i]
        if ai > 0
            r += - logmxp1(ai)
        else
            r = convert(T, Inf)
            break
        end
    end
    r
end

#=---------------
Gradient
---------------=#
function gradient{T<:AbstractFloat}(div::RKL, a::T, b::T)
    if a > 0 && b > 0
        u = - b/a + 1.0
    else
        u = convert(T, Inf)
    end
    return u
end

function gradient{T<:AbstractFloat}(div::RKL, a::T)
    if a > 0
        u = - 1.0 / a + 1.0
    else
        u = convert(T, Inf)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian{T<:AbstractFloat}(div::RKL, a::T, b::T)
    if a > 0 && b > 0
        u = b/a^2
    else
        u = convert(T, Inf)
    end
    return u
end

function hessian{T<:AbstractFloat}(div::RKL, a::T)
    if a > 0
        u = 1.0/a^2
    else
        u = convert(T, Inf)
    end
    u
end

################################################################################
## Modified Reverse Kullback-Leibler - MRKL
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate{T<:AbstractFloat}(div::MRKL, a::AbstractVector{T}, b::AbstractVector{T})
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    f0, f1, f2, uϑ = div.m
    r = zero(T)
    for i = eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        ui = ai/bi
        if ui >= uϑ
            r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui > 0 && ui <u₀
            r += -bi*log(ui) + ai - bi
        else
            r = convert(T, Inf)
            break
        end
    end
    return r
end

function evaluate{T<:AbstractFloat}(div::MRKL, a::AbstractVector{T})
    f0, f1, f2, uϑ = div.m
    r = zero(T)
    for i in eachindex(a)
        @inbounds ai = a[i]
        if ai >= uϑ
            r += f0 + f1*(ai-uϑ) + .5*f2*(ai-uϑ)^2
        elseif ai > 0 && ai < uϑ
            r += -log(ai) + ai - 1.0
        else
            r = convert(T, Inf)
            break
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient{T<:AbstractFloat}(div::MRKL, a::T)
    f0, f1, f2, uϑ = div.m
    ui = a
    if ui > u₀
        u = (f1 + f2*(ui-uϑ))
    elseif ui <= uϑ
        u = gradient(div.d, a)
    end
    return u
end

function gradient{T<:AbstractFloat}(div::MRKL, a::T, b::T)
    f0, f1, f2, uϑ = div.m
    ui = a/b
    if ui > u₀
        u = f1 + f2*(ui-uϑ)
    elseif ui <= uϑ
        u = gradient(div.d, a, b)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian{T<:AbstractFloat}(div::MRKL, a::T)
    f0, f1, f2, uϑ = div.m
    if a >= uϑ
       u  = f2
    else
       u = hessian(div.d, a)
    end
    return u
end

function hessian{T<:AbstractFloat}(div::MRKL, a::T, b::T)
    f0, f1, f2, uϑ = div.m
    ui = a/b
    if ui >= u₀
        u  = f2/b
    else
        u = hessian(rkl, a, b)
    end
    return u
end

################################################################################
## Fully Modified Reverse Kullback-Leibler - FMRKL
##
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate{T<:AbstractFloat}(div::FMRKL, a::AbstractVector{T}, b::AbstractVector{T})
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    @inbounds for i = eachindex(a, b)
        ai = a[i]
        bi = a[i]
        ui = ai/bi
        if ui >= uϑ
            r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui <= uφ
            r += (g0 + g1*(ui-uφ) + .5*g2*(ui-uφ)^2)*bi
        else
            r += -bi*log(ui) + ai - bi
        end
    end
    return r
end

function evaluate{T<:AbstractFloat}(div::FMRKL, a::AbstractVector{T})
    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    @inbounds for i = eachindex(a)
        ui = a[i]
        if ui >= uϑ
            r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)
        elseif ui <= uφ
            r += (g0 + g1*(ui-uφ) + .5*g2*(ui-uφ)^2)
        else
            r += -log(ui) + ui - 1.0
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient{T<:AbstractFloat}(div::FMRKL, a::T, b::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    @inbounds for i = eachindex(a, b)
        ai = a[i]
        bi = a[i]
        ui = ai/bi
        if ui >= uϑ
            u = f1 + f2*(ui-uϑ)
        elseif ui <= uφ
            u = g1 + g2*(ui-uφ)
        else
            u = 1-bi/ai
        end
    end
    return u
end

function gradient{T<:AbstractFloat}(div::FMRKL, a::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u = (f1 + f2*(a-uϑ))
    elseif a <= uφ
        u = (g1 + g2*(a-uφ))
    else
        u = 1.0 - 1.0/a
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian{T<:AbstractFloat}(div::FMRKL, a::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u  = f2
    elseif a <= uφ
        u  = g2
    else
        u = 1/(a*a)
    end
    return u
end

function hessian{T<:AbstractFloat}(div::FMRKL, a::T, b::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    ui = a/b
    if ui >= uϑ
        u  = f2/b
    elseif ui <= uφ
        u  = g2/b
    else
        u = b/(a*a)
    end
    return u
end


text:
import NCDatasets
const NC = NCDatasets
import SurfaceFluxes
const SF = SurfaceFluxes
const UF = SurfaceFluxes.UniversalFunctions
using Statistics
using StaticArrays
import Thermodynamics
import ArtifactWrappers
const AW = ArtifactWrappers
import UnPack
const TD = Thermodynamics

using CLIMAParameters: AbstractEarthParameterSet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

#! format: off
PyCLES_output_dataset = AW.ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "PyCLES_output",
    AW.ArtifactFile[
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/johlutwhohvr66wn38cdo7a6rluvz708.nc", filename = "Rico.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/zraeiftuzlgmykzhppqwrym2upqsiwyb.nc", filename = "Gabls.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/toyvhbwmow3nz5bfa145m5fmcb2qbfuz.nc", filename = "DYCOMS_RF01.nc",),
    AW.ArtifactFile(url = "https://caltech.box.com/shared/static/jci8l11qetlioab4cxf5myr1r492prk6.nc", filename = "Bomex.nc",),
    ],
)
#! format: on
PyCLES_output_dataset_path = AW.get_data_folder(PyCLES_output_dataset)

files = ["DYCOMS_RF01.nc", "Bomex.nc", "Rico.nc", "Gabls.nc"];
for f in files
    @info "Casename: $f"

    nt = NC.NCDataset(joinpath(PyCLES_output_dataset_path, f), "r") do data
        prof = data.group["profiles"]
        timeseries = data.group["timeseries"]
        z = Array(data["z_half"])
        u = Array(prof["u_mean"])
        v = Array(prof["v_mean"])
        qt = Array(prof["qt_mean"])
        ql = Array(prof["ql_mean"])
        ρ = Array(prof["rho"])
        p0 = Array(prof["p0"])
        b = Array(prof["buoyancy_mean"])
        θli = Array(prof["thetali_mean"]) # Care with variable names across cases ()?
        T = Array(prof["temperature_mean"]) # Care with variable names across cases ()?
        SHF = Array(timeseries["shf_surface_mean"]) # Care with variable names across cases ()?
        LHF = Array(timeseries["lhf_surface_mean"]) # Care with variable names across cases ()?
        (; z, u, v, qt, ql, ρ, p0, b, θli, T, SHF, LHF)
    end
    UnPack.@unpack z, u, v, qt, ql, ρ, p0, b, θli, T, SHF, LHF = nt

    function getval(X) # Consider only the last timestep
        X[:, end]
    end

    # Data at first interior node (x_in)
    # TODO Make sure that the first node ii = 1 is at the "surface"
    # for the tests to be consistent
    ii = 2
    FT = Float32

    shf = mean(getval(SHF))
    lhf = mean(getval(LHF))

    z_in = z[ii]
    u_in = getval(u)[ii]
    v_in = getval(v)[ii]
    b_in = getval(b)[ii]
    ρ_in = getval(ρ)[ii]
    θ_in = getval(θli)[ii]
    T_in = getval(T)[ii]
    qt_in = getval(qt)[ii]

    # Surface values for variables
    jj = 1
    u_sfc = FT(0)
    v_sfc = FT(0)
    b_sfc = getval(b)[jj]
    ρ_sfc = getval(ρ)[jj]
    qt_sfc = getval(qt)[jj]
    θ_sfc = getval(θli)[jj]
    T_sfc = getval(T)[jj]
    z_sfc = FT(0)

    # Roughness
    ## Roughness lengths
    z0m = FT(0.001)
    z0b = FT(0.001)

    ts_sfc = TD.PhaseEquil_ρθq(param_set, ρ_sfc, θ_sfc, qt_sfc)
    ts_in = TD.PhaseEquil_ρθq(param_set, ρ_in, θ_in, qt_in)

    u_in = SVector{2, FT}(u_in, v_in)
    u_sfc = SVector{2, FT}(u_sfc, v_sfc)

    state_sfc = SF.SurfaceValues(z_sfc, u_sfc, ts_sfc)
    state_in = SF.InteriorValues(z_in, u_in, ts_in)

    kwargs = (; state_in, state_sfc, z0m, z0b)

    # shf and lhf in dycoms are mising /wrong in the file
    if f == "DYCOMS_RF01.nc"
        sc = SF.Fluxes{FT}(; kwargs..., shf = FT(15), lhf = FT(115))
    elseif f == "Bomex.nc"
        sc = SF.FluxesAndFrictionVelocity{FT}(; kwargs..., shf = shf, lhf = lhf, ustar = FT(0.28))
    elseif f == "Rico.nc"
        sc = SF.Coefficients{FT}(; kwargs..., Cd = FT(0.001229), Ch = FT(0.001094))
    elseif f == "Gabls.nc"
        sc = SF.ValuesOnly{FT}(; kwargs...)
    end
    uf = UF.Businger()
    result = SF.surface_conditions(param_set, sc, uf)
end


text:
using WoodburyMatrices
using Test

@testset "Woodbury" begin
    d  = 2 .+ rand(5)
    dl = -rand(4)
    du = -rand(4)
    M = Tridiagonal(dl, d, du)
    F = lu(M)
    U = sprand(5,2,0.2)
    V = sprand(2,5,0.2)
    C = rand(2,2)
    W = Woodbury(F, U, C, V)

    src = rand(5, 8)
    @test W\src ≈ AxisAlgorithms.A_ldiv_B_md(W, src, 1)
    src = rand(5, 5, 5, 5)
    for dim = 1:4
        dest1 = mapslices(x->W\x, copy(src), dims=dim)
        dest2 = similar(src)
        AxisAlgorithms.A_ldiv_B_md!(dest2, W, src, dim)
        @test dest1 ≈ dest2
    end
end

max_stars_repo_name: timholy/FilterMD

max_issues_repo_name: timholy/FilterMD

max_forks_repo_name: timholy/FilterMD


text:
# if VERSION < v"1.1"
#     import Compat:eachcol
# else
    import DataFrames:eachcol
# end

import Base: iterate, length

struct JDFFileColIterator
    jdf::JDFFile
    cols::Vector{Symbol}
end

eachcol(jdf::JDFFile) = JDFFileColIterator(jdf, names(jdf))

Base.length(jdf::JDFFileColIterator) = length(jdf.cols)

Base.iterate(jdf::JDFFileColIterator, state = 1) = begin
    if state > length(jdf.cols)
        return nothing
    end

    res1 = sloadjdf(jdf.jdf, cols=[jdf.cols[state]])[!,1]
    return (res1, state + 1)
end


text:
using Gen
using Test
import Random

include("assignment.jl")
include("lightweight.jl")
include("basic.jl")
include("injective.jl")
include("inference.jl")
include("distribution.jl")

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
abstract type AbstractRHSFunction end
function build_RHS_function end

@inline A_mul_B!(α::T, A::SharedMatrix{T}, B::Vector{T}, β::T, C::Vector{T}) where T<:Union{Float64, Complex{Float64}} = BLAS.gemv!('N', α, A, B, β, C)
@inline A_mul_B!(α::T, A::Symmetric{T}, B::Vector{T}, β::T, C::Vector{T}) where T<:Union{Float64, Complex{Float64}} = BLAS.symv!(A.uplo, α, A.data, B, β, C)

@inline Base.similar(VA::VectorOfArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = similar(VA, T)
@inline Base.similar(VA::VectorOfArray, dims::Base.Dims{N}) where {N} = similar(VA)

############# Exact Approximation
struct ExactRHSFunction{TS<:AbstractSparseMatrix} <: AbstractRHSFunction
    H::TS
end
function build_RHS_function(A::Approx) where {Approx<:ExactApprox}
    ExactRHSFunction(build_Hamiltonian(A))
end

function (EF::ExactRHSFunction)(t, state::T, dstate::T) where {T<:AbstractVector{Complex{Float64}}}
    A_mul_B!(-1im, EF.H, state, zero(Complex{Float64}), dstate)
end
function (EF::ExactRHSFunction)(t, u::T, du::T) where {T<:VectorOfArray{Complex{Float64},2,Array{T1,1}}} where T1<:AbstractVector{Complex{Float64}}
    @inbounds for i in eachindex(u)
        A_mul_B!(-1im, EF.H, u[i], zero(Complex{Float64}), du[i])
    end
end

############# Structured Approximation

struct ClusteredRHSFunction{TH<:AbstractSparseMatrix, TIM<:Union{AbstractMatrix,AbstractSparseMatrix}, T3<:AbstractSparseMatrix} <: AbstractRHSFunction
    Js::NTuple{3, Float64}
    H::TH
    edgeSpinOperators::Matrix{T3}###row - spin position, col - spin projection
    average_values::VectorOfArray{Float64, 2, Array{Vector{Float64},1}}
    #average_values_arr::VectorOfArray{Float64, 3, Array{Matrix{Float64},1}}
    cl_field::VectorOfArray{Float64, 2, Array{Vector{Float64},1}}
    #cl_field_arr::VectorOfArray{Float64, 3, Array{Matrix{Float64},1}}
    cluster_temp_state::Vector{Complex{Float64}}
    IM::TIM
    num_of_espins::Int64
end
function build_RHS_function(A::Approx) where {Approx<:ClusteredApprox}
    H = build_Hamiltonian(A)
    Dh = H.m
    edge_spins, IM = build_InterClusterInteractions(A)
    average_values = VectorOfArray([zeros(Float64, size(IM,1)), zeros(Float64, size(IM,1)), zeros(Float64, size(IM,1))])
    #average_values_arr = VectorOfArray([reshape(average_values[i], (length(edge_spins), A.cluster_num)) for i in eachindex(average_values)])
    cl_field = VectorOfArray([zeros(Float64, size(IM,1)), zeros(Float64, size(IM,1)), zeros(Float64, size(IM,1))])
    #cl_field_arr = VectorOfArray([reshape(cl_field[i], (length(edge_spins), A.cluster_num)) for i in eachindex(cl_field)])
    edgeSpinOperators = Matrix{SharedSparseMatrixCSC{Complex{Float64}, Int64}}(3,length(edge_spins))
    for i in eachindex(edge_spins) for sigma in 1:3
        edgeSpinOperators[sigma, i] = build_Spin_Operator(Dh, edge_spins[i], sigma)
    end end
    cluster_temp_state = Vector{Complex{Float64}}(Dh)
    ClusteredRHSFunction(A.L.Js, H, edgeSpinOperators, average_values, cl_field, cluster_temp_state, IM, length(edge_spins))
end

function (ClF::ClusteredRHSFunction)(t, u::T, du::T) where T<:VectorOfArray{Complex{Float64},2,Array{Array{Complex{Float64},1},1}}
    dims = (ClF.num_of_espins, length(u))
    @inbounds for cluster in eachindex(u)
        cl_norm = norm(u[cluster])
        for espin in 1:ClF.num_of_espins for sigma in 1:3
            A_mul_B!(ClF.cluster_temp_state, ClF.edgeSpinOperators[sigma, espin], u[cluster])
            ClF.average_values[sigma][sub2ind(dims, espin, cluster)] = real(dot(u[cluster], ClF.cluster_temp_state))/cl_norm
        end end
    end
    @inbounds for sigma in 1:3
        A_mul_B!(ClF.Js[sigma], ClF.IM, ClF.average_values[sigma], zero(Float64), ClF.cl_field[sigma])
    end
    @inbounds for cluster in eachindex(u)
        A_mul_B!(du[cluster], ClF.H, u[cluster])
        for espin in 1:ClF.num_of_espins for sigma in 1:3
            A_mul_B!(ClF.cl_field[sigma][sub2ind(dims, espin, cluster)], ClF.edgeSpinOperators[sigma, espin], u[cluster], one(Complex{Float64}), du[cluster])
        end end
        du[:,cluster] .*= -1im
    end
end

############# PureClassical Approximation

struct PureClassicalRHSFunction{T<:AbstractMatrix} <: AbstractRHSFunction
    Js::NTuple{3,Float64}
    hs::NTuple{3, Float64}
    IM::T
    cl_field::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
end
function build_RHS_function(A::Approx) where Approx<:PureClassicalApprox
    cl_field = deepcopy(randomState(A))
    IM = build_ClassicalInteractions(A)
    PureClassicalRHSFunction{typeof(IM)}(A.L.Js, A.L.hs, IM, cl_field)
end
function build_RHS_function(A::PureClassicalApprox, factor::Float64)
    cl_field = deepcopy(randomState(A))
    IM = build_ClassicalInteractions(A, factor)
    PureClassicalRHSFunction{typeof(IM)}(A.L.Js, A.L.hs, IM, cl_field)
end

function (PCF::PureClassicalRHSFunction)(t, u::T, du::T) where T<:VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
    @inbounds for sigma in 1:3
        fill!(PCF.cl_field[sigma], PCF.hs[sigma])
        A_mul_B!(PCF.Js[sigma], PCF.IM, u[sigma], one(Float64), PCF.cl_field[sigma])
    end
    @. du[:,1] = u[2]*PCF.cl_field[3] - u[3]*PCF.cl_field[2]
    @. du[:,2] = u[3]*PCF.cl_field[1] - u[1]*PCF.cl_field[3]
    @. du[:,3] = u[1]*PCF.cl_field[2] - u[2]*PCF.cl_field[1]
end

############# Hybrid Approximation

struct HybridRHSFunction{TH<:AbstractSparseMatrix, T1<:AbstractMatrix, T2<:AbstractMatrix, T3<:AbstractSparseMatrix} <: AbstractRHSFunction
    Js::NTuple{3,Float64}
    hs::NTuple{3, Float64}
    H::TH
    IM::T1
    q2cl::T2
    cl2q::T2
    average_values::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}
    edgeSpinOperators::Matrix{T3}
    cl_field::ArrayPartition{Float64, NTuple{2,VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}
    q_temp_state::Vector{Complex{Float64}}
end
function build_RHS_function(A::HybridApprox)
    H = build_Hamiltonian(A)
    Dh = H.m
    IM = build_ClassicalInteractions(A)
    edge_spins, q2cl, cl2q = build_QuantumClassicalInteractions(A)
    average_values = random_cl_state(length(edge_spins))
    cl_field = ArrayPartition(random_cl_state(length(edge_spins)), random_cl_state(A.cl_num))
    edgeSpinOperators = Matrix{SharedSparseMatrixCSC{Complex{Float64}, Int64}}(length(edge_spins),3)
    for espin in eachindex(edge_spins)
        for sigma in 1:3
            edgeSpinOperators[espin, sigma] = build_Spin_Operator(Dh, edge_spins[espin], sigma)
        end
    end
    q_temp_state = Vector{Complex{Float64}}(Dh)
    HybridRHSFunction(A.L.Js, A.L.hs, H, IM, q2cl, cl2q, average_values, edgeSpinOperators, cl_field, q_temp_state)
end

function (HF::HybridRHSFunction)(t, u::T, du::T) where T<:ArrayPartition{Complex{Float64}, Tuple{Vector{Complex{Float64}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}
    qnorm = norm(u.x[1])
    @inbounds for sigma in 1:3
        for q in 1:size(HF.edgeSpinOperators,1)
            A_mul_B!(HF.q_temp_state, HF.edgeSpinOperators[q, sigma], u.x[1])
            HF.average_values[q, sigma] = real(dot(HF.q_temp_state, u.x[1]))/qnorm
        end
    end
    A_mul_B!(du.x[1], HF.H, u.x[1])
    @inbounds for sigma in 1:3
        A_mul_B!(HF.Js[sigma], HF.cl2q, u.x[2][sigma], zero(Float64), HF.cl_field.x[1][sigma])
        for q in 1:size(HF.edgeSpinOperators, 1)
            A_mul_B!(HF.cl_field.x[1][q,sigma], HF.edgeSpinOperators[q, sigma], u.x[1], one(Complex{Float64}), du.x[1])
        end
        fill!(HF.cl_field.x[2][sigma], HF.hs[sigma])
        A_mul_B!(HF.Js[sigma], HF.IM, u.x[2][sigma],  one(Float64), HF.cl_field.x[2][sigma])
        A_mul_B!(HF.Js[sigma], HF.q2cl, HF.average_values[sigma],  one(Float64), HF.cl_field.x[2][sigma])
    end
    du.x[1] .*= -1im
    @. du.x[2][:,1] = u.x[2][2]*HF.cl_field.x[2][3] - u.x[2][3]*HF.cl_field.x[2][2]
    @. du.x[2][:,2] = u.x[2][3]*HF.cl_field.x[2][1] - u.x[2][1]*HF.cl_field.x[2][3]
    @. du.x[2][:,3] = u.x[2][1]*HF.cl_field.x[2][2] - u.x[2][2]*HF.cl_field.x[2][1]
end

############# Composite Approximation

struct CompositeRHSFunction{RHSFunction1, RHSFunction2} <: AbstractRHSFunction
    RHS1::RHSFunction1
    RHS2::RHSFunction2
    q1tocl2::SharedMatrix{Float64}
    cl2toq1::SharedMatrix{Float64}
    cl1tocl2::SharedMatrix{Float64}
    cl2tocl1::SharedMatrix{Float64}
end

function build_RHS_function(A::CompositeApproximation{Approx1,Approx2}) where {Approx1<:HybridApprox, Approx2<:PureClassicalApprox}
    RHS1 = build_RHS_function(A.A1)
    RHS2 = build_RHS_function(A.A2, A.gamma^2)
    q1tocl2, cl2toq1, cl1tocl2, cl2tocl1 = build_InterLattice_Couplings(A, size(RHS1.edgeSpinOperators)[1])
    CompositeRHSFunction(RHS1,RHS2, q1tocl2, cl2toq1, cl1tocl2, cl2tocl1)
end

function (CF::CompositeRHSFunction)(t, u::T, du::T) where T<:ArrayPartition{Complex{Float64},Tuple{Array{Complex{Float64},1},VectorOfArray{Float64,2,Array{Array{Float64,1},1}},VectorOfArray{Float64,2,Array{Array{Float64,1},1}}}}
    ############## Unpacking
    HF = CF.RHS1
    PCF = CF.RHS2
    ############## Quantum averages
    qnorm = norm(u.x[1])
    @inbounds for sigma in 1:3
        for q in 1:size(HF.edgeSpinOperators,1)
            A_mul_B!(HF.q_temp_state, HF.edgeSpinOperators[q, sigma], u.x[1])
            HF.average_values[q, sigma] = real(dot(HF.q_temp_state, u.x[1]))/qnorm
        end
    end
    ############## Calculation of effective magnetic fields
    @inbounds for sigma in 1:3
        A_mul_B!(HF.Js[sigma], HF.cl2q, u.x[2][sigma], zero(Float64), HF.cl_field.x[1][sigma]) # quantum part
        fill!(HF.cl_field.x[2][sigma], HF.hs[sigma]) # init L1 fields
        fill!(PCF.cl_field[sigma], PCF.hs[sigma]) # init L2 fields
        A_mul_B!(HF.Js[sigma], HF.IM, u.x[2][sigma], one(Float64), HF.cl_field.x[2][sigma]) # cl2cl L1 fields
        A_mul_B!(PCF.Js[sigma], PCF.IM, u.x[3][sigma], one(Float64), PCF.cl_field[sigma]) # cl2cl L2 fields
        A_mul_B!(HF.Js[sigma], HF.q2cl, HF.average_values[sigma], one(Float64), HF.cl_field.x[2][sigma]) #cl21 L1 fields#
    end
    #### Cross fields
    A_mul_B!(PCF.Js[3], CF.cl2toq1, u.x[3][3], one(Float64), HF.cl_field.x[1][3])
    A_mul_B!(PCF.Js[3], CF.cl2tocl1, u.x[3][3], one(Float64), HF.cl_field.x[2][3])
    A_mul_B!(HF.Js[3], CF.q1tocl2, HF.average_values[3], one(Float64), PCF.cl_field[3])
    A_mul_B!(HF.Js[3], CF.cl1tocl2, u.x[2][3], one(Float64), PCF.cl_field[3])
    ############## Apply fields
    A_mul_B!(du.x[1], HF.H, u.x[1])
    @inbounds for sigma in 1:3
        for q in 1:size(HF.edgeSpinOperators, 1)
            A_mul_B!(HF.cl_field.x[1][q,sigma], HF.edgeSpinOperators[q, sigma], u.x[1], one(Complex{Float64}), du.x[1])
        end
    end
    du.x[1] .*= -1im
    @. du.x[2][:,1] = u.x[2][2]*HF.cl_field.x[2][3] - u.x[2][3]*HF.cl_field.x[2][2]
    @. du.x[2][:,2] = u.x[2][3]*HF.cl_field.x[2][1] - u.x[2][1]*HF.cl_field.x[2][3]
    @. du.x[2][:,3] = u.x[2][1]*HF.cl_field.x[2][2] - u.x[2][2]*HF.cl_field.x[2][1]
    @. du.x[3][:,1] = u.x[3][2]*PCF.cl_field[3] - u.x[3][3]*PCF.cl_field[2]
    @. du.x[3][:,2] = u.x[3][3]*PCF.cl_field[1] - u.x[3][1]*PCF.cl_field[3]
    @. du.x[3][:,3] = u.x[3][1]*PCF.cl_field[2] - u.x[3][2]*PCF.cl_field[1]
end


text:
/* Test boolean operators */

int main () {
  int x = 4;
  if (3 <= x && 4 != 2 && true) {
    printBool(true);
  } else {
    printString("apa");
  }

  printBool(true == true || dontCallMe(1));
  printBool(4.0 < -50.0 && dontCallMe(2));

  printBool(4 == x && true == !false && true);

  printBool(implies(false,false));
  printBool(implies(false,true));
  printBool(implies(true,false));
  printBool(implies(true,true));
  return 0 ;

}

boolean dontCallMe(int x) {
  printInt(x);
  return true;
}

void printBool(boolean b) {
  if (b) {
    printString("true");
  } else {
    printString("false");
 }
 return;
}

boolean implies(boolean x, boolean y) {
  return !x || x == y;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
module hhNeuron

mutable struct Neuron
    Vref::Float64
    gNa::Float64
    gK::Float64
    gL::Float64
    ENa::Float64
    EK::Float64
    EL::Float64
    Cm::Float64
    v::Float64
    m::Float64
    h::Float64
    n::Float64
end

function advance!(neuron, dt, I)
    # calculate the ion channel currents
    INa = neuron.gNa * neuron.m*neuron.m*neuron.m * neuron.h * (neuron.ENa - neuron.v)
    IK = neuron.gK * neuron.n*neuron.n*neuron.n*neuron.n * (neuron.EK - neuron.v)
    IL = neuron.gL * (neuron.EL - neuron.v)

    # calculate channel kinetics
    aM = 0.1 * (neuron.v - neuron.Vref -25) / ( 1 - exp(-(neuron.v-neuron.Vref-25)/10))
    bM = 4 * exp(-(neuron.v-neuron.Vref)/18)

    aH = 0.07 * exp(-(neuron.v-neuron.Vref)/20)
    bH = 1 / ( 1 + exp( -(neuron.v-neuron.Vref-30)/10 ) )

    aN = 0.01 * (neuron.v-neuron.Vref-10) / ( 1 - exp(-(neuron.v-neuron.Vref-10)/10) )
    bN = 0.125 * exp(-(neuron.v-neuron.Vref)/80)

    # calculate derivatives
    dv_dt = ( INa + IK + IL + I*1e-6 ) / neuron.Cm
    dm_dt = (1.0-neuron.m) * aM - neuron.m * bM
    dh_dt = (1.0-neuron.h) * aH - neuron.h * bH
    dn_dt = (1.0-neuron.n) * aN - neuron.n * bN

    # calculate next step
    neuron.v += dv_dt * dt
    neuron.m += dm_dt * dt
    neuron.h += dh_dt * dt
    neuron.n += dn_dt * dt

    return(neuron)
end

function simulation(dt::Float64=0.01, T::Int=1000)

    const numSteps = T/dt

    # initial membrane potential
    const Vinit = -70.0
    const Vref = -70.0

    # neuron morphology and membrane conductance
    const Smemb = 4000.0                                # surface area
    const Cmemb = 1.0                                   # membrane conductance
    const Cm = Cmemb * Smemb * 1e-8

    # channel parameter
    const gNa = 120.0 * Cmemb * Smemb * 1e-8            # sodium conductance
    const gK = 36.0 * Cmemb * Smemb * 1e-8              # potassium conductance
    const gL = 0.3 * Cmemb * Smemb * 1e-8               # leack conductance
    const ENa = 125.0                                   # sodium equilibrium potential
    const EK = -70.0                                    # potassium equilibrium potential
    const EL = -25.0                                    # leak equilibrium potential

    # initial values for m, h, n
    m = 0.0529
    h = 0.6
    n = 0.32

    v = []
    neuron = Neuron( Vref, gNa, gK, gL, ENa, EK, EL, Cm, Vinit, m, h, n )

    for i = 1:numSteps
        advance!(neuron, dt, 1000)
        push!(v, neuron.v)
    end

    return(v)
end

end

max_stars_repo_name: Trende/julia-neuron-models

max_issues_repo_name: Trende/julia-neuron-models

max_forks_repo_name: Trende/julia-neuron-models


text:
using Core: CodeInfo, Typeof
using Core.Compiler: InferenceState, widenconst, svec
using InteractiveUtils: typesof

worldcounter() = ccall(:jl_get_world_counter, UInt, ())

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

struct TypedMeta
  frame::InferenceState
  method::Method
  code::CodeInfo
  ret
end

function Base.show(io::IO, meta::TypedMeta)
  print(io, "Typed metadata for ")
  print(io, meta.method)
end

define_typeinf_code2() = isprecompiling() ||
@eval Core.Compiler function typeinf_code2(method::Method, @nospecialize(atypes), sparams::SimpleVector, run_optimizer::Bool, params::Params)
    if $(VERSION >= v"1.2-")
      code = specialize_method(method, atypes, sparams)
    else
      code = code_for_method(method, atypes, sparams, params.world)
    end
    code === nothing && return (nothing, Any)
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = InferenceResult(code)
    frame = InferenceState(result, false, params)
    frame === nothing && return (nothing, Any)
    if typeinf(frame) && run_optimizer
        opt = OptimizationState(frame)
        optimize(opt, result.result)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return frame
end

"""
    typed_meta(Tuple{...})

Same as [`@meta`](@ref), but represents the method after type inference. IR
constructed with typed metadata will have type annotations.

See also [`@typed_meta`](@ref).

    julia> IRTools.typed_meta(Tuple{typeof(gcd),Int,Int})
    Typed metadata for gcd(a::T, b::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at intfuncs.jl:31
"""
function typed_meta(T; world = worldcounter(), optimize = false)
  F = T.parameters[1]
  F isa DataType && (F.name.module === Core.Compiler ||
                     F <: Core.Builtin ||
                     F <: Core.Builtin) && return nothing
  _methods = Base._methods_by_ftype(T, -1, world)
  length(_methods) == 1 || return nothing
  type_signature, sps, method = first(_methods)
  params = Core.Compiler.Params(world)
  frame = Core.Compiler.typeinf_code2(method, type_signature, sps, optimize, params)
  ci = frame.src
  ci.inferred = true
  ci.method_for_inference_limit_heuristics = method
  if ci.ssavaluetypes == 0 # constant return; IRCode doesn't like this
    ci.ssavaluetypes = Any[Any]
  end
  return TypedMeta(frame, method, ci, widenconst(frame.result.result))
end

struct Meta
  method::Method
  code::CodeInfo
  nargs::Int
  sparams
end

function Base.show(io::IO, meta::Meta)
  print(io, "Metadata for ")
  print(io, meta.method)
end

# Workaround for what appears to be a Base bug
untvar(t::TypeVar) = t.ub
untvar(x) = x

"""
    meta(Tuple{...})

Construct metadata for a given method signature. Metadata can then be used to
construct [`IR`](@ref) or used to perform other reflection on the method.

See also [`@meta`](@ref), [`typed_meta`](@ref).

    julia> IRTools.meta(Tuple{typeof(gcd),Int,Int})
    Metadata for gcd(a::T, b::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at intfuncs.jl:31
"""
function meta(T; world = worldcounter())
  F = T.parameters[1]
  F == typeof(invoke) && return invoke_meta(T; world = world)
  F isa DataType && (F.name.module === Core.Compiler ||
                     F <: Core.Builtin ||
                     F <: Core.Builtin) && return nothing
  _methods = Base._methods_by_ftype(T, -1, world)
  length(_methods) == 0 && return nothing
  type_signature, sps, method = last(_methods)
  sps = svec(map(untvar, sps)...)
  @static if VERSION >= v"1.2-"
    mi = Core.Compiler.specialize_method(method, type_signature, sps)
    ci = Base.isgenerated(mi) ? Core.Compiler.get_staged(mi) : Base.uncompressed_ast(method)
  else
    mi = Core.Compiler.code_for_method(method, type_signature, sps, world, false)
    ci = Base.isgenerated(mi) ? Core.Compiler.get_staged(mi) : Base.uncompressed_ast(mi)
  end
  ci.method_for_inference_limit_heuristics = method
  if isdefined(ci, :edges)
    ci.edges = Core.MethodInstance[mi]
  end
  Base.Meta.partially_inline!(ci.code, [], method.sig, Any[sps...], 0, 0, :propagate)
  Meta(method, ci, method.nargs, sps)
end

function invoke_tweaks!(ci::CodeInfo)
  ci.slotnames = [:invoke, ci.slotnames[1], :T, ci.slotnames[2:end]...]
  ci.slotflags = [0x00, ci.slotflags[1], 0x00, ci.slotflags[2:end]...]
  ci.code = map(ci.code) do x
    prewalk(x) do x
      x isa SlotNumber ? SlotNumber(x.id == 1 ? 2 : x.id+2) : x
    end
  end
end

function invoke_meta(T; world)
  F = T.parameters[2]
  A = T.parameters[3]::Type{<:Type{<:Tuple}}
  T = Tuple{F,A.parameters[1].parameters...}
  m = meta(T, world = world)
  invoke_tweaks!(m.code)
  return Meta(m.method, m.code, m.nargs+2, m.sparams)
end

"""
    @meta f(args...)

Convenience macro for retrieving metadata without writing a full type signature.

    julia> IRTools.@meta gcd(10, 5)
    Metadata for gcd(a::T, b::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at intfuncs.jl:31
"""
macro meta(ex)
  isexpr(ex, :call) || error("@meta f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(meta(typesof($(esc.((f, args...))...))))
end

"""
    @typed_meta f(args...)

Convenience macro for retrieving typed metadata without writing a full type signature.

    julia> IRTools.@typed_meta gcd(10, 5)
    Typed metadata for gcd(a::T, b::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at intfuncs.jl:31
"""
macro typed_meta(ex)
  isexpr(ex, :call) || error("@meta f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(typed_meta(typesof($(esc.((f, args...))...))))
end

function code_ir(f, T)
  m = meta(Tuple{Typeof(f),T.parameters...})
  return IR(m)
end

function code_irm(ex)
  isexpr(ex, :call) || error("@code_ir f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(code_ir($(esc(f)), typesof($(esc.(args)...))))
end

"""
    @code_ir f(args...)

Convenience macro similar to `@code_lowered` or `@code_typed`. Retrieves the IR
for the given function call.

    julia> @code_ir gcd(10, 5)
    1: (%1, %2, %3)
      %4 = %2 == 0
      br 4 unless %4
    2: ...
"""
macro code_ir(ex)
  code_irm(ex)
end

codeinfo(m::Meta) = m.code
codeinfo(m::TypedMeta) = m.code

function argnames!(meta, names...)
  meta.code.slotnames = [names...]
end


text:
function seager_mr(mp,comp)
# computes the radius of a planet
# given the mass & composition of the planet
# according to the formulae in Seager et al. (2007)
k1=-0.20945
k2=0.0804
k3=0.394
if comp == "Fe"
  m1=4.34
  r1=2.23
end
if comp == "Earth"
  m1=6.41
  r1=3.19
#  r1=2.84
end
if comp == "MgSiO3"
  m1=7.38
  r1=3.58
end
if comp == "H2O"
  m1=8.16
  r1=4.73
end
if comp == "FeMgH2O"
#  m1 = 6.88
#  r1 = 4.02
  m1 = 6.41
  r1 = 3.63
end
ms=mp/m1
lgrp=k1 .+ log10.(ms)/3.0 .-k2*ms.^k3
return 10.0.^lgrp*r1
end


text:
module TestInterpolations
using Test
using AxisArrayHelpers
const AH = AxisArrayHelpers
import AxisKeys
const AK = AxisKeys
using AxisArrayConversion
const AC = AxisArrayConversion
using CoordinateTransformations

@testset "create_interpolate" begin
    data = (axes=(1:3,), values=[10,11,13])
    itp = @inferred create_interpolate(data)
    @test itp(1  ) ≈ 10
    @test itp(1.2) ≈ 10.2
    @test itp(2.2) ≈ 11.4
    @test_throws ArgumentError itp(0.9)

    datas = [
        AH.to(AK.KeyedArray, data),
        AH.to(Array, data),
        AH.to(NamedTuple, data),
    ]

    for data in datas
        itp = @inferred create_interpolate(data)
        @test itp(1  ) ≈ 10
        @test itp(1.2) ≈ 10.2
        @test itp(2.2) ≈ 11.4
        @test_throws ArgumentError itp(0.9)
    end

    data = (axes=(-2:2.0, ), values=[0,1,0,0,10.0])
    itp = @inferred create_interpolate(data, extrapolate=:replicate)
    @test itp(-2) === 0.0
    @test itp(-100) === 0.0
    @test itp(100) === 10.0
end

@testset "pullback" begin
    data = (axes=(1:3,), values=[2,0,10])
    pb = @inferred pullback(identity, (x=1:3,), data)
    AC.check_consistency(pb)
    @test pb.axes === (x=1:3,)
    @test pb.values == data.values

    data = (axes=(a=1:3,), values=[2,0,10.0])
    out = deepcopy(data)
    out.values .= NaN
    pb = @inferred pullback!(identity, out, data)
    @test pb === out
    @test out == data

    data = (axes=(a=1:3,), values=[2,0,10])
    pb_axes = (a=(1/2)*(1:3),)
    #pb_axes = (a=[2,2],)
    pb = @inferred pullback(pt -> 2pt, pb_axes, data, extrapolate=:replicate)
    AC.check_consistency(pb)
    @test pb.axes === pb_axes
    @test pb.values == data.values
end

@testset "pullback_axes" begin
    #A = [0 5 0; 1 0 0; 0 0 -1]
    #axes = (1:10, 2:5, 4:8)
    #N = length(axes)

    #sp = splitperm(A)
    #A2 = matrix_from_splitperm(sp)
    #@test A ≈ A2
    #Ainv = matrix_from_splitperm(inv_splitperm(sp))
    #@test A * Ainv ≈ LinearAlgebra.I

    #function permute(v, perm)
    #    permute!(copy(v), collect(perm))
    #end

    #v = randn(3)
    using Test
    @test (0:9, 0:18)      === AH.pullback_axes(Translation([1,2]), (x=1:10, y=2:20))
    @test_throws ArgumentError AH.pullback_axes(Translation([1,2,3]), (x=1:10, y=2:20))

    @test (2:2:20.0, 1:10.0)  === AH.pullback_axes(LinearMap([0 1; 1 0]), (1:10, 2:2:20))
    @test (0.5:0.5:5, 1:5.0)  === AH.pullback_axes(LinearMap([2 0; 0 1]), (1:10, 1:5))
    @test (-5:1:-1.0, 1:1:5.0) == AH.pullback_axes(LinearMap([-1 0; 0 1]), (1:5, 1:5))

    for lin in [
            LinearMap([1 0 0; 0 1 0; 0 0 1]),
            LinearMap([randn() 0 0; 0 randn() 0; 0 0 randn()]),
            LinearMap([0 1 0; -1 0 0; 0 0 2]),
            LinearMap([0 1 0; 0 0 1; 1 0 0]),
            LinearMap([0 randn() 0; 0 0 randn(); randn() 0 0]),
        ]
        trans = Translation(randn(3))
        trafo = trans ∘ lin
        axs = (1:5, -2:2:4, 1:10)
        axs_pb  = AH.pullback_axes(trafo, axs)
        axs_pb1 = AH.pullback_axes(trans, axs)
        axs_pb2 = AH.pullback_axes(lin, axs_pb1)
        @test all(length.(axs_pb2) .== length.(axs_pb))
        @test all(axs_pb2 .≈ axs_pb)
    end
end

@testset "common_grid" begin
    grid1 = (1:10,)
    grid2 = (2:20,)
    @test AH.common_grid([grid1, grid2]) == (2:10,)
    grid1 = (1:10,-10:2:10)
    grid2 = (2:20,-8:2:4)
    @test AH.common_grid([grid1, grid2]) == (2:10,-8:2:4)
end

@testset "common_ground" begin
    data1 = (axes=(a=1:10,), values=1:10)
    data2 = (axes=(a=0:4,), values=0:10:40)
    out1, out2 = AH.common_ground((data1, data2))
    @test out1 == (axes = (a = [1, 2, 3, 4],), values = [1.0, 2.0, 3.0, 4.0])
    @test out2 == (axes = (a = [1, 2, 3, 4],), values = [10.0, 20.0, 30.0, 40.0])
end

@testset "flip_decreasing_axes" begin
    @test AH.flip_decreasing_axes((axes=(1:3,), values=[1,4,2])) == (axes=(1:3,), values=[1,4,2])
    @test AH.flip_decreasing_axes((axes=([0,-2, -4],), values=[1,4,2])) == (axes=([-4, -2, 0],), values=[2,4,1])
    nt1 = (axes=(x=[-1, -2], y=[1,2,3]), values = [1 2 3; 4 5 6])
    nt2 = (axes = (x = [-2, -1], y = [1, 2, 3]), values = [4 5 6; 1 2 3])
    @test AH.flip_decreasing_axes(nt1) == nt2
end
end#module


text:
include("Technology.jl")

mutable struct SystemState <: ScenarioNode
    year::Int64
    capacityMix::Array{Technology,1}
    installedCapacity::Float64
    installedLimitByTechnologyType::Dict{Symbol,Real} #MW
    installedByTechnologyType::Dict{Symbol,Real} #MW
    emissionPrice::Real

    CalculateInstalledCapacity::Function
    TimeSeries::Function

    function SystemState(year,capacityMix)
        self = new()
        self.year = year
        self.capacityMix = capacityMix
        self.installedCapacity = 0

        self.TimeSeries = function(day::Int64, hour::Int64, systemLoad::Load, yearlyTechnologyMix::Array{Technology,1})
            systemLoad.TimeSeries(day,hour)
            for technology in yearlyTechnologyMix
                if typeof(technology) == RenewablePlants
                    technology.TimeSeries(day,hour)
                end
            end
        end

        self.CalculateInstalledCapacity = function()
            for technology in self.capacityMix
                self.installedCapacity += technology.designProperties.installedCapacity
            end
        end

        return self
    end
end

mutable struct SystemHistory <: Buffer
    systemCapacityMix::Array{Array{Technology,1},1}
    systemElectricityPrice::Array{Float64,1}
    systemDecisionSequence::Array{Array{Array{Float64,1},1},1}
    capacityMixProjection::Array{Array,1}
    InvestmentHistory::Array{Symbol,1}

    function SystemHistory()
        self = new()
        self.systemCapacityMix = Array{Array{Technology,1},1}(undef,0)
        self.systemElectricityPrice = Array{Float64,1}(undef,0)
        self.systemDecisionSequence = Array{Array{Array{Float64,1},1},1}(undef,0)
        self.capacityMixProjection = Array{Array,1}(undef,0)
        self.InvestmentHistory = Array{Symbol,1}(undef,0)

        return self
    end

end


text:
using Interpolations
using Plots
# Given CDF generate independent random numbers
function rCDF(x,y,n)
    ru = rand(n)
    a = LinearInterpolation(y, x, extrapolation_bc = Flat())
    xout = a(ru)
    return xout
end
# Aux function to be used for mountain plot
function breakCDF(x,y)
    xp =  x .>= 0
    xn =  x .< 0
    yp =  1 .-y[xp]
    yn =  y[xn]
    return x[xp], yp, x[xn] ,yn
end
# Giving vector of random numbers return CDF
function getCDF(x)
    sx = sort(x)
    n  = length(x)
    y  = (1:n)/(n+1)
    return sx,y
end
# CDF plot
function mountainPlot(x)
    x,y = getCDF(x)
    xp,yp,xn,yn = breakCDF(x,y)
    plot(xp,yp,color="blue",yaxis=:log);plot!(xn,yn,color="blue",yaxis=:log)
end

function getrnd(r,n)
    x,y = getCDF(r)
    xout = rCDF(x,y,n)
    return xout
end

max_stars_repo_name: silvaac/SMCarticle

max_issues_repo_name: silvaac/SMCarticle

max_forks_repo_name: silvaac/SMCarticle


text:
# Meachnical problems

# force XX.5 -> XX + 1
round5(x) = floor(Int, x + 0.5)

# 1
function pressure_vessel(x)
    ## Pressure Vessal Problems
    x1 = round5(x[1]);
    x2 = round5(x[2]);
    x3 = x[3];
    x4 = x[4];
    z1 = 0.0625*x1;
    z2 = 0.0625*x2;
    ## objective function
    f = zeros(2);
    f[1] = 1.7781 .* z1 .* x3 .^ 2+0.6224 .* z1 .* x2 .* x4+3.1661 .* z1 .^ 2 .* x4+19.84 .* z1 .^ 2 .* x3;
    f[2] = -pi .* x3 .^ 2 .* x4-(4/3) .* pi .* x3 .^ 3;
    ## constraints
    g = zeros(2);
    h = zeros(1);
    g[1] = 0.00954 .* x3-z2;
    g[2] = 0.0193 .* x3-z1;

    return f, g, h
end

# 2
function vibrating_platform(x)
    ## Vibrating Platform
    d1 = x[1];
    d2 = x[2];
    d3 = x[3];
    b = x[4];
    L = x[5];
    rho1 = 100; rho2 = 2770; rho3 = 7780;
    E1 = 1.6; E2 = 70; E3 = 200;
    c1 = 500; c2 = 1500; c3 = 800;
    mu = 2*b .* (rho1 .* d1+rho2 .* (d2-d1)+rho3 .* (d3-d2));
    EI = (2*b ./ 3) .* (E1 .* d1 .^ 3+E2 .* (d2 .^ 3-d1 .^ 3)+rho3 .* (d3-d2));
    ## objective function
    f = zeros(2);
    f[1] = (-pi) ./ (2*L) .^ 2 .* (abs(EI ./ mu)) .^ 0.5;
    f[2] = 2*b .* L .* (c1 .* d1+c2 .* (d2-d1)+c3 .* (d3-d2));
    ## constraints
    g = zeros(5);
    h = zeros(1);
    g[1] = mu .* L -2800;
    g[2] = d1-d2;
    g[3] = d2-d1-0.15;
    g[4] = d2-d3;
    g[5] = d3-d2-0.01;

    return f, g, h
end

# 3
function two_bar_Truss_design_problems(x)
    ## Two bar Truss design Problems
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    ## objective function
    f = zeros(2);
    f[1] = x1 .* (16+x3 .^ 2) .^ (0.5)+x2 .* (1+x3 .^ 2) .^ (0.5);
    f[2] = (20 .* (16+x3 .^ 2) .^ (0.5)) ./ (x3 .* x1);
    ## constraints
    g = zeros(3);
    h = zeros(1);
    g[1] = f[1]-0.1;
    g[2] = f[2]-1e5;
    g[3] = (80 .* (1+x3 .^ 2) .^ (0.5)) ./ (x3 .* x2)-1e5;

    return f, g, h
end

# 4
function weldan_beam_design(x)
    ## Weldan Beam Design Problem
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    P = 6000;
    L = 14;
    E = 30e6;
    tmax = 13600;
    sigmax = 30000;
    G = 12e6;
    Pc = (4.013 .* E .* ((x3 .^ 2+x4 .^ 6) ./ 36) .^ 0.5) ./ (L .^ 2) .* (1-x3 ./ (2 .* L) .* (E ./ (4*G))^(0.5));
    sigma = (6 .* P .* L) ./ (x4 .* x3 .^ 2);
    J = 2 .* (sqrt(2) .* x1 .* x2 .* (x2 .^ 2/12+((x1+x3) ./ 2) .^ 2));
    R = sqrt(x2 .^ 2/4+((x1+x3) ./ 2) .^ 2);
    M = P .* (L+x2 ./ 2);
    tho1 = P ./ (sqrt(2) .* x1 .* x2);
    tho2 = M .* R ./ J;
    tho = sqrt(tho1 .^ 2+2 .* tho1 .* tho2 .* x2 ./ (2*R)+tho2 .^ 2);
    ## objective function
    f = zeros(2);
    f[1] = 1.10471 .* x1 .^ 2 .* x2+0.04811 .* x3 .* x4 .* (14+x2);
    f[2] = (4 .* P .* L .^ 3) ./ (E .* x4 .* x3 .^ 3);
    ## constraints
    g = zeros(4);
    h = zeros(1);
    g[1] = tho -tmax;
    g[2] = sigma - sigmax;
    g[3] = x1-x4;
    g[4] = P - Pc;

    return f, g, h
end

# 5
function disc_brake_design(x)
    ## Disc Brake Design Problem
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    ## objective function
    f = zeros(2);
    f[1] = 4.9e-5 .* (x2 .^ 2-x1 .^ 2) .* (x4-1);
    f[2] = 9.82e6 .* (x2 .^ 2-x1 .^ 2) ./ (x3 .* x4 .* (x2 .^ 3-x1 .^ 3));
    ## constraints
    g = zeros(4);
    h = zeros(1);
    g[1] = 20 - (x2 - x1);
    g[2] = x3 ./ (3.14 .* (x2 .^ 2-x1 .^ 2)) - 0.4;
    g[3] = 2.22e-3 .* x3 .* (x2 .^ 3-x1 .^ 3) ./ (x2 .^ 2-x1 .^ 2) .^ 2-1;
    g[4] = 900 - 2.66e-2 .* x3 .* x4 .* (x2 .^ 3-x1 .^ 3) ./ (x2 .^ 2-x1 .^ 2);

    return f, g, h
end

# 6
function speed_reducer_design(x)
    ## Speed Reducer Design Problem
    x1 = x[1];
    x2 = x[2];
    x3 = round5(x[3]);
    x4 = x[4];
    x5 = x[5];
    x6 = x[6];
    x7 = x[7];
    ## objective function
    f = zeros(2);
    f[1] = 0.7854 * x1 * x2 ^ 2 * (10 *(x3 ^ 2)/3+14.933 * x3-43.0934) - 1.508 * x1 * (x6 ^ 2+x7 ^ 2)+7.477 * (x6 ^ 3+x7 ^ 3)+0.7854 * (x4 * x6 ^ 2+x5 * x7 ^ 2);
    f[2] = sqrt((745 .* x4 ./ (x2 .* x3)) .^ 2+1.69e7) ./ (0.1 .* x6 .^ 3);
    ## constraints
    g = zeros(11);
    h = zeros(1);
    g[1] = 1 ./ (x1 .* x2 .^ 2 .* x3)-1/27;
    g[2] = 1 ./ (x1 .* x2 .^ 2 .* x3 .^ 2)-1/397.5;
    g[3] = x4 .^ 3 ./ (x2 .* x3 .* x6 .^ 4)-1/1.93;
    g[4] = x5 .^ 3 ./ (x2 .* x3 .* x7 .^ 4)-1/1.93;
    g[5] = x2 .* x3-40;
    g[6] = x1 ./ x2-12;
    g[7] = -x1 ./ x2+5;
    g[8] = 1.9-x4+1.5 .* x6;
    g[9] = 1.9-x5+1.1 .* x7;
    g[10] = f[2]-1300;
    g[11] = sqrt((745 .* x5 ./ (x2 .* x3)) .^ 2+1.575e8) ./ (0.1*x7 .^ 3)-850;

    return f, g, h
end

# 7
function gear_train_design(x)
    ## Gear Train Design Problems
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    ## objective function
    f = zeros(2);
    f[1] = abs(6.931-x3 .* x4 ./ (x1 .* x2));
    f[2] = maximum(x);
    ## constraints
    g = zeros(1);
    h = zeros(1);
    g[1] = f[1] ./ 6.931-0.5;

    return f, g, h
end

# 8
function car_side_impact_design(x)
    ## Car Side Impact Desifn Problem
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    x5 = x[5];
    x6 = x[6];
    x7 = x[7];
    VMBP = 10.58-0.674 .* x1 .* x2-0.67275 .* x2;
    VFD  = 16.45-0.489 .* x3 .* x7-0.843 .* x5 .* x6;
    ## objective function
    f = zeros(3);
    f[1] = 1.98+4.9 .* x1 .* 6.67 .* x2+6.98 .* x3+4.01 .* x4+1.78 .* x5+1e-5 .* x6+2.73 .* x7;
    f[2] = 4.72-0.5 .* x4-0.19 .* x2 .* x3;
    f[3] = 0.5 .* (VMBP+VFD);
    ## constraints
    g = zeros(9);
    h = zeros(1);
    g[1] = -1+1.16-0.3717 .* x2 .* x4-0.0092928 .* x3;
    g[2] = -0.32+0.261-0.0159 .* x1 .* x2-0.06486 .* x1-0.019 .* x2 .* x7+0.0144 .* x2 .* x5+0.0154464 .* x6;
    g[3] = -0.32+0.74-0.61 .* x2-0.031296 .* x3-0.031872 .* x7+0.227 .* x2 .^ 2;
    g[4] = -0.32+0.214+0.00817 .* x5-0.045195 .* x1-0.0135168 .* x1+0.03099 .* x2 .* x6-0.018 .* x2 .* x7+0.007176 .* x3+0.023232 .* x3-0.00364 .* x5 .* x6-0.018 .* x2 .^ 2;
    g[5] = -32+33.86+2.95 .* x3-5.057 .* x1 .* x2-3.795 .* x2-3.4431 .* x7+1.45728;
    g[6] = -32+28.98+3.818 .* x3-4.2 .* x1 .* x2+1.27296 .* x6-2.68065 .* x7;
    g[7] = -32+46.36-9.9 .* x2-4.4505 .* x1;
    g[8] = f[2]-4;
    g[9] = VMBP - 9.9;

    return f, g, h
end

# 9
function four_bar_plane_truss(x)
    ## Four Bar Plane Truss
    F = 10; E = 2e5; L = 200; sig = 10;
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    ## objective function
    f = zeros(2);
    f[1] = L .* (2*x1+sqrt(2) .* x2+sqrt(2) .* x3+x4);
    f[2] = F .* L ./ E .* (2 ./ x1+2 .* sqrt(2) ./ x2-2 .* sqrt(2) ./ x3+2 ./ x4);
    ## constraints
    g = zeros(1);
    h = zeros(1);
    return f, g, h
end

# 10
function two_bar_plane_truss(x)
    ## Two bar plane Truss
    x1 = x[1];
    x2 = x[2];
    rho = 0.283;
    h = 100;
    P = 104;
    E = 3e7;
    rho0 = 2e4;
    ## objective function
    f = zeros(2);
    f[1] = 2 .* rho .* h .* x2 .* sqrt(1+x1 .^ 2);
    f[2] = rho .* h .* (1+x1 .^ 2) .^ 1.5 .* (1+x1 .^ 4) .^ 0.5 ./ (2 .* sqrt(2) .* E .* x1 .^ 2 .* x2);
    ## constraints
    g = zeros(2);
    h = zeros(2);
    g[1] = P .* (1+x1) .* (1+x1 .^ 2) .^ 0.5 ./ (2 .* sqrt(2) .* x1 .* x2)-rho0;
    g[2] = P .* (-x1+1) .* (1+x1 .^ 2) .^ 0.5 ./ (2 .* sqrt(2) .* x1 .* x2)-rho0;

    return f, g, h
end

# 11
function water_tesource_management(x)
    ## Water Resource Management Problem
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    ## objectives
    f = zeros(5);
    f[1] = 106780.37  .*  (x2 + x3) + 61704.67 ;
    f[2] = 3000  .*  x1 ;
    f[3] = 305700  .*  2289  .*  x2  ./  ((0.06 .* 2289) ^ 0.65) ;
    f[4] = 250  .*  2289  .*  exp(-39.75 .* x2+9.9 .* x3+2.74) ;
    f[5] = 25  .*  (1.39  ./ (x1 .* x2) + 4940 .* x3 -80) ;
    ## Constraints
    g = zeros(7);
    h = zeros(1);
    g[1] = 1 - (0.00139 ./ (x1 .* x2)+4.94 .* x3-0.08);
    g[2] = 1 - (0.000306 ./ (x1 .* x2)+1.082 .* x3-0.0986);
    g[3] = 50000 - (12.307 ./ (x1 .* x2) + 49408.24 .* x3+4051.02);
    g[4] = 16000 - (2.098 ./ (x1 .* x2)+8046.33 .* x3-696.71);
    g[5] = 10000 - (2.138 ./ (x1 .* x2)+7883.39 .* x3-705.04);
    g[6] = 2000 - (0.417 .* x1 .* x2 + 1721.26 .* x3-136.54);
    g[7] = 550 - (0.164 ./ (x1 .* x2)+631.13 .* x3-54.48);
    g = -g;

    return f, g, h
end

# 12
function simply_supported_I_beam_design(x)
    ## Simply supported I-beam design
    x1 = x[1];
    x2 = x[2];
    x3 = x[3];
    x4 = x[4];
    P = 600;
    L = 200;
    E = 2e4;
    ## objectives
    f = zeros(2);
    f[1] = 2  .*  x2  .*  x4 + x3  .*  (x1- 2 .* x4);
    f[2] = P .* L .^ 3 ./ (48 .* E .* (x3  .* ((x1 - 2 .* x4) .^ 3)+2 .* x2 .* x4 .* (4 .* x4 .* x4+3 .* x1 .* (x1-2 .* x4))) ./ 12);
    ## constraints
    g = zeros(1);
    h = zeros(1);
    g[1] = -16 + 180000*x1 ./ (x3 .* ((x1 - 2*x4) .^ 3)+2*x2 .* x4 .* (4 .* x4 .* x4+3*x1 .* (x1-2*x4))) + 15000 .* x2 ./ ((x1-2 .* x4) .* x3 .^ 3 + 2 .* x4 .* x2 .^ 3);

    return f, g, h
end

# 13
function gear_box_design(x)
    ## Gear box design
    x1 = x[1];
    x2 = x[2];
    x3 = round5(x[3]);
    x4 = x[4];
    x5 = x[5];
    x6 = x[6];
    x7 = x[7];
    ## objective function
    f = zeros(3);
    f[1] = 0.7854 .* x1 .* x2 .^ 2 .* (10 .* x3 .^ 2/3+14.933 .* x3-43.0934) - 1.508 .* x1 .* (x6 .^ 2+x7 .^ 2)+7.477 .* (x6 .^ 3+x7 .^ 3)+0.7854 .* (x4 .* x6 .^ 2+x5 .* x7 .^ 2);
    f[2] = sqrt((745 .* x4 ./ (x2 .* x3)) .^ 2+1.69e7) ./ (0.1 .* x6 .^ 3);
    f[3] = sqrt((745 .* x5 ./ (x2 .* x3)) .^ 2+1.575e8) ./ (0.1*x7 .^ 3);
    ## constraints
    g = zeros(11);
    h = zeros(1);
    g[1] = 1 ./ (x1 .* x2 .^ 2 .* x3)-1/27;
    g[2] = 1 ./ (x1 .* x2 .^ 2 .* x3 .^ 2)-1/397.5;
    g[3] = x4 .^ 3 ./ (x2 .* x3 .* x6 .^ 4)-1/1.93;
    g[4] = x5 .^ 3 ./ (x2 .* x3 .* x7 .^ 4)-1/1.93;
    g[5] = x2 .* x3-40;
    g[6] = x1 ./ x2-12;
    g[7] = -x1 ./ x2+5;
    g[8] = 1.9-x4+1.5 .* x6;
    g[9] = 1.9-x5+1.1 .* x7;
    g[10] = f[2]-1300;
    g[11] = f[3]-1100;

    return f, g, h
end

# 14
function multiple_disk_clutch_brake_design(x)
    ## Multiple-disk clutch brake design problem
    Mf = 3; Ms = 40; Iz = 55; n = 250; Tmax = 15; s = 1.5; delta = 0.5;
    Vsrmax = 10; rho = 0.0000078; pmax = 1; mu = 0.6; Lmax = 30; delR = 20;
    Rsr = 2 ./ 3 .* (x[2] .^ 3-x[1] .^ 3) ./ (x[2] .^ 2 .* x[1] .^ 2);
    Vsr = pi .* Rsr .* n ./ 30;
    A   = pi .* (x[2] .^ 2-x[1] .^ 2);
    Prz = x[4] ./ A;
    w   = pi .* n ./ 30;
    Mh  = 2/3 .* mu .* x[4] .* x[5] .* (x[2] .^ 3-x[1] .^ 3) ./ (x[2] .^ 2-x[1] .^ 2);
    T   = Iz .* w ./ (Mh+Mf);
    ## objective function
    f = zeros(2);
    f[1] = pi .* (x[2] .^ 2-x[1] .^ 2) .* x[3] .* (x[5]+1) .* rho;
    f[2] = T;
    ## constraints
    g = zeros(8)
    g[1] = -x[2]+x[1]+delR;
    g[2] = (x[5]+1) .* (x[3]+delta)-Lmax;
    g[3] = Prz-pmax;
    g[4] = Prz .* Vsr-pmax .* Vsrmax;
    g[5] = Vsr-Vsrmax;
    g[6] = T-Tmax;
    g[7] = s .* Ms-Mh;
    g[8] = -T;
    h      = zeros(1);

    return f, g, h
end

# 15
function spring_design(x)
    ## Spring design Problem
    x1 = round5(x[1]);
    x2 = x[2];
    d  = [0.009,0.0095,0.0104,0.0118,0.0128,0.0132,0.014,
          0.015, 0.0162, 0.0173, 0.018, 0.020, 0.023, 0.025,
          0.028, 0.032, 0.035, 0.041, 0.047, 0.054, 0.063,
          0.072, 0.080, 0.092, 0.0105, 0.120, 0.135, 0.148,
          0.162, 0.177, 0.192, 0.207, 0.225, 0.244, 0.263,
          0.283, 0.307, 0.331, 0.362,0.394,0.4375,0.500];
    x3 = d[max(1,min(42,round5(x[3])))];# x3 = x3(:);
    ## constants
    cf = (4 .* x2 ./ x3-1) ./ (4 .* x2 ./ x3-4)+0.615 .* x3 ./ x2;
    K  = (11.5 .* 10 .^ 6 .* x3 .^ 4) ./ (8 .* x1 .* x2 .^ 3);
    lf = 1000 ./ K + 1.05 .* (x1+2) .* x3;
    sigp = 300 ./ K;
    ## objective function
    f = zeros(2);
    f[1] = (pi .^ 2 .* x2 .* x3 .^ 2 .* (x1+2)) ./ 4;
    f[2] = (8000 .* cf .* x2) ./ (pi .* x3 .^ 3);
    ## constraints
    g = zeros(8)
    g[1] = (8000 .* cf .* x2) ./ (pi .* x3 .^ 3)-189000;
    g[2] = lf-14;
    g[3] = 0.2-x3;
    g[4] = x2-3;
    g[5] = 3-x2 ./ x3;
    g[6] = sigp - 6;
    g[7] = sigp+700 ./ K+1.05 .* (x1+2) .* x3-lf;
    g[8] = 1.25-700 ./ K;
    h = zeros(1);

    return f, g, h
end

# 16
function cantilever_beam_design(x)
    ## Cantilever beam design problem
    x1 = x[1];
    x2 = x[2];

    P = 1;
    E = 207000000;
    Sy = 300000;
    delta_max = 0.005;
    rho = 7800;

    ## objectives
    f = zeros(2);
    f[1] = 0.25  .*  rho  .*  pi  .*  x2  .*  x1 .^ 2;
    f[2] = (64  .*  P  .*  x2 .^ 3) ./ (3  .*  E  .*  pi  .*  x1 .^ 4);

    ## constraints
    g = zeros(2);
    h = zeros(1);
    g[1] = -Sy + (32  .*  P  .*  x2) ./ (pi  .*  x1 .^ 3);
    g[2] = -delta_max + (64  .*  P  .*  x2 .^ 3) ./ (3  .*  E  .*  pi  .*  x1 .^ 4);

    return f, g, h
end

# 17
function bulk_carriers_design(x)
    ## Bulk carriers design problem
    L = x[1];
    B = x[2];
    D = x[3];
    T = x[4];
    V_k = x[5];
    C_B = x[6];

    a = 4977.06 .* C_B .^ 2 - 8105.61 .* C_B + 4456.51;
    b = -10847.2 .* C_B .^ 2 + 12817 .* C_B - 6960.32;
    F_n = 0.5144 ./ (9.8065  .*  L) .^ 0.5;
    P = ((1.025 .* L .* B .* T .* C_B) .^ (2/3) .* V_k .^ 3) ./ (a + b .* F_n);

    W_s = 0.034 .* L .^ 1.7 .* B .^ 0.6 .* D .^ 0.4 .* C_B .^ 0.5;
    W_o = L .^ 0.8 .* B .^ 0.6 .* D .^ 0.3 .* C_B .^ 0.1;
    W_m = 0.17 .* P .^ 0.9;
    ls = W_s+W_o+W_m;

    D_wt = 1.025 .* L .* B .* T .* C_B-ls;
    F_c = 0.19 .* 24 .* P ./ 1000 + 0.2;
    D_cwt = D_wt - F_c .* ((5000 .* V_k) ./ 24 + 5)-2 .* D_wt .^ 0.5;
    R_trp = 350 ./ ((5000 .* V_k) ./ 24 + 2 .* (D_cwt ./ 8000 + 0.5));
    ac = D_cwt .* R_trp;
    S_d = 5000 .* V_k ./ 24;

    C_c = 0.2 .* 1.3  .*  (2000 .* W_s .^ 0.85 + 3500 .* W_o + 2400 .* P .^ 0.8);
    C_r = 40000 .* D_wt .^ 0.3;
    C_v = (1.05 .* 100 .* F_c .* S_d + 6.3 .* D_wt .^ 0.8) .* R_trp;

    ## objectives
    f = zeros(3);
    f[1]= (C_c + C_r + C_v) ./ ac;
    f[2] = ls;
    f[3] = -ac;
    ## constraints
    g = zeros(9);
    h = zeros(1);
    g[1] = L ./ B - 6;
    g[2]= 15 - L ./ D;
    g[3] = 19 - L ./ T;
    g[4] = 0.45 .* D_wt .^ 0.31 - T;
    g[5] = 0.7 .* D + 0.7 - T;
    g[6] = 0.32 - F_n;
    g[7] = 0.53 .* T + ((0.085 .* C_B - 0.002) .* B .^ 2) ./ (T .* C_B)-(1 + 0.52 .* D) - 0.07 .* B;
    g[8] = D_wt - 3000;
    g[9] = 500000 - D_wt;
    g = -g;

    return f, g, h
end

# 18
function front_rail_design(x)
    ## Front rail design problem
    hh = x[1];
    w = x[2];
    t = x[3];

    Ea = 14496.5;
    Fa = 234.9;
    E = -70973.4 + 958.656  .*  w+ 614.173  .*  hh - 3.827  .*  w  .*  hh + 57.023  .*  w  .*  t + 63.274  .*  hh  .*  t-3.582  .*  w .^ 2 - 1.4842  .*  hh .^ 2 - 1890.174  .*  t .^ 2;
    F = 111.854 - 20.210  .*  w + 7.560  .*  hh - 0.025  .*  w  .*  hh + 2.731  .*  w  .*  t - 1.479  .*  hh  .*  t + 0.165  .*  w .^ 2;

    ## objectives
    f = zeros(2);
    f[1] = Ea ./ E;
    f[2] = F ./ Fa;

    ## constraints
    g = zeros(3);
    h = zeros(1);
    g[1] = (hh - 136) .*  (146 - hh);
    g[2] = (w - 58)  .*  (66 - w);
    g[3] = (t - 1.4)  .*  (2.2 - t);
    g = -g;

    return f, g, h
end

# 19
function multi_product_batch_plant(x)
    ## Multi-product batch plant
    ## constant
    S = [2 3 4;
         4 6 3];
    t = [ 8 20 8;
         16 4 4];
    H = 6000; alp = 250; beta = 0.6;
    Q1 = 40000; Q2 = 20000;
    ## decision Variable
    N1 = round5(x[1]); N2 = round5(x[2]); N3 = round5(x[3]);
    V1 = x[4]; V2 = x[5]; V3 = x[6];
    TL1 = x[7]; TL2 = x[8];
    B1 = x[9]; B2 = x[10];
    ## objective function
    f = zeros(3);
    f[1] = alp .* (N1 .* V1 .^ beta+N2 .* V2 .^ beta+N3 .* V3 .^ beta);
    f[2] = 65 .* (Q1 ./ B1+Q2 ./ B2)+0.08 .* Q1+0.1 .* Q2;
    f[3] = Q1 .* TL1 ./ B1+Q2 .* TL2 ./ B2;
    ## constraints
    g = zeros(10)
    g[1] = Q1 .* TL1 ./ B1+Q2 .* TL2 ./ B2-H;
    g[2] = S[1,1] .* B1+S[2,1] .* B2-V1;
    g[3] = S[1,2] .* B1+S[2,2] .* B2-V2;
    g[4] = S[1,3] .* B1+S[2,3] .* B2-V3;
    g[5] = t[1,1]-N1 .* TL1;
    g[6] = t[1,2]-N2 .* TL1;
    g[7] = t[1,3]-N3 .* TL1;
    g[8] = t[2,1]-N1 .* TL2;
    g[9] = t[2,2]-N2 .* TL2;
    g[10] = t[2,3]-N3 .* TL2;
    h = zeros(1);

    return f, g, h
end

# 20
function hydro_static_thrust_bearing_design(x)
    ## Hydro-static thrust bearing design problem
    R = x[1]; Ro = x[2];  mu = x[3]; Q = x[4];
    gamma = 0.0307; C = 0.5; n = -3.55; C1 = 10.04;
    Ws = 101000; Pmax = 1000; delTmax = 50; hmin = 0.001;
    gg = 386.4; N = 750;
    P    = (log10(log10(8.122*1e6 .* mu+0.8))-C1) ./ n;
    delT = 2 .* (10 .^ P-560);
    Ef   = 9336 .* Q .* gamma .* C .* delT;
    h    = (2 .* pi .* N ./ 60) .^ 2 .* 2 .* pi .* mu ./ Ef .* (R .^ 4 ./ 4-Ro .^ 4 ./ 4)-1e-5;
    Po   = (6 .* mu .* Q ./ (pi .* h .^ 3)) .* log(R ./ Ro);
    W    = pi .* Po ./ 2 .* (R .^ 2-Ro .^ 2) ./ (log(R ./ Ro)-1e-5);
    ##  objective function
    f = zeros(2);
    f[1] = (Q .* Po ./ 0.7+Ef) ./ 12;
    f[2] = gamma ./ (gg .* Po) .* (Q ./ (2 .* pi .* R .* h));
    ##  constraints
    g = zeros(7)
    g[1] = Ws-W;
    g[2] = Po-Pmax;
    g[3] = delT-delTmax;
    g[4] = hmin-h;
    g[5] = Ro-R;
    g[6] = f[2]-0.001;
    g[7] = W ./ (pi .* (R .^ 2-Ro .^ 2)+1e-5)-5000;
    h      = zeros(1);

    return f, g, h
end

# 21
function crash_energy_management_for_high_speed_train(x)
    ## Crash energy management for high-speed train
    x1 = x[1]; x2 = x[2]; x3 = x[3];
    x4 = x[4]; x5 = x[5]; x6 = x[6];
    ## objective function
    f = zeros(2);
    f[1] = 1.3667145844797-0.00904459793976106*x1-0.0016193573938033*x2-0.00758531275221425*x3-0.00440727360327102*x4-0.00572216860791644*x5-0.00936039926190721*x6+2.62510221107328*10^(-6)*(x1 .^ 2)+4.92982681358861*10^(-7)*(x2 .^ 2)+2.25524989067108*10^(-6)*(x3 .^ 2)+
    1.84605439400301*10^(-6)*(x4 .^ 2)+2.17175358243416*10^(-6)*(x5 .^ 2)+3.90158043948054*10^(-6)*(x6 .^ 2)+4.55276994245781*10^(-7)*x1 .* x2-6.37013576290982*10^(-7)*x1 .* x3+8.26736480446359*10^(-7)*x1 .* x4+5.66352809442276*10^(-8)*x1 .* x5-3.20213897443278*10^(-7)*x1 .* x6+
    1.18015467772812*10^(-8)*x2 .* x3+9.25820391546515*10^(-8)*x2 .* x4-1.05705364119837*10^(-7)*x2 .* x5-4.74797783014687*10^(-7)*x2 .* x6-5.02319867013788*10^(-7)*x3 .* x4+9.54284258085225*10^(-7)*x3 .* x5+1.80533309229454*10^(-7)*x3 .* x6-1.07938022118477*10^(-6)*x4 .* x5-
    1.81370642220182*10^(-7)*x4 .* x6-2.24238851688047*10^(-7)*x5 .* x6;
    f[2] = -1.19896668942683+3.04107017009774*x1+1.23535701600191*x2+2.13882039381528*x3+2.33495178382303*x4+2.68632494801975*x5+3.43918953617606*x6-7.89144544980703*10^(-4)*(x1 .^ 2)-2.06085185698215*10^(-4)*(x2 .^ 2)-7.15269900037858*10^(-4)*(x3 .^ 2)-7.8449237573837*10^(-4)*(x4 .^ 2)-
    9.31396896237177*10^(-4)*(x5 .^ 2)-1.40826531972195*10^(-3)*(x6 .^ 2)-1.60434988248392*10^(-4)*x1 .* x2+2.0824655419411*10^(-4)*x1 .* x3-3.0530659653553*10^(-4)*x1 .* x4-8.10145973591615*10^(-5)*x1 .* x5+6.94728759651311*10^(-5)*x1 .* x6+1.18015467772812*10^(-8)*x2 .* x3+
    9.25820391546515*10^(-8)*x2 .* x4-1.05705364119837*10^(-7)*x2 .* x5+1.69935290196781*10^(-4)*x2 .* x6+2.32421829190088*10^(-5)*x3 .* x4-2.0808624041163476*10^(-4)*x3 .* x5+1.75576341867273*10^(-5)*x3 .* x6+2.68422081654044*10^(-4)*x4 .* x5+4.39852066801981*10^(-5)*x4 .* x6+
    2.96785446021357*10^(-5)*x5 .* x6;
    ## constraints
    g = zeros(4);
    h = zeros(1);
    g[1] = f[1]-5;
    g[2] = -f[1];
    g[3] = f[2] - 28;
    g[4] = -f[2];

    return f, g, h
end


text:
# This file is a part of JuliaPackageTemplate.jl, licensed under the MIT License (MIT).

using JuliaPackageTemplate
using Test

@testset "hello_world" begin
    @test JuliaPackageTemplate.hello_world() == 42
end


text:
"""
credits in making this file goes to the julia development team:
    David Varela      @00vareladavid
    David Widmann      @devmotion
    Eliott Saba        @staticfloat
    Fredrik Ekre       @fredrikekre
    Ian Butterworth    @IanButterworth
   Kristoffer Carlsson @KristofferC
                       @ianshmean
   Tim Holy            @timholy

"""

module SubdirTests
#=
import Pkg # ensure we are using the correct Pkg

using Pkg, UUIDs, Test
using Pkg.REPLMode: pkgstr
using Pkg.Types: PackageSpec
using Pkg: stdout_f, stderr_f
=#
using Utils

# turn /-paths into \-paths on Windows
function safe_path(path::String)
    if Sys.iswindows() && !isabspath(path)
        path = joinpath(split(path, "/")...)
    end
    return path
end

# turn /-paths into \-paths on Windows
function safe_path(path::String)
    if Sys.iswindows() && !isabspath(path)
        path = joinpath(split(path, "/")...)
    end
    return path
end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function closer(a::Aircraft, b::Aircraft, p::Aircraft)
	if (a.dyn.x - p.dyn.x)^2 + (a.dyn.y - p.dyn.y)^2 <
		(b.dyn.x - p.dyn.x)^2 + (b.dyn.y - p.dyn.y)^2
		return true
	else
		return false
	end
end

function farther(a::Aircraft, b::Aircraft, p::Aircraft)
	if (a.dyn.x - p.dyn.x)^2 + (a.dyn.y - p.dyn.y)^2 >
		(b.dyn.x - p.dyn.x)^2 + (b.dyn.y - p.dyn.y)^2
		return true
	else
		return false
	end
end

"""
	sort the Aircraft in an array
	using quick sort, returns ascending order
"""
function sort_agents!(p::Aircraft, agents_arr::Vector{Aircraft}, lo::Int64, hi::Int64)
    i, j = lo, hi
    while i < hi
        pivot = agents_arr[(lo + hi) >>> 1]
        while i <= j
            while closer(agents_arr[i], pivot, p); i += 1; end
            while farther(agents_arr[j], pivot, p); j -= 1; end
            if i <= j
                agents_arr[i], agents_arr[j] = agents_arr[j], agents_arr[i]
                i, j = i + 1, j - 1
            end
        end
        if lo < j; sort_agents!(p, agents_arr, lo, j); end
        lo, j = i, hi
    end
end

function get_section_id(ac::Aircraft, intruder::Aircraft, num_sections::Int64)
	ac_pos_angle = norm_angle(norm_angle(atan2(intruder.dyn.y - ac.dyn.y, intruder.dyn.x - ac.dyn.x)) - norm_angle(ac.dyn.heading))
	# ac_pos_angle = norm_angle(atan2(intruder.dyn.y - ac.dyn.y, intruder.dyn.x - ac.dyn.x))
	section = Int(floor(ac_pos_angle / (2 * pi / num_sections)) + 1)
	return section
end

function get_state(ac, intruder)
	intruder_pos_angle = atan2(intruder.dyn.y - ac.dyn.y, intruder.dyn.x - ac.dyn.x)

	return [ft2m(sqrt((ac.dyn.x - intruder.dyn.x)^2 + (ac.dyn.y - intruder.dyn.y)^2)),
			norm_angle(intruder_pos_angle - ac.dyn.heading),
			norm_angle(intruder.dyn.heading - ac.dyn.heading),
			ft2m(ac.dyn.v),
			ft2m(intruder.dyn.v)]
end

function get_full_and_low_fi_state(correction_type::Symbol, ac::Aircraft, All_AC::Vector{Aircraft}, num_sections::Int64)

	intruders = get_observation!(ac, All_AC, ac.sensor) # all intruders are within the sensing range

	if correction_type == :closest || policy
		# Method using num_sections closest intruders
		full_state = zeros(num_sections * STATE_DIM + AUG_STATE_DIM)
		hi = length(intruders)
		lo = (hi == 0) ? 0 : 1
		if hi > 0
			# # 1. Using multiple close intruders
			# # Ascending order in terms of distance from ego_agent
			sort_agents!(ac, intruders, lo, hi)
			# # 2. Using multiple dangerous intruders
			# # Ascending order in terms of maximum q_vals wrt ego_agent
			# # The minimum maximum q_vals => the most dangerous
			# sort_agents_q_vals!(env, intruders, lo, hi)

			for i in 1 : num_sections
				if i <= length(intruders)
					full_state[STATE_DIM * (i - 1) + 1 : STATE_DIM * i] .= get_state(ac, intruders[i])
				else
					full_state[STATE_DIM * (i - 1) + 1 : STATE_DIM * i] .= terminal_state
				end
			end
		else # no intruder
			full_state[1 : STATE_DIM * num_sections] .= TERM_VAR
		end

	elseif correction_type == :sector
		# # Method using num_sections sections
		section_ac_arr = Dict{Int64, Vector{Aircraft}}()
		for isection in 1:num_sections
			section_ac_arr[isection] = Aircraft[]
		end

		for int in intruders
			if int.id != ac.id
				section_id = get_section_id(ac, int, num_sections)
				push!(section_ac_arr[section_id], int)
			end
		end

		full_state = zeros(num_sections * STATE_DIM + AUG_STATE_DIM)

		# sort ac in each section
		for isection in 1:num_sections
			hi = length(section_ac_arr[isection])
			lo = (hi == 0) ? 0 : 1
			if hi > 0
				sort_agents!(ac, section_ac_arr[isection], lo, hi)
				# resulting section_ac_arr is in ascending order,
				# grab ac that is closest to the ego ac for each section
				# to update the full_state
				full_state[STATE_DIM * (isection - 1) + 1 : STATE_DIM * isection] .=
					get_state(ac, section_ac_arr[isection][1]) # [m]
			else
				# no intruder in the section
				# set the corresponding state var to terminal values
				full_state[STATE_DIM * (isection - 1) + 1 : STATE_DIM * isection] .=
					[TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR]
			end
		end
		#
	end # if correction_type

	# append the destination information to the full_state
	# 1. deviation from the destination, measured by rad
	deviation = norm_angle(ac.dyn.heading -
		norm_angle(atan2(ac.y_dest - ac.dyn.y, ac.x_dest - ac.dyn.x)))
	if deviation > π
		deviation = deviation - 2 * π # ∈ [-π, π]
		# deviation = 2 * π - deviation # ∈ [0, π], deprecated
	end

	# 2. current distance to destination
	dist_to_dest = sqrt((ac.dyn.y - ac.y_dest)^2 + (ac.dyn.x - ac.x_dest)^2) # [ft]
	# 3. previous distance to destination
	prev_dist_to_dest = ac.dist_to_dest # [ft]
	## update ego_agent_dist_to_dest
	ac.dist_to_dest = dist_to_dest # [ft]

	# update state entries
	full_state[end - AUG_STATE_DIM + 1 : end] .= [deviation, ft2m(dist_to_dest), ft2m(prev_dist_to_dest)] # [m]

	# find the lowfi_state from the full_state
	closest_intruder_section_id = indmin(full_state[[(i - 1) * STATE_DIM + IND_rho for i in 1:num_sections]])
	lowfi_state = full_state[STATE_DIM * (closest_intruder_section_id - 1) + 1 : STATE_DIM * closest_intruder_section_id]

	return full_state, lowfi_state
end

function normalize_obs!(ac::Aircraft, num_sections::Int64, obs::Vector{Float64})
	for isection in 0 : num_sections - 1
		if obs[isection * STATE_DIM + IND_rho] == TERM_VAR
			obs[isection * STATE_DIM + 1 : isection * STATE_DIM + STATE_DIM] *= 0 # if terminal, set empty
			# obs[isection * STATE_DIM + 1 : isection * STATE_DIM + STATE_DIM] = 1 # if terminal, set to 1
		else
			obs[isection * STATE_DIM + IND_rho] /= SENSING_RANGE
			obs[isection * STATE_DIM + IND_theta] /= (2 * pi)
			obs[isection * STATE_DIM + IND_phi] /= (2 * pi)
			# obs[isection * STATE_DIM + IND_v_own] = (obs[isection * STATE_DIM + IND_v_own] - V_MIN) / (V_MAX - V_MIN + 1e-12)
			obs[isection * STATE_DIM + IND_v_own] /= V_MAX
			obs[isection * STATE_DIM + IND_v_int] /= V_MAX
		end
	end
	obs[IND_deviation] /= pi
	obs[IND_dist_to_dest] /= ft2m(ac.nominal_route_len)
	obs[IND_prev_dist_to_dest] /= ft2m(ac.nominal_route_len)
end

"""
	update the observation of an Aircraft object with Nearest_Intruder_Tracker
"""
function get_observation!(ac::Aircraft, All_AC::Array{Aircraft, 1}, ::Nearest_Intruder_Tracker)
	# search for the closest intruder from ownship
	intruders = Aircraft[]
	closest_intruder = nothing
	closest_dist = Inf
	num_intruders_at_this_step = 0
	for other_ac in All_AC
		if other_ac.id != ac.id
			dist = sqrt((ac.dyn.x - other_ac.dyn.x)^2 + (ac.dyn.y - other_ac.dyn.y)^2)
			if dist <= ac.sensor.sensing_range
				push!(intruders, other_ac)
				num_intruders_at_this_step += 1
				ac.cumulative_num_ac_encountered += 1
				if dist < closest_dist
					closest_dist = dist
					closest_intruder = other_ac
				end
			end
		end
	end

	push!(ac.num_ac_encountered_arr, num_intruders_at_this_step)

	if closest_intruder == nothing
		ac.sensor.obs = []
	else
		ac.sensor.obs = [construct_obs(closest_dist, ac, closest_intruder)]
	end

	return intruders
end

"""
	Update the observation of an Aircraft object with Multi_Threat_Tracker
"""
function get_observation!(ac::Aircraft, All_AC::Array{Aircraft, 1}, ::Multi_Threat_Tracker)
	ac.sensor.obs = [] # dump previous obs
	intruders = Aircraft[]
	for other_ac in All_AC
		if other_ac.id != ac.id
			dist = sqrt((ac.dyn.x - other_ac.dyn.x)^2 + (ac.dyn.y - other_ac.dyn.y)^2)
		    if dist <= ac.sensor.sensing_range
				obs = construct_obs(dist, ac, other_ac)
		    	push!(ac.sensor.obs, obs)
		    	push!(intruders, other_ac)
		    end
		end
	end
	ac.cumulative_num_ac_encountered += length(intruders)
	push!(ac.num_ac_encountered_arr, length(intruders))
	return intruders
end

"""
	Helper: construct the Observation object
"""
function construct_obs(dist::Float64, ac::Aircraft, other_ac::Aircraft)
	rho = dist
	intruder_pos_angle = atan2(other_ac.dyn.y - ac.dyn.y, other_ac.dyn.x - ac.dyn.x)
	theta = intruder_pos_angle - ac.dyn.heading
	psi = other_ac.dyn.heading - ac.dyn.heading
	# constrain theta and psi
	if theta > pi
		theta = theta - 2 * pi
	elseif theta < - pi
		theta = theta + 2 * pi
	end
	if psi > pi
		psi = psi - 2 * pi
	elseif psi < - pi
		psi = psi + 2 * pi
	end

	sos = ac.dyn.v
	soi = other_ac.dyn.v
	tau = 0 # coaltitude
	pRA = ac.cas.Advisory_to_Ind_dict[ac.advisory]

	obs = Observation(rho, theta, psi, sos, soi, tau, pRA, intruder_pos_angle)
	return obs
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
# Tests for @snn_kw constructors

coretype(T::UnionAll) = coretype(T.body)
coretype(T::DataType) = T
paramnames(T) = coretype(T).parameters
function test_typeparams(Model; args=())
    Model = coretype(Model)
    n = Model(args...)
    for idx in 1:length(fieldnames(Model))
        fieldtypes(Model)[idx] isa TypeVar || continue
        field, Tf = fieldnames(Model)[idx], fieldtypes(Model)[idx].name

        # TODO: Test where vector is not <: DenseArray
        if Tf == :FT
            @test getfield(n, field) isa Float32
            _n = Model(args...;FT=Float64)
            @test getfield(_n, field) isa Float64
        elseif Tf == :VFT
            @test getfield(n, field) isa Vector{Float32}
            _n = Model(args...;VFT=Vector{Float64})
            @test getfield(_n, field) isa Vector{Float64}
        elseif Tf == :MFT
            @test getfield(n, field) isa Matrix{Float32}
            _n = Model(args...;MFT=Matrix{Float64})
            @test getfield(_n, field) isa Matrix{Float64}
        elseif Tf == :VBT
            @test getfield(n, field) isa Vector{Bool}
            _n = Model(args...;VBT=Vector{Any})
            @test getfield(_n, field) isa Vector{Any}
            @test getfield(_n, field)[1] isa Bool
        end
    end
end

@testset "Constructors" begin
@testset "Type parameters" begin
    for Model in (SNN.HH, SNN.IF, SNN.IF2, SNN.IZ, SNN.NoisyIF, SNN.Poisson, SNN.Rate)
        test_typeparams(Model)
    end
    test_typeparams(SNN.RateSynapse; args=(SNN.Rate(),SNN.Rate()))
    test_typeparams(SNN.PINningSynapse; args=(SNN.Rate(),SNN.Rate()))
    test_typeparams(SNN.FLSynapse; args=(SNN.Rate(),SNN.Rate()))
    test_typeparams(SNN.SpikingSynapse; args=(SNN.IF(),SNN.IF(),:ge))
end # Type parameters
end # Constructors