text:
abstract type GridType end
struct Inner <: GridType end
struct Outer <: GridType end

"""
3D grid with the following configuration. x and y are periodic coordinates,
uniformly spaced. The u coordinate is split into an "inner" and "outer" portion;
we therefore will have an inner and an outer grid. The outer grid is further
decomposed into several u-domains, and in each of these domains the u-direction
uses Gauss-Lobatto points.

Since the boundary conditions for the nested system are always specified at u=0,
the inner grid necessarily starts at u=0 and finishes at u=u_outer_min.
"""
Base.@kwdef struct SpecCartGrid3D{T<:Real}
    x_min            :: T
    x_max            :: T
    x_nodes          :: Int
    y_min            :: T
    y_max            :: T
    y_nodes          :: Int
    u_outer_min      :: T
    u_outer_max      :: T
    u_outer_domains  :: Int     = 1
    u_outer_nodes    :: Int # number of points per domain
    u_inner_nodes    :: Int
    fd_order         :: Int     = 4
    filter_γ         :: T       = 8.0
    sigma_diss       :: T       = 0.2
end

function Jecco.Atlas(grid::SpecCartGrid3D{T}) where {T}
    u_inner_coord = GaussLobatto{1}("u", zero(T), grid.u_outer_min, grid.u_inner_nodes)

    N_outer_sys = grid.u_outer_domains
    delta_udom  = (grid.u_outer_max - grid.u_outer_min) / N_outer_sys

    u_outer_coords =
        [GaussLobatto{1}("u", grid.u_outer_min + (i-1)*delta_udom,
                          grid.u_outer_min + i*delta_udom, grid.u_outer_nodes)
         for i in 1:N_outer_sys]

    xcoord  = Cartesian{2}("x", grid.x_min, grid.x_max, grid.x_nodes, endpoint=false)
    ycoord  = Cartesian{3}("y", grid.y_min, grid.y_max, grid.y_nodes, endpoint=false)

    inner_chart  = Chart(u_inner_coord, xcoord, ycoord)
    outer_charts = [Chart(u_outer_coords[i], xcoord, ycoord)
                     for i in 1:N_outer_sys]

    Atlas([inner_chart; outer_charts])
end

struct Filters{F1}
    exp_filter    :: F1
end
function Filters(filter_gamma::T, KO_order::Int, sigma_diss::T,
                 Nu::Int, Nx::Int, Ny::Int) where {T}
    exp_filter    = Exp_Filter{1}(filter_gamma, Nu, Nx, Ny)
    Filters{typeof(exp_filter)}(exp_filter)
end

struct System{GT,Cu,Cx,Cy,TDu,TDx,TDy,TI,TF,TDKOx,TDKOy}
    gridtype    :: GT
    ucoord      :: Cu
    xcoord      :: Cx
    ycoord      :: Cy
    Du          :: TDu
    Duu         :: TDu
    Dx          :: TDx
    Dxx         :: TDx
    Dy          :: TDy
    Dyy         :: TDy
    uinterp     :: TI
    filters     :: TF
    DKOx        :: TDKOx
    DKOy        :: TDKOy
end

function System(gridtype::GT, ucoord::GaussLobattoCoord,
                xcoord::CartesianCoord, ycoord::CartesianCoord, ord::Int,
                filter_gamma::T, sigma_diss::T) where {GT<:GridType,T<:Real}
    Du  = ChebDeriv{1}(1, ucoord.min, ucoord.max, ucoord.nodes)
    Duu = ChebDeriv{1}(2, ucoord.min, ucoord.max, ucoord.nodes)

    Dx  = CenteredDiff{2}(1, ord, Jecco.delta(xcoord), xcoord.nodes)
    Dxx = CenteredDiff{2}(2, ord, Jecco.delta(xcoord), xcoord.nodes)

    Dy  = CenteredDiff{3}(1, ord, Jecco.delta(ycoord), ycoord.nodes)
    Dyy = CenteredDiff{3}(2, ord, Jecco.delta(ycoord), ycoord.nodes)

    uinterp = ChebInterpolator(ucoord.min, ucoord.max, ucoord.nodes)

    KO_order = ord + 1
    filters  = Filters(filter_gamma, KO_order, sigma_diss, ucoord.nodes,
                       xcoord.nodes, ycoord.nodes)

    DKOx  = KO_Centered{2}(KO_order, sigma_diss, Jecco.delta(xcoord), xcoord.nodes)
    DKOy  = KO_Centered{3}(KO_order, sigma_diss, Jecco.delta(ycoord), ycoord.nodes)

    System{GT, typeof(ucoord), typeof(xcoord), typeof(ycoord), typeof(Du),
           typeof(Dx), typeof(Dy), typeof(uinterp),
           typeof(filters),
           typeof(DKOx), typeof(DKOy)}(gridtype, ucoord, xcoord, ycoord,
                                       Du, Duu, Dx, Dxx, Dy, Dyy, uinterp, filters,
                                       DKOx, DKOy)
end

Base.size(sys::System) = (sys.ucoord.nodes, sys.xcoord.nodes, sys.ycoord.nodes)

struct SystemPartition{N,A} <: AbstractPartition{N,A}
    x :: NTuple{N,A}
end

function SystemPartition(x_::AbstractVector{S}) where {S<:System}
    x = Tuple(x_)
    N = length(x)
    SystemPartition{N,eltype(x)}(x)
end

"""
    SystemPartition(grid::SpecCartGrid3D)

Create a `SystemPartition`, which is a `Tuple` of `System` where each element
corresponds to a different u-domain
"""
function SystemPartition(grid::SpecCartGrid3D{T}) where {T}
    u_inner_coord = GaussLobatto{1}("u", zero(T), grid.u_outer_min, grid.u_inner_nodes)

    N_outer_sys = grid.u_outer_domains
    delta_udom  = (grid.u_outer_max - grid.u_outer_min) / N_outer_sys

    u_outer_coords =
        [GaussLobatto{1}("u", grid.u_outer_min + (i-1)*delta_udom,
                          grid.u_outer_min + i*delta_udom, grid.u_outer_nodes)
         for i in 1:N_outer_sys]

    xcoord  = Cartesian{2}("x", grid.x_min, grid.x_max, grid.x_nodes, endpoint=false)
    ycoord  = Cartesian{3}("y", grid.y_min, grid.y_max, grid.y_nodes, endpoint=false)

    inner_system = System(Inner(), u_inner_coord, xcoord, ycoord, grid.fd_order,
                          grid.filter_γ, grid.sigma_diss)

    outer_systems = [System(Outer(), u_outer_coords[i], xcoord, ycoord, grid.fd_order,
                            grid.filter_γ, grid.sigma_diss) for i in 1:N_outer_sys]

    SystemPartition([inner_system; outer_systems])
end

"""
    Boundary(grid::SpecCartGrid3D)

Create a `Boundary` struct with arrays of `size = (1,p.x_nodes,p.y_nodes)`
"""
function Boundary(grid::SpecCartGrid3D{T}) where {T}
    Nx = grid.x_nodes
    Ny = grid.y_nodes
    Boundary{T}(undef, Nx, Ny)
end

"""
    Gauge(grid::SpecCartGrid3D)

Create a `Gauge` struct with arrays of `size = (1,grid.x_nodes,grid.y_nodes)`
"""
function Gauge(grid::SpecCartGrid3D{T}) where {T}
    Nx = grid.x_nodes
    Ny = grid.y_nodes
    Gauge{T}(undef, Nx, Ny)
end

"""
    BulkEvolvedPartition(grid::SpecCartGrid3D)

Returns `BulkPartition`, of `length = 1 + grid.u_outer_domains`, with
elements of type `BulkEvolved`. The first `BulkEvolved` has arrays of `size =
(grid.u_inner_nodes, grid.x_nodes, grid.y_nodes)`, and the remaining ones have
`size = (grid.u_outer_nodes, grid.x_nodes, grid.y_nodes)`
"""
function BulkEvolvedPartition(grid::SpecCartGrid3D{T}) where {T}
    Nx = grid.x_nodes
    Ny = grid.y_nodes
    Nu_in  = grid.u_inner_nodes
    Nu_out = grid.u_outer_nodes
    N_outer_sys = grid.u_outer_domains

    bulk_in   = BulkEvolved{T}(undef, Nu_in, Nx, Ny)
    bulk_out  = [BulkEvolved{T}(undef, Nu_out, Nx, Ny) for i in 1:N_outer_sys]

    BulkPartition((bulk_in, bulk_out...))
end

"""
    BulkConstrainedPartition(grid::SpecCartGrid3D)

Returns `BulkPartition`, of `length = 1 + grid.u_outer_domains`, with
elements of type `BulkConstrained`. The first `BulkConstrained` has arrays of `size =
(grid.u_inner_nodes, grid.x_nodes, grid.y_nodes)`, and the remaining ones have
`size = (grid.u_outer_nodes, grid.x_nodes, grid.y_nodes)`
"""
function BulkConstrainedPartition(grid::SpecCartGrid3D{T}) where {T}
    Nx = grid.x_nodes
    Ny = grid.y_nodes
    Nu_in  = grid.u_inner_nodes
    Nu_out = grid.u_outer_nodes
    N_outer_sys = grid.u_outer_domains

    bulk_in   = BulkConstrained{T}(undef, Nu_in, Nx, Ny)
    bulk_out  = [BulkConstrained{T}(undef, Nu_out, Nx, Ny) for i in 1:N_outer_sys]

    BulkPartition((bulk_in, bulk_out...))
end
"""
    BulkConstrainedPartition(systems::SystemPartition)

Returns `BulkPartition`, of `length = length(systems)`
"""
function BulkConstrainedPartition(systems::SystemPartition)
    T     = Jecco.coord_eltype(systems[1].ucoord)
    bulks = [BulkConstrained{T}(undef, size(sys)...) for sys in systems]

    BulkPartition(bulks...)
end

function BulkDerivPartition(grid::SpecCartGrid3D{T}) where {T}
    Nx = grid.x_nodes
    Ny = grid.y_nodes
    Nu_in  = grid.u_inner_nodes
    Nu_out = grid.u_outer_nodes
    N_outer_sys = grid.u_outer_domains

    bulk_in   = BulkDeriv{T}(undef, Nu_in, Nx, Ny)
    bulk_out  = [BulkDeriv{T}(undef, Nu_out, Nx, Ny) for i in 1:N_outer_sys]

    BulkPartition((bulk_in, bulk_out...))
end
function BulkDerivPartition(systems::SystemPartition)
    T     = Jecco.coord_eltype(systems[1].ucoord)
    bulks = [BulkDeriv{T}(undef, size(sys)...) for sys in systems]

    BulkPartition(bulks...)
end

function HorizonCache(sys::System, ord::Int)
    _, Nx, Ny = size(sys)
    hx = Jecco.delta(sys.xcoord)
    hy = Jecco.delta(sys.ycoord)
    T1 = Jecco.coord_eltype(sys.ucoord)
    M  = Nx * Ny

    # there is no method for lu when acting on SparseMatrixCSC{BigFloat,Int64}, so
    # let's use instead type Float64 for such cases
    if T1 == BigFloat
        T2 = Float64
    else
        T2 = T1
    end

    bulkhorizon = BulkHorizon{T1}(Nx, Ny)

    axx    = Vector{T2}(undef, M)
    ayy    = Vector{T2}(undef, M)
    axy    = Vector{T2}(undef, M)
    bx     = Vector{T2}(undef, M)
    by     = Vector{T2}(undef, M)
    cc     = Vector{T2}(undef, M)
    b_vec  = Vector{T2}(undef, M)

    Dx_    = CenteredDiff{1}(1, ord, hx, Nx)
    Dxx_   = CenteredDiff{1}(2, ord, hx, Nx)
    Dy_    = CenteredDiff{2}(1, ord, hy, Ny)
    Dyy_   = CenteredDiff{2}(2, ord, hy, Ny)

    #=
    use the Kronecker product (kron) to build 2-dimensional derivation matrices
    from 1-dimensional ones. see for instance:

    https://en.wikipedia.org/wiki/Kronecker_product
    https://arxiv.org/pdf/1801.01483.pdf (section 5)
    =#
    Dx_2D  = T2.(kron(I(Ny), SparseMatrixCSC(Dx_)))
    Dxx_2D = T2.(kron(I(Ny), SparseMatrixCSC(Dxx_)))
    Dy_2D  = T2.(kron(SparseMatrixCSC(Dy_), I(Nx)))
    Dyy_2D = T2.(kron(SparseMatrixCSC(Dyy_), I(Nx)))
    Dxy_2D = T2.(Dx_2D * Dy_2D)

    _Dx_2D  = copy(Dx_2D)
    _Dxx_2D = copy(Dxx_2D)
    _Dy_2D  = copy(Dy_2D)
    _Dyy_2D = copy(Dyy_2D)
    _Dxy_2D = copy(Dxy_2D)

    HorizonCache{T1,T2,typeof(Dx_2D)}(bulkhorizon, axx, ayy, axy, bx, by, cc, b_vec,
                                  Dx_2D,  Dy_2D, Dxx_2D,  Dyy_2D,  Dxy_2D,
                                  _Dx_2D, _Dy_2D, _Dxx_2D, _Dyy_2D, _Dxy_2D)
end


text:
"""
T5 fileformat
See EUROCONTROL NEST Manual for T5 fileformat description

T5 DataFrame Column Names

Field               Type    Size    Comment
:FLIGHT_ID           Int64           Flight ID
:SECTOR_NAME         String
:ENTRY_DATETIME      DateTime
:EXIT_DATETIME       DateTime
:ENTRY_FL            Float64         Flight level in decimal
:EXIT_FL             Float64         Flight level in decimal
:ENTRY_SEGMENT_NAME  String          Separator between point name is '_'
                                    '##' as a prefix means that the flight was
                                    already present in the sector before it was
                                    opened
:EXIT_SEGMENT_NAME   String          Separator between point name is '_', might
                                    be different from ENTRY_SEGMENT_NAME
                                    '##' as a prefix means that the flight was
                                    still present in the sector after it was
                                    closed.
:TOT_DISTANCE_IN_SEGMENT_M Float64   For that FLIGHT_ID in that SECTOR_NAME,
                                    expressed in meters.
:TOT_TIME_IN_SEGMENT_S Float64       For that FLIGHT_ID in that SECTOR_NAME,
                                    expressed in seconds.

Example T5 file format
213765625 EHAMCR1 1513288140.000000 1513288226.000000 0.000000 35.000000 0.000000 1.669195 EHAM_*AM80 !AAEW_!AAEX 4.749612 86.000000
213765625 EHAAFIR 1513288140.000000 1513288980.000000 0.000000 255.000000 0.000000 14.478077 EHAM_*AM80 !AAEb_SONEB 84.060844 840.000000
213765625 EHALL 1513288140.000000 1513288980.000000 0.000000 255.000000 0.000000 14.478077 EHAM_*AM80 !AAEb_SONEB 84.060844 840.000000
213765625 EH_DN 1513288140.000000 1513288982.000000 0.000000 255.361450 0.000000 0.247570 EHAM_*AM80 SONEB_OLDOD 84.308411 842.000000
"""
module T5

export read

using Format
using CSV
using Dates
using DataFrames

const fileformat = Dict(1=>Int64, 2=>String, 3=>Float64, 4=>Float64, 5=>Float64,
6=>Float64, 7=>Float64, 8=>Float64, 9=>String, 10=>String, 11=>Float64,
12=>Float64)

const header_format = ["FLIGHT_ID", "SECTOR_NAME", "ENTRY_DATETIME",
"EXIT_DATETIME", "ENTRY_FL", "EXIT_FL", "ENTRY_DISTANCE_M", "EXIT_DISTANCE_M",
"ENTRY_SEGMENT_NAME", "EXIT_SEGMENT_NAME", "TOT_DISTANCE_IN_SEGMENT_M",
"TOT_TIME_IN_SEGMENT_S"]

function read(file)
    df = CSV.read(file, delim=" ", types=fileformat, header=header_format,
    copycols=true)
    reformat!(df);
    return df
end

function reformat!(df)
    #Date conversion
    df[:,:TEMP] = df[:,:ENTRY_DATETIME]
    select!(df, Not(:ENTRY_DATETIME))
    df[:,:ENTRY_DATETIME] =  unix2datetime.(df[:,:TEMP])
    select!(df, Not(:TEMP))

    df[:,:TEMP] = df[:,:EXIT_DATETIME]
    select!(df, Not(:EXIT_DATETIME))
    df[:,:EXIT_DATETIME] =  unix2datetime.(df[:,:TEMP])
    select!(df, Not(:TEMP))

    df[:,:ENTRY_DISTANCE_M] = df[:,:ENTRY_DISTANCE_M] * 1852.0
    df[:,:EXIT_DISTANCE_M] = df[:,:EXIT_DISTANCE_M] * 1852.0
    df[:,:TOT_DISTANCE_IN_SEGMENT_M] = df[:,:TOT_DISTANCE_IN_SEGMENT_M] * 1852.0

    select!(df,[:FLIGHT_ID, :SECTOR_NAME, :ENTRY_DATETIME, :EXIT_DATETIME,
    :ENTRY_FL, :EXIT_FL, :ENTRY_SEGMENT_NAME, :EXIT_SEGMENT_NAME,
    :TOT_DISTANCE_IN_SEGMENT_M, :TOT_TIME_IN_SEGMENT_S])
end

end # module

ext: jl


text:
module ReinforcementLearningExperiments

using Reexport

@reexport using ReinforcementLearning

const EXPERIMENTS_DIR = joinpath(@__DIR__, "experiments")
# for f in readdir(EXPERIMENTS_DIR)
#     include(joinpath(EXPERIMENTS_DIR, f))
# end
include(joinpath(EXPERIMENTS_DIR, "JuliaRL_BasicDQN_CartPole.jl"))
include(joinpath(EXPERIMENTS_DIR, "JuliaRL_DQN_CartPole.jl"))
include(joinpath(EXPERIMENTS_DIR, "JuliaRL_PrioritizedDQN_CartPole.jl"))
include(joinpath(EXPERIMENTS_DIR, "JuliaRL_QRDQN_CartPole.jl"))

# dynamic loading environments
function __init__() end

end # module


text:
# syntax: proto3
mutable struct Method <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Method(; kwargs...)
        obj = new(meta(Method), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Method
const __meta_Method = Ref{ProtoMeta}()
function meta(::Type{Method})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Method)
            __meta_Method[] = target = ProtoMeta(Method)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :request_type_url => AbstractString, :request_streaming => Bool, :response_type_url => AbstractString, :response_streaming => Bool, :options => Base.Vector{Option}, :syntax => Int32]
            meta(target, Method, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Method[]
    end
end
function Base.getproperty(obj::Method, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :request_type_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :request_streaming
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :response_type_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :response_streaming
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    elseif name === :syntax
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct Mixin <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Mixin(; kwargs...)
        obj = new(meta(Mixin), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Mixin
const __meta_Mixin = Ref{ProtoMeta}()
function meta(::Type{Mixin})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Mixin)
            __meta_Mixin[] = target = ProtoMeta(Mixin)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :root => AbstractString]
            meta(target, Mixin, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Mixin[]
    end
end
function Base.getproperty(obj::Mixin, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :root
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct Api <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Api(; kwargs...)
        obj = new(meta(Api), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Api
const __meta_Api = Ref{ProtoMeta}()
function meta(::Type{Api})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Api)
            __meta_Api[] = target = ProtoMeta(Api)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :methods => Base.Vector{Method}, :options => Base.Vector{Option}, :version => AbstractString, :source_context => SourceContext, :mixins => Base.Vector{Mixin}, :syntax => Int32]
            meta(target, Api, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Api[]
    end
end
function Base.getproperty(obj::Api, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :methods
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Method}
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    elseif name === :version
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :source_context
        return (obj.__protobuf_jl_internal_values[name])::SourceContext
    elseif name === :mixins
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Mixin}
    elseif name === :syntax
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export Api, Method, Mixin


text:
# Autogenerated wrapper script for LibArchive_jll for x86_64-linux-musl
export bsdcat, bsdcpio, bsdtar, libarchive

using acl_jll
using Attr_jll
using Bzip2_jll
using Expat_jll
using Libiconv_jll
using Lz4_jll
using OpenSSL_jll
using XZ_jll
using Zlib_jll
using Zstd_jll
JLLWrappers.@generate_wrapper_header("LibArchive")
JLLWrappers.@declare_executable_product(bsdcat)
JLLWrappers.@declare_executable_product(bsdcpio)
JLLWrappers.@declare_executable_product(bsdtar)
JLLWrappers.@declare_library_product(libarchive, "libarchive.so.13")
function __init__()
    JLLWrappers.@generate_init_header(acl_jll, Attr_jll, Bzip2_jll, Expat_jll, Libiconv_jll, Lz4_jll, OpenSSL_jll, XZ_jll, Zlib_jll, Zstd_jll)
    JLLWrappers.@init_executable_product(
        bsdcat,
        "bin/bsdcat",
    )

    JLLWrappers.@init_executable_product(
        bsdcpio,
        "bin/bsdcpio",
    )

    JLLWrappers.@init_executable_product(
        bsdtar,
        "bin/bsdtar",
    )

    JLLWrappers.@init_library_product(
        libarchive,
        "lib/libarchive.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JuMP, Base.Test, AmplNLWriter

# Example with no objective (#50)

# solver = AmplNLSolver(Ipopt.amplexe, ["print_level=0"])

@testset "example: jump_no_obj" begin
    m = Model(solver=solver)
    @variable(m, 0 <= yp <= 1, Int)
    @variable(m, 0 <= l <= 1000.0)
    @variable(m, 0 <= f <= 1000.0)
    @NLconstraint(m, .087 * l >= f ^ 2)
    @constraint(m, l <= yp * 1000.0)

    @test solve(m) == :Optimal
    @test getobjectivevalue(m) == 0
end


text:
# ---
# title: 92. Reverse Linked List II
# id: problem92
# author: Qling
# date: 2021-03-18
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-linked-list-ii/description/>
# hidden: true
# ---
#
# Reverse a linked list from position _m_ to _n_. Do it in one-pass.
#
# **Note:  **1 ≤ _m_ ≤ _n_ ≤ length of list.
#
# **Example:**
#
#
#
#     Input: 1->2->3->4->5->NULL, _m_ = 2, _n_ = 4
#     Output: 1->4->3->2->5->NULL
#
#
#
## @lc code=start
using LeetCode

function reverse_between(head::ListNode{Int}, m::Int, n::Int)::ListNode
    (isnothing(head) || m == n) && return head
    p = dummy = ListNode()
    next!(dummy, head)
    for _ = 1:m-1
        p = next(p)
    end

    tail = next(p)

    for _ = 1:n-m
        tmp = next(p)
        next!(p, next(tail))
        next!(tail, next(next(tail)))
        next!(next(p), tmp)
    end

    return next(dummy)
end
## @lc code=end


text:
module BA_5B

using NBodyIPs.FastPolys
using StaticArrays

import NBodyIPs.BAInvariants: invariants, invariants_d, invariants_ed, tdegrees, rtdegrees

# TODO: THIS IS A HACK TO GET the `ninvariants` function right, but it cannot
#       yet be used to actually produce basis functions
# tdegrees(::Val{5}) = ntuple(_->10, 10), ntuple(_->144, 10)

tdegrees(::Val{5}) =
   (1, 1, 2, 2, 2, 3, 3, 3, 4, 4), (0, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12)

# length-angle degrees (small hack in the last primary invariant)
rtdegrees(::Val{5}) =
   ( (1, 0), (0, 1), (2, 0), (0, 2), (1, 1), (3, 0), (0, 3), (0, 3), (0, 4), (4, 4) ),
   ( (0, 0), (0, 2), (2, 1), (2, 1), (1, 2), (1, 2), (0, 3), (0, 4), (4, 0), (3, 1), (2, 2), (2, 2), (2, 2), (1, 3), (1, 3), (1, 3), (1, 4), (1, 4), (2, 3), (0, 5), (2, 3), (4, 1), (3, 2), (2, 3), (2, 3), (1, 4), (1, 4), (0, 5), (2, 4), (2, 4), (1, 5), (4, 2), (1, 5), (2, 4), (2, 4), (0, 6), (3, 3), (2, 4), (0, 6), (4, 2), (3, 3), (1, 5), (1, 5), (3, 3), (4, 2), (1, 5), (5, 2), (1, 6), (2, 5), (2, 5), (3, 4), (4, 3), (5, 2), (4, 3), (2, 5), (2, 5), (3, 4), (4, 3), (2, 5), (3, 4), (3, 4), (2, 5), (3, 4), (0, 7), (1, 6), (1, 6), (1, 6), (2, 6), (3, 5), (1, 7), (4, 4), (2, 6), (3, 5), (2, 6), (4, 4), (6, 2), (3, 5), (3, 5), (5, 3), (3, 5), (1, 7), (2, 6), (6, 2), (5, 3), (3, 5), (4, 4), (1, 7), (4, 4), (4, 4), (1, 7), (0, 8), (5, 4), (3, 6), (5, 4), (4, 5), (1, 8), (2, 7), (5, 4), (2, 7), (6, 3), (2, 7), (5, 4), (6, 3), (4, 5), (3, 6), (4, 5), (1, 8), (2, 7), (2, 7), (7, 2), (2, 7), (3, 6), (3, 6), (5, 4), (0, 9), (1, 8), (3, 6), (3, 7), (4, 6), (5, 5), (3, 7), (3, 7), (3, 7), (1, 9), (5, 5), (2, 8), (2, 8), (6, 4), (4, 6), (5, 5), (3, 7), (6, 4), (6, 5), (1, 10), (5, 6), (2, 9), (3, 8), (5, 6), (2, 9), (0, 11), (3, 9), (5, 7), (6, 6), (3, 9) )

const G_BA_5B = [
[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ]
,[ 3, 1, 4, 2, 6, 10, 8, 7, 5, 9 ]
,[ 2, 4, 1, 3, 9, 5, 8, 7, 10, 6 ]
,[ 4, 1, 2, 3, 7, 9, 10, 5, 6, 8 ]
,[ 1, 3, 4, 2, 6, 7, 5, 10, 8, 9 ]
,[ 3, 4, 2, 1, 10, 8, 6, 9, 7, 5 ]
,[ 2, 1, 3, 4, 5, 8, 9, 6, 7, 10 ]
,[ 4, 2, 3, 1, 9, 10, 7, 8, 5, 6 ]
,[ 1, 4, 2, 3, 7, 5, 6, 9, 10, 8 ]
,[ 3, 2, 1, 4, 8, 6, 10, 5, 9, 7 ]
,[ 2, 3, 4, 1, 8, 9, 5, 10, 6, 7 ]
,[ 4, 3, 1, 2, 10, 7, 9, 6, 8, 5 ]
,[ 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 ]
,[ 3, 4, 1, 2, 10, 6, 8, 7, 9, 5 ]
,[ 2, 1, 4, 3, 5, 9, 8, 7, 6, 10 ]
,[ 4, 2, 1, 3, 9, 7, 10, 5, 8, 6 ]
,[ 1, 4, 3, 2, 7, 6, 5, 10, 9, 8 ]
,[ 3, 2, 4, 1, 8, 10, 6, 9, 5, 7 ]
,[ 2, 3, 1, 4, 8, 5, 9, 6, 10, 7 ]
,[ 4, 3, 2, 1, 10, 9, 7, 8, 6, 5 ]
,[ 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 ]
,[ 3, 1, 2, 4, 6, 8, 10, 5, 7, 9 ]
,[ 2, 4, 3, 1, 9, 8, 5, 10, 7, 6 ]
,[ 4, 1, 3, 2, 7, 10, 9, 6, 5, 8 ]
,]
simplex_permutations(x::SVector{10}) = [x[G_BA_5B[i]] for i=1:24]
# Primary invariants for BA_5B
 # : definitions at the beginning of the file
const P1_1 = (1,2,3,4,)

const P2_1 = (5,9,6,8,7,10,)

const P3_1 = (1,2,3,4,)

const P4_1 = (5,9,6,8,7,10,)

const P5_1 = (1,2,3,2,1,4,2,4,1,3,4,3,)
const P5_2 = (5,9,6,8,7,10,5,9,6,8,7,10,)

const P6_1 = (1,2,3,4,)

const P7_1 = (5,5,6,7,)
const P7_2 = (6,8,8,9,)
const P7_3 = (7,9,10,10,)

const P8_1 = (5,9,6,8,7,10,)

const P9_1 = (5,9,6,8,7,10,)

const P10_1 = (1,2,3,4,5,5,5,5,6,7,8,9,6,6,6 ,7,8,10,7,7 ,9,10,8,8 ,9,10,9 ,10)
const P10_2 = (1,2,3,4,6,7,8,9,5,5,5,5,7,8,10,6,6,6 ,9,10,7,7 ,9,10,8,8 ,10,9)

# Irreducible secondaries for group BA_5B
 # : definitions at the beginning of the file
const IS1_1 = (5,5,6,8,5,7,5,9,6,6,7,8,)
const IS1_2 = (6,9,10,9,7,10,8,10,7,8,9,10,)

const IS2_1 = (1,2,1,2,1,3,)
const IS2_2 = (2,4,3,3,4,4,)
const IS2_3 = (10,6,9,7,8,5,)

const IS3_1 = (1,2,1,2,1,3,)
const IS3_2 = (2,4,3,3,4,4,)
const IS3_3 = (5,9,6,8,7,10,)

const IS4_1 = (8,7,7,10,9,6,6,8,10,5,5,9,)
const IS4_2 = (1,2,3,2,1,4,2,4,1,3,4,3,)

const IS5_1 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS5_2 = (5,7,6,8,7,6,5,8,6,5,5,9,6,7,5,8,7,5,5,8,5,6,6,9,)
const IS5_3 = (8,9,7,10,9,10,6,9,10,8,7,10,8,10,7,9,10,9,6,10,9,8,7,10,)

const IS6_1 = (5,9,6,8,7,10,5,9,6,8,7,10,6,10,5,8,7,9,6,10,5,8,7,9,)
const IS6_2 = (6,5,10,9,5,7,8,10,7,6,9,8,5,6,9,10,6,7,8,9,7,5,10,8,)

const IS7_1 = (1,2,3,2,1,4,2,4,1,3,4,3,)
const IS7_2 = (2,1,1,3,2,1,1,2,3,1,1,2,)
const IS7_3 = (3,4,4,4,4,3,3,3,4,2,2,4,)

const IS8_1 = (1,2,3,2,1,4,2,4,1,3,4,3,)
const IS8_2 = (2,4,1,3,4,3,1,2,3,2,1,4,)
const IS8_3 = (10,6,9,7,8,5,10,6,9,7,8,5,)

const IS9_1 = (1,2,1,2,1,3,)
const IS9_2 = (2,4,3,3,4,4,)
const IS9_3 = (5,6,6,7,7,5,)
const IS9_4 = (10,9,9,8,8,10,)

const IS10_1 = (1,2,1,2,1,3,1,2,1,2,1,3,1,3,1,2,1,2,1,3,1,2,1,2,)
const IS10_2 = (2,4,3,3,4,4,2,4,3,3,4,4,3,4,2,3,4,4,3,4,2,3,4,4,)
const IS10_3 = (6,5,9,7,5,5,8,6,7,6,8,5,5,5,9,7,6,6,8,5,7,5,8,6,)
const IS10_4 = (10,6,10,9,8,7,10,10,9,7,9,8,9,6,10,10,8,7,9,9,10,7,10,8,)

const IS11_1 = (5,9,6,8,7,10,)
const IS11_2 = (1,2,1,2,1,3,)
const IS11_3 = (2,4,3,3,4,4,)

const IS12_1 = (10,6,9,7,8,5,10,6,9,7,8,5,)
const IS12_2 = (1,2,3,2,1,4,2,4,1,3,4,3,)
const IS12_3 = (5,9,6,8,7,10,5,9,6,8,7,10,)

const IS13_1 = (8,7,7,10,9,6,6,8,10,5,5,9,8,7,7,9,10,5,5,8,9,6,6,10,)
const IS13_2 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS13_3 = (5,9,6,8,7,10,5,9,6,8,7,10,6,10,5,8,7,9,6,10,5,8,7,9,)

const IS14_1 = (8,7,7,10,9,6,6,8,10,5,5,9,8,7,7,9,10,5,5,8,9,6,6,10,)
const IS14_2 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS14_3 = (9,10,5,6,10,8,7,5,8,9,6,7,10,9,6,5,9,8,7,6,8,10,5,7,)

const IS15_1 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS15_2 = (2,1,1,3,2,1,1,2,3,1,1,2,2,1,1,2,3,1,1,2,2,1,1,3,)
const IS15_3 = (3,4,4,4,4,3,3,3,4,2,2,4,3,4,4,4,4,2,2,3,4,3,3,4,)
const IS15_4 = (5,9,6,8,7,10,5,9,6,8,7,10,6,10,5,8,7,9,6,10,5,8,7,9,)

const IS16_1 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS16_2 = (2,4,1,3,4,3,1,2,3,2,1,4,3,4,1,2,4,2,1,3,2,3,1,4,)
const IS16_3 = (8,6,7,7,8,5,6,6,9,5,5,5,8,5,7,7,8,5,5,5,9,6,6,6,)
const IS16_4 = (10,7,9,10,9,6,10,8,10,7,8,9,9,7,10,9,10,6,9,8,10,7,8,10,)

const IS17_1 = (5,9,6,8,7,10,)
const IS17_2 = (1,2,1,2,1,3,)
const IS17_3 = (2,4,3,3,4,4,)
const IS17_4 = (10,6,9,7,8,5,)

const IS18_1 = (1,2,1,2,1,3,1,2,1,2,1,3,1,3,1,2,1,2,1,3,1,2,1,2,)
const IS18_2 = (2,4,3,3,4,4,2,4,3,3,4,4,3,4,2,3,4,4,3,4,2,3,4,4,)
const IS18_3 = (5,5,6,7,5,5,5,6,6,6,7,5,5,5,5,7,6,6,6,5,5,5,7,6,)
const IS18_4 = (6,6,9,8,7,7,8,9,7,7,8,8,6,6,9,8,7,7,8,9,7,7,8,8,)
const IS18_5 = (10,9,10,9,8,10,10,10,9,8,9,10,9,10,10,10,8,9,9,10,10,8,10,9,)

const IS19_1 = (9,10,5,6,10,8,7,5,8,9,6,7,10,9,6,5,9,8,7,6,8,10,5,7,)
const IS19_2 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS19_3 = (5,5,6,8,5,7,5,9,6,6,7,8,5,6,5,8,6,7,6,9,5,5,7,8,)
const IS19_4 = (6,9,10,9,7,10,8,10,7,8,9,10,6,10,9,10,7,9,8,10,7,8,10,9,)

const IS20_1 = (10,6,9,7,8,5,10,6,9,7,8,5,9,5,10,7,8,6,9,5,10,7,8,6,)
const IS20_2 = (1,2,3,2,1,4,2,4,1,3,4,3,1,3,2,3,1,4,3,4,1,2,4,2,)
const IS20_3 = (5,7,6,8,7,6,5,8,6,5,5,9,6,7,5,8,7,5,5,8,5,6,6,9,)
const IS20_4 = (8,9,7,10,9,10,6,9,10,8,7,10,8,10,7,9,10,9,6,10,9,8,7,10,)

const IS21_1 = (5,6,6,7,7,5,5,6,7,5,7,6,)
const IS21_2 = (10,9,9,8,8,10,10,9,8,10,8,9,)
const IS21_3 = (6,5,10,9,5,7,8,7,6,9,10,8,)

# Primary invariants for BA_5B
 # : definitions of the types at the beginning of the file
const P1 = Val((P1_1,))
const P2 = Val((P2_1,))
const P3 = Val((P3_1,))
const P4 = Val((P4_1,))
const P5 = Val((P5_1,P5_2,))
const P6 = Val((P6_1,))
const P7 = Val((P7_1,P7_2,P7_3,))
const P8 = Val((P8_1,))
const P9 = Val((P9_1,))
const P10 = Val((P10_1,P10_2))

# Irreducible secondaries for group BA_5B
 # : definitions of the types at the beginning of the file
const IS1 = Val((IS1_1,IS1_2,))
const IS2 = Val((IS2_1,IS2_2,IS2_3,))
const IS3 = Val((IS3_1,IS3_2,IS3_3,))
const IS4 = Val((IS4_1,IS4_2,))
const IS5 = Val((IS5_1,IS5_2,IS5_3,))
const IS6 = Val((IS6_1,IS6_2,))
const IS7 = Val((IS7_1,IS7_2,IS7_3,))
const IS8 = Val((IS8_1,IS8_2,IS8_3,))
const IS9 = Val((IS9_1,IS9_2,IS9_3,IS9_4,))
const IS10 = Val((IS10_1,IS10_2,IS10_3,IS10_4,))
const IS11 = Val((IS11_1,IS11_2,IS11_3,))
const IS12 = Val((IS12_1,IS12_2,IS12_3,))
const IS13 = Val((IS13_1,IS13_2,IS13_3,))
const IS14 = Val((IS14_1,IS14_2,IS14_3,))
const IS15 = Val((IS15_1,IS15_2,IS15_3,IS15_4,))
const IS16 = Val((IS16_1,IS16_2,IS16_3,IS16_4,))
const IS17 = Val((IS17_1,IS17_2,IS17_3,IS17_4,))
const IS18 = Val((IS18_1,IS18_2,IS18_3,IS18_4,IS18_5,))
const IS19 = Val((IS19_1,IS19_2,IS19_3,IS19_4,))
const IS20 = Val((IS20_1,IS20_2,IS20_3,IS20_4,))
const IS21 = Val((IS21_1,IS21_2,IS21_3,))

function invariants(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for BA_5B
 # : what goes in the function for the evaluation
P1 = fpoly((x1,) , BA_5B.P1)
P2 = fpoly((x1,) , BA_5B.P2)
P3 = fpoly((x2,) , BA_5B.P3)
P4 = fpoly((x2,) , BA_5B.P4)
P5 = fpoly((x1,x1,) , BA_5B.P5)
P6 = fpoly((x3,) , BA_5B.P6)
P7 = fpoly((x1,x1,x1,) , BA_5B.P7)
P8 = fpoly((x3,) , BA_5B.P8)
P9 = fpoly((x4,) , BA_5B.P9)
P10 = fpoly((x3,x1) , BA_5B.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for group BA_5B
 # : what goes in the function for the evaluation
IS1 = fpoly((x1,x1,) , BA_5B.IS1)
IS2 = fpoly((x1,x1,x1,) , BA_5B.IS2)
IS3 = fpoly((x1,x1,x1,) , BA_5B.IS3)
IS4 = fpoly((x2,x1,) , BA_5B.IS4)
IS5 = fpoly((x1,x1,x1,) , BA_5B.IS5)
IS6 = fpoly((x2,x1,) , BA_5B.IS6)
IS7 = fpoly((x2,x1,x1,) , BA_5B.IS7)
IS8 = fpoly((x2,x1,x1,) , BA_5B.IS8)
IS9 = fpoly((x1,x1,x1,x1,) , BA_5B.IS9)
IS10 = fpoly((x1,x1,x1,x1,) , BA_5B.IS10)
IS11 = fpoly((x2,x1,x1,) , BA_5B.IS11)
IS12 = fpoly((x2,x1,x1,) , BA_5B.IS12)
IS13 = fpoly((x2,x1,x1,) , BA_5B.IS13)
IS14 = fpoly((x2,x1,x1,) , BA_5B.IS14)
IS15 = fpoly((x2,x1,x1,x1,) , BA_5B.IS15)
IS16 = fpoly((x2,x1,x1,x1,) , BA_5B.IS16)
IS17 = fpoly((x2,x1,x1,x1,) , BA_5B.IS17)
IS18 = fpoly((x1,x1,x1,x1,x1,) , BA_5B.IS18)
IS19 = fpoly((x2,x1,x1,x1,) , BA_5B.IS19)
IS20 = fpoly((x2,x1,x1,x1,) , BA_5B.IS20)
IS21 = fpoly((x2,x2,x1,) , BA_5B.IS21)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

SEC1  = 1
SEC2  = IS1
SEC3  = IS2
SEC4  = IS3
SEC5  = IS4
SEC6  = IS5
SEC7  = IS6
SEC8  = IS1^2
SEC9  = IS7
SEC10  = IS8
SEC11  = IS9
SEC12  = IS10
SEC13  = IS11
SEC14  = IS12
SEC15  = IS13
SEC16  = IS14
SEC17  = IS1*IS5
SEC18  = IS1*IS4
SEC19  = IS1*IS3
SEC20  = IS1*IS6
SEC21  = IS1*IS2
SEC22  = IS15
SEC23  = IS16
SEC24  = IS17
SEC25  = IS18
SEC26  = IS19
SEC27  = IS20
SEC28  = IS21
SEC29  = IS4^2
SEC30  = IS3*IS6
SEC31  = IS5*IS6
SEC32  = IS1*IS7
SEC33  = IS4*IS6
SEC34  = IS2*IS6
SEC35  = IS1*IS9
SEC36  = IS6^2
SEC37  = IS1*IS8
SEC38  = IS1*IS11
SEC39  = IS1^3
SEC40  = IS3^2
SEC41  = IS3*IS5
SEC42  = IS1*IS12
SEC43  = IS1*IS14
SEC44  = IS3*IS4
SEC45  = IS2*IS3
SEC46  = IS1*IS13
SEC47  = IS5*IS7
SEC48  = IS6*IS13
SEC49  = IS4*IS13
SEC50  = IS5*IS13
SEC51  = IS5*IS10
SEC52  = IS3*IS11
SEC53  = IS3*IS8
SEC54  = IS3*IS9
SEC55  = IS1^2*IS3
SEC56  = IS6*IS11
SEC57  = IS6*IS8
SEC58  = IS2*IS11
SEC59  = IS1*IS18
SEC60  = IS4*IS11
SEC61  = IS1*IS16
SEC62  = IS6*IS9
SEC63  = IS3*IS12
SEC64  = IS1*IS21
SEC65  = IS6*IS12
SEC66  = IS6*IS14
SEC67  = IS1*IS20
SEC68  = IS12*IS13
SEC69  = IS10*IS14
SEC70  = IS6*IS20
SEC71  = IS10^2
SEC72  = IS12*IS14
SEC73  = IS4*IS17
SEC74  = IS3*IS21
SEC75  = IS2*IS17
SEC76  = IS7*IS10
SEC77  = IS9*IS14
SEC78  = IS1*IS3*IS5
SEC79  = IS7*IS12
SEC80  = IS11*IS13
SEC81  = IS1^2*IS13
SEC82  = IS6*IS17
SEC83  = IS7*IS9
SEC84  = IS5*IS15
SEC85  = IS1*IS3*IS4
SEC86  = IS9*IS10
SEC87  = IS5*IS21
SEC88  = IS11^2
SEC89  = IS9^2
SEC90  = IS4*IS21
SEC91  = IS6*IS21
SEC92  = IS8*IS17
SEC93  = IS14*IS18
SEC94  = IS1*IS4*IS7
SEC95  = IS1*IS3*IS11
SEC96  = IS1*IS6*IS14
SEC97  = IS1^3*IS3
SEC98  = IS7*IS19
SEC99  = IS3*IS6^2
SEC100  = IS7*IS18
SEC101  = IS10*IS21
SEC102  = IS1*IS3*IS8
SEC103  = IS11*IS15
SEC104  = IS1*IS6*IS7
SEC105  = IS1*IS5*IS10
SEC106  = IS1^2*IS15
SEC107  = IS12*IS21
SEC108  = IS1*IS5*IS12
SEC109  = IS13*IS20
SEC110  = IS8*IS15
SEC111  = IS9*IS21
SEC112  = IS4^2*IS5
SEC113  = IS1*IS4*IS10
SEC114  = IS10*IS16
SEC115  = IS1^2*IS21
SEC116  = IS1*IS6*IS12
SEC117  = IS1*IS5*IS11
SEC118  = IS1*IS9*IS13
SEC119  = IS1*IS10^2
SEC120  = IS16*IS18
SEC121  = IS4*IS6*IS10
SEC122  = IS4*IS6*IS11
SEC123  = IS1^3*IS8
SEC124  = IS1*IS6*IS20
SEC125  = IS15*IS20
SEC126  = IS1*IS3*IS21
SEC127  = IS5*IS6*IS12
SEC128  = IS1*IS7*IS9
SEC129  = IS17*IS18
SEC130  = IS2*IS4*IS11
SEC131  = IS1*IS2*IS20
SEC132  = IS4^2*IS7
SEC133  = IS2*IS3*IS17
SEC134  = IS1^2*IS6*IS12
SEC135  = IS5*IS8*IS13
SEC136  = IS4^2*IS21
SEC137  = IS1*IS9*IS19
SEC138  = IS3*IS4*IS17
SEC139  = IS1^3*IS17
SEC140  = IS1^3*IS21
SEC141  = IS5*IS11*IS21
SEC142  = IS1*IS4*IS6*IS7
SEC143  = IS2*IS11*IS18
SEC144  = IS5*IS13*IS20

return (@SVector [P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,]), (@SVector [SEC1,SEC2,SEC3,SEC4,SEC5,SEC6,SEC7,SEC8,SEC9,SEC10,SEC11,SEC12,SEC13,SEC14,SEC15,SEC16,SEC17,SEC18,SEC19,SEC20,SEC21,SEC22,SEC23,SEC24,SEC25,SEC26,SEC27,SEC28,SEC29,SEC30,SEC31,SEC32,SEC33,SEC34,SEC35,SEC36,SEC37,SEC38,SEC39,SEC40,SEC41,SEC42,SEC43,SEC44,SEC45,SEC46,SEC47,SEC48,SEC49,SEC50,SEC51,SEC52,SEC53,SEC54,SEC55,SEC56,SEC57,SEC58,SEC59,SEC60,SEC61,SEC62,SEC63,SEC64,SEC65,SEC66,SEC67,SEC68,SEC69,SEC70,SEC71,SEC72,SEC73,SEC74,SEC75,SEC76,SEC77,SEC78,SEC79,SEC80,SEC81,SEC82,SEC83,SEC84,SEC85,SEC86,SEC87,SEC88,SEC89,SEC90,SEC91,SEC92,SEC93,SEC94,SEC95,SEC96,SEC97,SEC98,SEC99,SEC100,SEC101,SEC102,SEC103,SEC104,SEC105,SEC106,SEC107,SEC108,SEC109,SEC110,SEC111,SEC112,SEC113,SEC114,SEC115,SEC116,SEC117,SEC118,SEC119,SEC120,SEC121,SEC122,SEC123,SEC124,SEC125,SEC126,SEC127,SEC128,SEC129,SEC130,SEC131,SEC132,SEC133,SEC134,SEC135,SEC136,SEC137,SEC138,SEC139,SEC140,SEC141,SEC142,SEC143,SEC144,])
 end

function invariants_d(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1

   dx1 = @SVector ones(10)
   dx2 = 2 * x1
   dx3 = 3 * x2
   dx4 = 4 * x3
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for BA_5B
 # : what goes in the function for the derivatives
dP1 = fpoly_d((x1,),(dx1,) , BA_5B.P1)
dP2 = fpoly_d((x1,),(dx1,) , BA_5B.P2)
dP3 = fpoly_d((x2,),(dx2,) , BA_5B.P3)
dP4 = fpoly_d((x2,),(dx2,) , BA_5B.P4)
dP5 = fpoly_d((x1,x1,),(dx1,dx1,) , BA_5B.P5)
dP6 = fpoly_d((x3,),(dx3,) , BA_5B.P6)
dP7 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.P7)
dP8 = fpoly_d((x3,),(dx3,) , BA_5B.P8)
dP9 = fpoly_d((x4,),(dx4,) , BA_5B.P9)
dP10 = fpoly_d((x3,x1,),(dx3,dx1,) , BA_5B.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for group BA_5B
 # : what goes in the function for the evaluation
IS1 = fpoly((x1,x1,) , BA_5B.IS1)
IS2 = fpoly((x1,x1,x1,) , BA_5B.IS2)
IS3 = fpoly((x1,x1,x1,) , BA_5B.IS3)
IS4 = fpoly((x2,x1,) , BA_5B.IS4)
IS5 = fpoly((x1,x1,x1,) , BA_5B.IS5)
IS6 = fpoly((x2,x1,) , BA_5B.IS6)
IS7 = fpoly((x2,x1,x1,) , BA_5B.IS7)
IS8 = fpoly((x2,x1,x1,) , BA_5B.IS8)
IS9 = fpoly((x1,x1,x1,x1,) , BA_5B.IS9)
IS10 = fpoly((x1,x1,x1,x1,) , BA_5B.IS10)
IS11 = fpoly((x2,x1,x1,) , BA_5B.IS11)
IS12 = fpoly((x2,x1,x1,) , BA_5B.IS12)
IS13 = fpoly((x2,x1,x1,) , BA_5B.IS13)
IS14 = fpoly((x2,x1,x1,) , BA_5B.IS14)
IS15 = fpoly((x2,x1,x1,x1,) , BA_5B.IS15)
IS16 = fpoly((x2,x1,x1,x1,) , BA_5B.IS16)
IS17 = fpoly((x2,x1,x1,x1,) , BA_5B.IS17)
IS18 = fpoly((x1,x1,x1,x1,x1,) , BA_5B.IS18)
IS19 = fpoly((x2,x1,x1,x1,) , BA_5B.IS19)
IS20 = fpoly((x2,x1,x1,x1,) , BA_5B.IS20)
IS21 = fpoly((x2,x2,x1,) , BA_5B.IS21)

# Irreducible secondaries for group BA_5B
 # : what goes in the function for the derivatives
dIS1 = fpoly_d((x1,x1,),(dx1,dx1,) , BA_5B.IS1)
dIS2 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS2)
dIS3 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS3)
dIS4 = fpoly_d((x2,x1,),(dx2,dx1,) , BA_5B.IS4)
dIS5 = fpoly_d((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS5)
dIS6 = fpoly_d((x2,x1,),(dx2,dx1,) , BA_5B.IS6)
dIS7 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS7)
dIS8 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS8)
dIS9 = fpoly_d((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , BA_5B.IS9)
dIS10 = fpoly_d((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , BA_5B.IS10)
dIS11 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS11)
dIS12 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS12)
dIS13 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS13)
dIS14 = fpoly_d((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS14)
dIS15 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS15)
dIS16 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS16)
dIS17 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS17)
dIS18 = fpoly_d((x1,x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,dx1,) , BA_5B.IS18)
dIS19 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS19)
dIS20 = fpoly_d((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS20)
dIS21 = fpoly_d((x2,x2,x1,),(dx2,dx2,dx1,) , BA_5B.IS21)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

dSEC1   = @SVector zeros(10)
dSEC2  = dIS1
dSEC3  = dIS2
dSEC4  = dIS3
dSEC5  = dIS4
dSEC6  = dIS5
dSEC7  = dIS6
dSEC8  =  + dIS1*2IS1
dSEC9  = dIS7
dSEC10  = dIS8
dSEC11  = dIS9
dSEC12  = dIS10
dSEC13  = dIS11
dSEC14  = dIS12
dSEC15  = dIS13
dSEC16  = dIS14
dSEC17  =  + dIS1*IS5 + dIS5*IS1
dSEC18  =  + dIS1*IS4 + dIS4*IS1
dSEC19  =  + dIS1*IS3 + dIS3*IS1
dSEC20  =  + dIS1*IS6 + dIS6*IS1
dSEC21  =  + dIS1*IS2 + dIS2*IS1
dSEC22  = dIS15
dSEC23  = dIS16
dSEC24  = dIS17
dSEC25  = dIS18
dSEC26  = dIS19
dSEC27  = dIS20
dSEC28  = dIS21
dSEC29  =  + dIS4*2IS4
dSEC30  =  + dIS3*IS6 + dIS6*IS3
dSEC31  =  + dIS5*IS6 + dIS6*IS5
dSEC32  =  + dIS1*IS7 + dIS7*IS1
dSEC33  =  + dIS4*IS6 + dIS6*IS4
dSEC34  =  + dIS2*IS6 + dIS6*IS2
dSEC35  =  + dIS1*IS9 + dIS9*IS1
dSEC36  =  + dIS6*2IS6
dSEC37  =  + dIS1*IS8 + dIS8*IS1
dSEC38  =  + dIS1*IS11 + dIS11*IS1
dSEC39  =  + dIS1*3 * IS1 ^ 2
dSEC40  =  + dIS3*2IS3
dSEC41  =  + dIS3*IS5 + dIS5*IS3
dSEC42  =  + dIS1*IS12 + dIS12*IS1
dSEC43  =  + dIS1*IS14 + dIS14*IS1
dSEC44  =  + dIS3*IS4 + dIS4*IS3
dSEC45  =  + dIS2*IS3 + dIS3*IS2
dSEC46  =  + dIS1*IS13 + dIS13*IS1
dSEC47  =  + dIS5*IS7 + dIS7*IS5
dSEC48  =  + dIS6*IS13 + dIS13*IS6
dSEC49  =  + dIS4*IS13 + dIS13*IS4
dSEC50  =  + dIS5*IS13 + dIS13*IS5
dSEC51  =  + dIS5*IS10 + dIS10*IS5
dSEC52  =  + dIS3*IS11 + dIS11*IS3
dSEC53  =  + dIS3*IS8 + dIS8*IS3
dSEC54  =  + dIS3*IS9 + dIS9*IS3
dSEC55  =  + dIS1*(2IS1) * IS3 + dIS3*IS1 ^ 2
dSEC56  =  + dIS6*IS11 + dIS11*IS6
dSEC57  =  + dIS6*IS8 + dIS8*IS6
dSEC58  =  + dIS2*IS11 + dIS11*IS2
dSEC59  =  + dIS1*IS18 + dIS18*IS1
dSEC60  =  + dIS4*IS11 + dIS11*IS4
dSEC61  =  + dIS1*IS16 + dIS16*IS1
dSEC62  =  + dIS6*IS9 + dIS9*IS6
dSEC63  =  + dIS3*IS12 + dIS12*IS3
dSEC64  =  + dIS1*IS21 + dIS21*IS1
dSEC65  =  + dIS6*IS12 + dIS12*IS6
dSEC66  =  + dIS6*IS14 + dIS14*IS6
dSEC67  =  + dIS1*IS20 + dIS20*IS1
dSEC68  =  + dIS12*IS13 + dIS13*IS12
dSEC69  =  + dIS10*IS14 + dIS14*IS10
dSEC70  =  + dIS6*IS20 + dIS20*IS6
dSEC71  =  + dIS10*2IS10
dSEC72  =  + dIS12*IS14 + dIS14*IS12
dSEC73  =  + dIS4*IS17 + dIS17*IS4
dSEC74  =  + dIS3*IS21 + dIS21*IS3
dSEC75  =  + dIS2*IS17 + dIS17*IS2
dSEC76  =  + dIS7*IS10 + dIS10*IS7
dSEC77  =  + dIS9*IS14 + dIS14*IS9
dSEC78  =  + dIS1*IS3 * IS5 + dIS3*IS1 * IS5 + dIS5*IS1 * IS3
dSEC79  =  + dIS7*IS12 + dIS12*IS7
dSEC80  =  + dIS11*IS13 + dIS13*IS11
dSEC81  =  + dIS1*(2IS1) * IS13 + dIS13*IS1 ^ 2
dSEC82  =  + dIS6*IS17 + dIS17*IS6
dSEC83  =  + dIS7*IS9 + dIS9*IS7
dSEC84  =  + dIS5*IS15 + dIS15*IS5
dSEC85  =  + dIS1*IS3 * IS4 + dIS3*IS1 * IS4 + dIS4*IS1 * IS3
dSEC86  =  + dIS9*IS10 + dIS10*IS9
dSEC87  =  + dIS5*IS21 + dIS21*IS5
dSEC88  =  + dIS11*2IS11
dSEC89  =  + dIS9*2IS9
dSEC90  =  + dIS4*IS21 + dIS21*IS4
dSEC91  =  + dIS6*IS21 + dIS21*IS6
dSEC92  =  + dIS8*IS17 + dIS17*IS8
dSEC93  =  + dIS14*IS18 + dIS18*IS14
dSEC94  =  + dIS1*IS4 * IS7 + dIS4*IS1 * IS7 + dIS7*IS1 * IS4
dSEC95  =  + dIS1*IS3 * IS11 + dIS3*IS1 * IS11 + dIS11*IS1 * IS3
dSEC96  =  + dIS1*IS6 * IS14 + dIS6*IS1 * IS14 + dIS14*IS1 * IS6
dSEC97  =  + dIS1*(3 * IS1 ^ 2) * IS3 + dIS3*IS1 ^ 3
dSEC98  =  + dIS7*IS19 + dIS19*IS7
dSEC99  =  + dIS3*IS6 ^ 2 + dIS6*IS3 * (2IS6)
dSEC100  =  + dIS7*IS18 + dIS18*IS7
dSEC101  =  + dIS10*IS21 + dIS21*IS10
dSEC102  =  + dIS1*IS3 * IS8 + dIS3*IS1 * IS8 + dIS8*IS1 * IS3
dSEC103  =  + dIS11*IS15 + dIS15*IS11
dSEC104  =  + dIS1*IS6 * IS7 + dIS6*IS1 * IS7 + dIS7*IS1 * IS6
dSEC105  =  + dIS1*IS5 * IS10 + dIS5*IS1 * IS10 + dIS10*IS1 * IS5
dSEC106  =  + dIS1*(2IS1) * IS15 + dIS15*IS1 ^ 2
dSEC107  =  + dIS12*IS21 + dIS21*IS12
dSEC108  =  + dIS1*IS5 * IS12 + dIS5*IS1 * IS12 + dIS12*IS1 * IS5
dSEC109  =  + dIS13*IS20 + dIS20*IS13
dSEC110  =  + dIS8*IS15 + dIS15*IS8
dSEC111  =  + dIS9*IS21 + dIS21*IS9
dSEC112  =  + dIS4*(2IS4) * IS5 + dIS5*IS4 ^ 2
dSEC113  =  + dIS1*IS4 * IS10 + dIS4*IS1 * IS10 + dIS10*IS1 * IS4
dSEC114  =  + dIS10*IS16 + dIS16*IS10
dSEC115  =  + dIS1*(2IS1) * IS21 + dIS21*IS1 ^ 2
dSEC116  =  + dIS1*IS6 * IS12 + dIS6*IS1 * IS12 + dIS12*IS1 * IS6
dSEC117  =  + dIS1*IS5 * IS11 + dIS5*IS1 * IS11 + dIS11*IS1 * IS5
dSEC118  =  + dIS1*IS9 * IS13 + dIS9*IS1 * IS13 + dIS13*IS1 * IS9
dSEC119  =  + dIS1*IS10 ^ 2 + dIS10*IS1 * (2IS10)
dSEC120  =  + dIS16*IS18 + dIS18*IS16
dSEC121  =  + dIS4*IS6 * IS10 + dIS6*IS4 * IS10 + dIS10*IS4 * IS6
dSEC122  =  + dIS4*IS6 * IS11 + dIS6*IS4 * IS11 + dIS11*IS4 * IS6
dSEC123  =  + dIS1*(3 * IS1 ^ 2) * IS8 + dIS8*IS1 ^ 3
dSEC124  =  + dIS1*IS6 * IS20 + dIS6*IS1 * IS20 + dIS20*IS1 * IS6
dSEC125  =  + dIS15*IS20 + dIS20*IS15
dSEC126  =  + dIS1*IS3 * IS21 + dIS3*IS1 * IS21 + dIS21*IS1 * IS3
dSEC127  =  + dIS5*IS6 * IS12 + dIS6*IS5 * IS12 + dIS12*IS5 * IS6
dSEC128  =  + dIS1*IS7 * IS9 + dIS7*IS1 * IS9 + dIS9*IS1 * IS7
dSEC129  =  + dIS17*IS18 + dIS18*IS17
dSEC130  =  + dIS2*IS4 * IS11 + dIS4*IS2 * IS11 + dIS11*IS2 * IS4
dSEC131  =  + dIS1*IS2 * IS20 + dIS2*IS1 * IS20 + dIS20*IS1 * IS2
dSEC132  =  + dIS4*(2IS4) * IS7 + dIS7*IS4 ^ 2
dSEC133  =  + dIS2*IS3 * IS17 + dIS3*IS2 * IS17 + dIS17*IS2 * IS3
dSEC134  =  + dIS1*(2IS1) * IS6 * IS12 + dIS6*IS1 ^ 2 * IS12 + dIS12*IS1 ^ 2 * IS6
dSEC135  =  + dIS5*IS8 * IS13 + dIS8*IS5 * IS13 + dIS13*IS5 * IS8
dSEC136  =  + dIS4*(2IS4) * IS21 + dIS21*IS4 ^ 2
dSEC137  =  + dIS1*IS9 * IS19 + dIS9*IS1 * IS19 + dIS19*IS1 * IS9
dSEC138  =  + dIS3*IS4 * IS17 + dIS4*IS3 * IS17 + dIS17*IS3 * IS4
dSEC139  =  + dIS1*(3 * IS1 ^ 2) * IS17 + dIS17*IS1 ^ 3
dSEC140  =  + dIS1*(3 * IS1 ^ 2) * IS21 + dIS21*IS1 ^ 3
dSEC141  =  + dIS5*IS11 * IS21 + dIS11*IS5 * IS21 + dIS21*IS5 * IS11
dSEC142  =  + dIS1*IS4 * IS6 * IS7 + dIS4*IS1 * IS6 * IS7 + dIS6*IS1 * IS4 * IS7 + dIS7*IS1 * IS4 * IS6
dSEC143  =  + dIS2*IS11 * IS18 + dIS11*IS2 * IS18 + dIS18*IS2 * IS11
dSEC144  =  + dIS5*IS13 * IS20 + dIS13*IS5 * IS20 + dIS20*IS5 * IS13

return (dP1,dP2,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,), (dSEC1,dSEC2,dSEC3,dSEC4,dSEC5,dSEC6,dSEC7,dSEC8,dSEC9,dSEC10,dSEC11,dSEC12,dSEC13,dSEC14,dSEC15,dSEC16,dSEC17,dSEC18,dSEC19,dSEC20,dSEC21,dSEC22,dSEC23,dSEC24,dSEC25,dSEC26,dSEC27,dSEC28,dSEC29,dSEC30,dSEC31,dSEC32,dSEC33,dSEC34,dSEC35,dSEC36,dSEC37,dSEC38,dSEC39,dSEC40,dSEC41,dSEC42,dSEC43,dSEC44,dSEC45,dSEC46,dSEC47,dSEC48,dSEC49,dSEC50,dSEC51,dSEC52,dSEC53,dSEC54,dSEC55,dSEC56,dSEC57,dSEC58,dSEC59,dSEC60,dSEC61,dSEC62,dSEC63,dSEC64,dSEC65,dSEC66,dSEC67,dSEC68,dSEC69,dSEC70,dSEC71,dSEC72,dSEC73,dSEC74,dSEC75,dSEC76,dSEC77,dSEC78,dSEC79,dSEC80,dSEC81,dSEC82,dSEC83,dSEC84,dSEC85,dSEC86,dSEC87,dSEC88,dSEC89,dSEC90,dSEC91,dSEC92,dSEC93,dSEC94,dSEC95,dSEC96,dSEC97,dSEC98,dSEC99,dSEC100,dSEC101,dSEC102,dSEC103,dSEC104,dSEC105,dSEC106,dSEC107,dSEC108,dSEC109,dSEC110,dSEC111,dSEC112,dSEC113,dSEC114,dSEC115,dSEC116,dSEC117,dSEC118,dSEC119,dSEC120,dSEC121,dSEC122,dSEC123,dSEC124,dSEC125,dSEC126,dSEC127,dSEC128,dSEC129,dSEC130,dSEC131,dSEC132,dSEC133,dSEC134,dSEC135,dSEC136,dSEC137,dSEC138,dSEC139,dSEC140,dSEC141,dSEC142,dSEC143,dSEC144,)
 end

function invariants_ed(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1

   dx1 = @SVector ones(10)
   dx2 = 2 * x1
   dx3 = 3 * x2
   dx4 = 4 * x3
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

# Primary invariants for BA_5B
 # : what goes in the function for the evaluation and derivatives
P1, dP1 = fpoly_ed((x1,),(dx1,) , BA_5B.P1)
P2, dP2 = fpoly_ed((x1,),(dx1,) , BA_5B.P2)
P3, dP3 = fpoly_ed((x2,),(dx2,) , BA_5B.P3)
P4, dP4 = fpoly_ed((x2,),(dx2,) , BA_5B.P4)
P5, dP5 = fpoly_ed((x1,x1,),(dx1,dx1,) , BA_5B.P5)
P6, dP6 = fpoly_ed((x3,),(dx3,) , BA_5B.P6)
P7, dP7 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.P7)
P8, dP8 = fpoly_ed((x3,),(dx3,) , BA_5B.P8)
P9, dP9 = fpoly_ed((x4,),(dx4,) , BA_5B.P9)
P10, dP10 = fpoly_ed((x3,x1),(dx3,dx1) , BA_5B.P10)

   #------------------------------------------------
   # Irreducible secondaries
   #------------------------------------------------

# Irreducible secondaries for group BA_5B
 # : what goes in the function for the evaluation and derivatives
IS1, dIS1 = fpoly_ed((x1,x1,),(dx1,dx1,) , BA_5B.IS1)
IS2, dIS2 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS2)
IS3, dIS3 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS3)
IS4, dIS4 = fpoly_ed((x2,x1,),(dx2,dx1,) , BA_5B.IS4)
IS5, dIS5 = fpoly_ed((x1,x1,x1,),(dx1,dx1,dx1,) , BA_5B.IS5)
IS6, dIS6 = fpoly_ed((x2,x1,),(dx2,dx1,) , BA_5B.IS6)
IS7, dIS7 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS7)
IS8, dIS8 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS8)
IS9, dIS9 = fpoly_ed((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , BA_5B.IS9)
IS10, dIS10 = fpoly_ed((x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,) , BA_5B.IS10)
IS11, dIS11 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS11)
IS12, dIS12 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS12)
IS13, dIS13 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS13)
IS14, dIS14 = fpoly_ed((x2,x1,x1,),(dx2,dx1,dx1,) , BA_5B.IS14)
IS15, dIS15 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS15)
IS16, dIS16 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS16)
IS17, dIS17 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS17)
IS18, dIS18 = fpoly_ed((x1,x1,x1,x1,x1,),(dx1,dx1,dx1,dx1,dx1,) , BA_5B.IS18)
IS19, dIS19 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS19)
IS20, dIS20 = fpoly_ed((x2,x1,x1,x1,),(dx2,dx1,dx1,dx1,) , BA_5B.IS20)
IS21, dIS21 = fpoly_ed((x2,x2,x1,),(dx2,dx2,dx1,) , BA_5B.IS21)

   #------------------------------------------------
   # All secondaries
   #------------------------------------------------

SEC1  = 1
SEC2  = IS1
SEC3  = IS2
SEC4  = IS3
SEC5  = IS4
SEC6  = IS5
SEC7  = IS6
SEC8  = IS1^2
SEC9  = IS7
SEC10  = IS8
SEC11  = IS9
SEC12  = IS10
SEC13  = IS11
SEC14  = IS12
SEC15  = IS13
SEC16  = IS14
SEC17  = IS1*IS5
SEC18  = IS1*IS4
SEC19  = IS1*IS3
SEC20  = IS1*IS6
SEC21  = IS1*IS2
SEC22  = IS15
SEC23  = IS16
SEC24  = IS17
SEC25  = IS18
SEC26  = IS19
SEC27  = IS20
SEC28  = IS21
SEC29  = IS4^2
SEC30  = IS3*IS6
SEC31  = IS5*IS6
SEC32  = IS1*IS7
SEC33  = IS4*IS6
SEC34  = IS2*IS6
SEC35  = IS1*IS9
SEC36  = IS6^2
SEC37  = IS1*IS8
SEC38  = IS1*IS11
SEC39  = IS1^3
SEC40  = IS3^2
SEC41  = IS3*IS5
SEC42  = IS1*IS12
SEC43  = IS1*IS14
SEC44  = IS3*IS4
SEC45  = IS2*IS3
SEC46  = IS1*IS13
SEC47  = IS5*IS7
SEC48  = IS6*IS13
SEC49  = IS4*IS13
SEC50  = IS5*IS13
SEC51  = IS5*IS10
SEC52  = IS3*IS11
SEC53  = IS3*IS8
SEC54  = IS3*IS9
SEC55  = IS1^2*IS3
SEC56  = IS6*IS11
SEC57  = IS6*IS8
SEC58  = IS2*IS11
SEC59  = IS1*IS18
SEC60  = IS4*IS11
SEC61  = IS1*IS16
SEC62  = IS6*IS9
SEC63  = IS3*IS12
SEC64  = IS1*IS21
SEC65  = IS6*IS12
SEC66  = IS6*IS14
SEC67  = IS1*IS20
SEC68  = IS12*IS13
SEC69  = IS10*IS14
SEC70  = IS6*IS20
SEC71  = IS10^2
SEC72  = IS12*IS14
SEC73  = IS4*IS17
SEC74  = IS3*IS21
SEC75  = IS2*IS17
SEC76  = IS7*IS10
SEC77  = IS9*IS14
SEC78  = IS1*IS3*IS5
SEC79  = IS7*IS12
SEC80  = IS11*IS13
SEC81  = IS1^2*IS13
SEC82  = IS6*IS17
SEC83  = IS7*IS9
SEC84  = IS5*IS15
SEC85  = IS1*IS3*IS4
SEC86  = IS9*IS10
SEC87  = IS5*IS21
SEC88  = IS11^2
SEC89  = IS9^2
SEC90  = IS4*IS21
SEC91  = IS6*IS21
SEC92  = IS8*IS17
SEC93  = IS14*IS18
SEC94  = IS1*IS4*IS7
SEC95  = IS1*IS3*IS11
SEC96  = IS1*IS6*IS14
SEC97  = IS1^3*IS3
SEC98  = IS7*IS19
SEC99  = IS3*IS6^2
SEC100  = IS7*IS18
SEC101  = IS10*IS21
SEC102  = IS1*IS3*IS8
SEC103  = IS11*IS15
SEC104  = IS1*IS6*IS7
SEC105  = IS1*IS5*IS10
SEC106  = IS1^2*IS15
SEC107  = IS12*IS21
SEC108  = IS1*IS5*IS12
SEC109  = IS13*IS20
SEC110  = IS8*IS15
SEC111  = IS9*IS21
SEC112  = IS4^2*IS5
SEC113  = IS1*IS4*IS10
SEC114  = IS10*IS16
SEC115  = IS1^2*IS21
SEC116  = IS1*IS6*IS12
SEC117  = IS1*IS5*IS11
SEC118  = IS1*IS9*IS13
SEC119  = IS1*IS10^2
SEC120  = IS16*IS18
SEC121  = IS4*IS6*IS10
SEC122  = IS4*IS6*IS11
SEC123  = IS1^3*IS8
SEC124  = IS1*IS6*IS20
SEC125  = IS15*IS20
SEC126  = IS1*IS3*IS21
SEC127  = IS5*IS6*IS12
SEC128  = IS1*IS7*IS9
SEC129  = IS17*IS18
SEC130  = IS2*IS4*IS11
SEC131  = IS1*IS2*IS20
SEC132  = IS4^2*IS7
SEC133  = IS2*IS3*IS17
SEC134  = IS1^2*IS6*IS12
SEC135  = IS5*IS8*IS13
SEC136  = IS4^2*IS21
SEC137  = IS1*IS9*IS19
SEC138  = IS3*IS4*IS17
SEC139  = IS1^3*IS17
SEC140  = IS1^3*IS21
SEC141  = IS5*IS11*IS21
SEC142  = IS1*IS4*IS6*IS7
SEC143  = IS2*IS11*IS18
SEC144  = IS5*IS13*IS20

dSEC1   = @SVector zeros(10)
dSEC2  = dIS1
dSEC3  = dIS2
dSEC4  = dIS3
dSEC5  = dIS4
dSEC6  = dIS5
dSEC7  = dIS6
dSEC8  =  + dIS1*2IS1
dSEC9  = dIS7
dSEC10  = dIS8
dSEC11  = dIS9
dSEC12  = dIS10
dSEC13  = dIS11
dSEC14  = dIS12
dSEC15  = dIS13
dSEC16  = dIS14
dSEC17  =  + dIS1*IS5 + dIS5*IS1
dSEC18  =  + dIS1*IS4 + dIS4*IS1
dSEC19  =  + dIS1*IS3 + dIS3*IS1
dSEC20  =  + dIS1*IS6 + dIS6*IS1
dSEC21  =  + dIS1*IS2 + dIS2*IS1
dSEC22  = dIS15
dSEC23  = dIS16
dSEC24  = dIS17
dSEC25  = dIS18
dSEC26  = dIS19
dSEC27  = dIS20
dSEC28  = dIS21
dSEC29  =  + dIS4*2IS4
dSEC30  =  + dIS3*IS6 + dIS6*IS3
dSEC31  =  + dIS5*IS6 + dIS6*IS5
dSEC32  =  + dIS1*IS7 + dIS7*IS1
dSEC33  =  + dIS4*IS6 + dIS6*IS4
dSEC34  =  + dIS2*IS6 + dIS6*IS2
dSEC35  =  + dIS1*IS9 + dIS9*IS1
dSEC36  =  + dIS6*2IS6
dSEC37  =  + dIS1*IS8 + dIS8*IS1
dSEC38  =  + dIS1*IS11 + dIS11*IS1
dSEC39  =  + dIS1*3 * IS1 ^ 2
dSEC40  =  + dIS3*2IS3
dSEC41  =  + dIS3*IS5 + dIS5*IS3
dSEC42  =  + dIS1*IS12 + dIS12*IS1
dSEC43  =  + dIS1*IS14 + dIS14*IS1
dSEC44  =  + dIS3*IS4 + dIS4*IS3
dSEC45  =  + dIS2*IS3 + dIS3*IS2
dSEC46  =  + dIS1*IS13 + dIS13*IS1
dSEC47  =  + dIS5*IS7 + dIS7*IS5
dSEC48  =  + dIS6*IS13 + dIS13*IS6
dSEC49  =  + dIS4*IS13 + dIS13*IS4
dSEC50  =  + dIS5*IS13 + dIS13*IS5
dSEC51  =  + dIS5*IS10 + dIS10*IS5
dSEC52  =  + dIS3*IS11 + dIS11*IS3
dSEC53  =  + dIS3*IS8 + dIS8*IS3
dSEC54  =  + dIS3*IS9 + dIS9*IS3
dSEC55  =  + dIS1*(2IS1) * IS3 + dIS3*IS1 ^ 2
dSEC56  =  + dIS6*IS11 + dIS11*IS6
dSEC57  =  + dIS6*IS8 + dIS8*IS6
dSEC58  =  + dIS2*IS11 + dIS11*IS2
dSEC59  =  + dIS1*IS18 + dIS18*IS1
dSEC60  =  + dIS4*IS11 + dIS11*IS4
dSEC61  =  + dIS1*IS16 + dIS16*IS1
dSEC62  =  + dIS6*IS9 + dIS9*IS6
dSEC63  =  + dIS3*IS12 + dIS12*IS3
dSEC64  =  + dIS1*IS21 + dIS21*IS1
dSEC65  =  + dIS6*IS12 + dIS12*IS6
dSEC66  =  + dIS6*IS14 + dIS14*IS6
dSEC67  =  + dIS1*IS20 + dIS20*IS1
dSEC68  =  + dIS12*IS13 + dIS13*IS12
dSEC69  =  + dIS10*IS14 + dIS14*IS10
dSEC70  =  + dIS6*IS20 + dIS20*IS6
dSEC71  =  + dIS10*2IS10
dSEC72  =  + dIS12*IS14 + dIS14*IS12
dSEC73  =  + dIS4*IS17 + dIS17*IS4
dSEC74  =  + dIS3*IS21 + dIS21*IS3
dSEC75  =  + dIS2*IS17 + dIS17*IS2
dSEC76  =  + dIS7*IS10 + dIS10*IS7
dSEC77  =  + dIS9*IS14 + dIS14*IS9
dSEC78  =  + dIS1*IS3 * IS5 + dIS3*IS1 * IS5 + dIS5*IS1 * IS3
dSEC79  =  + dIS7*IS12 + dIS12*IS7
dSEC80  =  + dIS11*IS13 + dIS13*IS11
dSEC81  =  + dIS1*(2IS1) * IS13 + dIS13*IS1 ^ 2
dSEC82  =  + dIS6*IS17 + dIS17*IS6
dSEC83  =  + dIS7*IS9 + dIS9*IS7
dSEC84  =  + dIS5*IS15 + dIS15*IS5
dSEC85  =  + dIS1*IS3 * IS4 + dIS3*IS1 * IS4 + dIS4*IS1 * IS3
dSEC86  =  + dIS9*IS10 + dIS10*IS9
dSEC87  =  + dIS5*IS21 + dIS21*IS5
dSEC88  =  + dIS11*2IS11
dSEC89  =  + dIS9*2IS9
dSEC90  =  + dIS4*IS21 + dIS21*IS4
dSEC91  =  + dIS6*IS21 + dIS21*IS6
dSEC92  =  + dIS8*IS17 + dIS17*IS8
dSEC93  =  + dIS14*IS18 + dIS18*IS14
dSEC94  =  + dIS1*IS4 * IS7 + dIS4*IS1 * IS7 + dIS7*IS1 * IS4
dSEC95  =  + dIS1*IS3 * IS11 + dIS3*IS1 * IS11 + dIS11*IS1 * IS3
dSEC96  =  + dIS1*IS6 * IS14 + dIS6*IS1 * IS14 + dIS14*IS1 * IS6
dSEC97  =  + dIS1*(3 * IS1 ^ 2) * IS3 + dIS3*IS1 ^ 3
dSEC98  =  + dIS7*IS19 + dIS19*IS7
dSEC99  =  + dIS3*IS6 ^ 2 + dIS6*IS3 * (2IS6)
dSEC100  =  + dIS7*IS18 + dIS18*IS7
dSEC101  =  + dIS10*IS21 + dIS21*IS10
dSEC102  =  + dIS1*IS3 * IS8 + dIS3*IS1 * IS8 + dIS8*IS1 * IS3
dSEC103  =  + dIS11*IS15 + dIS15*IS11
dSEC104  =  + dIS1*IS6 * IS7 + dIS6*IS1 * IS7 + dIS7*IS1 * IS6
dSEC105  =  + dIS1*IS5 * IS10 + dIS5*IS1 * IS10 + dIS10*IS1 * IS5
dSEC106  =  + dIS1*(2IS1) * IS15 + dIS15*IS1 ^ 2
dSEC107  =  + dIS12*IS21 + dIS21*IS12
dSEC108  =  + dIS1*IS5 * IS12 + dIS5*IS1 * IS12 + dIS12*IS1 * IS5
dSEC109  =  + dIS13*IS20 + dIS20*IS13
dSEC110  =  + dIS8*IS15 + dIS15*IS8
dSEC111  =  + dIS9*IS21 + dIS21*IS9
dSEC112  =  + dIS4*(2IS4) * IS5 + dIS5*IS4 ^ 2
dSEC113  =  + dIS1*IS4 * IS10 + dIS4*IS1 * IS10 + dIS10*IS1 * IS4
dSEC114  =  + dIS10*IS16 + dIS16*IS10
dSEC115  =  + dIS1*(2IS1) * IS21 + dIS21*IS1 ^ 2
dSEC116  =  + dIS1*IS6 * IS12 + dIS6*IS1 * IS12 + dIS12*IS1 * IS6
dSEC117  =  + dIS1*IS5 * IS11 + dIS5*IS1 * IS11 + dIS11*IS1 * IS5
dSEC118  =  + dIS1*IS9 * IS13 + dIS9*IS1 * IS13 + dIS13*IS1 * IS9
dSEC119  =  + dIS1*IS10 ^ 2 + dIS10*IS1 * (2IS10)
dSEC120  =  + dIS16*IS18 + dIS18*IS16
dSEC121  =  + dIS4*IS6 * IS10 + dIS6*IS4 * IS10 + dIS10*IS4 * IS6
dSEC122  =  + dIS4*IS6 * IS11 + dIS6*IS4 * IS11 + dIS11*IS4 * IS6
dSEC123  =  + dIS1*(3 * IS1 ^ 2) * IS8 + dIS8*IS1 ^ 3
dSEC124  =  + dIS1*IS6 * IS20 + dIS6*IS1 * IS20 + dIS20*IS1 * IS6
dSEC125  =  + dIS15*IS20 + dIS20*IS15
dSEC126  =  + dIS1*IS3 * IS21 + dIS3*IS1 * IS21 + dIS21*IS1 * IS3
dSEC127  =  + dIS5*IS6 * IS12 + dIS6*IS5 * IS12 + dIS12*IS5 * IS6
dSEC128  =  + dIS1*IS7 * IS9 + dIS7*IS1 * IS9 + dIS9*IS1 * IS7
dSEC129  =  + dIS17*IS18 + dIS18*IS17
dSEC130  =  + dIS2*IS4 * IS11 + dIS4*IS2 * IS11 + dIS11*IS2 * IS4
dSEC131  =  + dIS1*IS2 * IS20 + dIS2*IS1 * IS20 + dIS20*IS1 * IS2
dSEC132  =  + dIS4*(2IS4) * IS7 + dIS7*IS4 ^ 2
dSEC133  =  + dIS2*IS3 * IS17 + dIS3*IS2 * IS17 + dIS17*IS2 * IS3
dSEC134  =  + dIS1*(2IS1) * IS6 * IS12 + dIS6*IS1 ^ 2 * IS12 + dIS12*IS1 ^ 2 * IS6
dSEC135  =  + dIS5*IS8 * IS13 + dIS8*IS5 * IS13 + dIS13*IS5 * IS8
dSEC136  =  + dIS4*(2IS4) * IS21 + dIS21*IS4 ^ 2
dSEC137  =  + dIS1*IS9 * IS19 + dIS9*IS1 * IS19 + dIS19*IS1 * IS9
dSEC138  =  + dIS3*IS4 * IS17 + dIS4*IS3 * IS17 + dIS17*IS3 * IS4
dSEC139  =  + dIS1*(3 * IS1 ^ 2) * IS17 + dIS17*IS1 ^ 3
dSEC140  =  + dIS1*(3 * IS1 ^ 2) * IS21 + dIS21*IS1 ^ 3
dSEC141  =  + dIS5*IS11 * IS21 + dIS11*IS5 * IS21 + dIS21*IS5 * IS11
dSEC142  =  + dIS1*IS4 * IS6 * IS7 + dIS4*IS1 * IS6 * IS7 + dIS6*IS1 * IS4 * IS7 + dIS7*IS1 * IS4 * IS6
dSEC143  =  + dIS2*IS11 * IS18 + dIS11*IS2 * IS18 + dIS18*IS2 * IS11
dSEC144  =  + dIS5*IS13 * IS20 + dIS13*IS5 * IS20 + dIS20*IS5 * IS13

return (@SVector [P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,]), (@SVector [SEC1,SEC2,SEC3,SEC4,SEC5,SEC6,SEC7,SEC8,SEC9,SEC10,SEC11,SEC12,SEC13,SEC14,SEC15,SEC16,SEC17,SEC18,SEC19,SEC20,SEC21,SEC22,SEC23,SEC24,SEC25,SEC26,SEC27,SEC28,SEC29,SEC30,SEC31,SEC32,SEC33,SEC34,SEC35,SEC36,SEC37,SEC38,SEC39,SEC40,SEC41,SEC42,SEC43,SEC44,SEC45,SEC46,SEC47,SEC48,SEC49,SEC50,SEC51,SEC52,SEC53,SEC54,SEC55,SEC56,SEC57,SEC58,SEC59,SEC60,SEC61,SEC62,SEC63,SEC64,SEC65,SEC66,SEC67,SEC68,SEC69,SEC70,SEC71,SEC72,SEC73,SEC74,SEC75,SEC76,SEC77,SEC78,SEC79,SEC80,SEC81,SEC82,SEC83,SEC84,SEC85,SEC86,SEC87,SEC88,SEC89,SEC90,SEC91,SEC92,SEC93,SEC94,SEC95,SEC96,SEC97,SEC98,SEC99,SEC100,SEC101,SEC102,SEC103,SEC104,SEC105,SEC106,SEC107,SEC108,SEC109,SEC110,SEC111,SEC112,SEC113,SEC114,SEC115,SEC116,SEC117,SEC118,SEC119,SEC120,SEC121,SEC122,SEC123,SEC124,SEC125,SEC126,SEC127,SEC128,SEC129,SEC130,SEC131,SEC132,SEC133,SEC134,SEC135,SEC136,SEC137,SEC138,SEC139,SEC140,SEC141,SEC142,SEC143,SEC144,]), (@SVector [dP1,dP2,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,]), (@SVector [dSEC1,dSEC2,dSEC3,dSEC4,dSEC5,dSEC6,dSEC7,dSEC8,dSEC9,dSEC10,dSEC11,dSEC12,dSEC13,dSEC14,dSEC15,dSEC16,dSEC17,dSEC18,dSEC19,dSEC20,dSEC21,dSEC22,dSEC23,dSEC24,dSEC25,dSEC26,dSEC27,dSEC28,dSEC29,dSEC30,dSEC31,dSEC32,dSEC33,dSEC34,dSEC35,dSEC36,dSEC37,dSEC38,dSEC39,dSEC40,dSEC41,dSEC42,dSEC43,dSEC44,dSEC45,dSEC46,dSEC47,dSEC48,dSEC49,dSEC50,dSEC51,dSEC52,dSEC53,dSEC54,dSEC55,dSEC56,dSEC57,dSEC58,dSEC59,dSEC60,dSEC61,dSEC62,dSEC63,dSEC64,dSEC65,dSEC66,dSEC67,dSEC68,dSEC69,dSEC70,dSEC71,dSEC72,dSEC73,dSEC74,dSEC75,dSEC76,dSEC77,dSEC78,dSEC79,dSEC80,dSEC81,dSEC82,dSEC83,dSEC84,dSEC85,dSEC86,dSEC87,dSEC88,dSEC89,dSEC90,dSEC91,dSEC92,dSEC93,dSEC94,dSEC95,dSEC96,dSEC97,dSEC98,dSEC99,dSEC100,dSEC101,dSEC102,dSEC103,dSEC104,dSEC105,dSEC106,dSEC107,dSEC108,dSEC109,dSEC110,dSEC111,dSEC112,dSEC113,dSEC114,dSEC115,dSEC116,dSEC117,dSEC118,dSEC119,dSEC120,dSEC121,dSEC122,dSEC123,dSEC124,dSEC125,dSEC126,dSEC127,dSEC128,dSEC129,dSEC130,dSEC131,dSEC132,dSEC133,dSEC134,dSEC135,dSEC136,dSEC137,dSEC138,dSEC139,dSEC140,dSEC141,dSEC142,dSEC143,dSEC144,])
 end

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
struct PICO{T <: Real, N <: Int} <: AbstractMeltRate
    γT :: T                       #(calibrated) heat exchange velocity
    λ1 :: T                       #liquidus slope
    λ2 :: T                       #liquidus intercept
    λ3 :: T                       #liquidus pressure coefficient
    L :: T                        #latent heat of fusion for ice
    c :: T                        #specific heat capacity of ocean water
    βs :: T                       #haline contraction coefficient
    βt :: T                       #thermal expansion coefficient
    g  :: T                       #gravity
    ρi :: T                       #ice density
    ρw :: T                       #water density
    C  :: T                       #overturning strength
    nbox :: N                     #number of boxes
    S0 :: T                       #input salinity
    T0 :: T                       #input temperature
    ice_front_mask::Array{Bool,2} #specify ice front location
    use_box_mean_depth::Bool      #flag to specify whether to use the mean depth in a box for the calculation of pk
end

"""
    function PICO(; <kwargs>)

Construct a pico structure to prescribe the melt rate in WAVI

Keyword arguments
=================
- γT: calibrated heat exchange velocity (units m/s). Defaults to 1.44*1e-5 according to Favier 2019 for 10 boxes (10.5194/gmd-12-2255-2019)
- λ1: liquidus slope (units ∘C)
- λ2: liquidus intercept (units ∘C)
- λ3: liquidus pressure coefficient (units K/m)
- L: latent heat of fusion for ice (units: J/kg)
- c: specific heat capacity of water (units J/kg/K)
- βs: haline contraction coefficient (units psu^(-1))
- βt: thermal expansion coefficient (units ∘C^{-1})
- g: gravitational acceleration (units m / s^2)
- ρi: ice density (units kg/m^3)
- ρw: water density (units kg/m^3)
- C: overturning circulation parameter (units m^6 /s /kg)
- S0: input salinity in box 0  (units PSU)
- T0: input temperature in box 0 (units ∘C)
- ice_front_mask: mask specifying the location of the ice front
- use_box_mean_depth: flag to specify whether or not to use the mean depth of the draft in each box in the computation of pressure in box
"""
function PICO(;
                γT = 1.44*1e-5,
                λ1 = -5.73e-2,
                λ2 = 8.32e-4,
                λ3 = 7.61e-4,
                L = 3.35e5,
                c = 3.974e3,
                βs = 7.7e-4,
                βt = 7.5e-5,
                g = 9.81,
                ρi = 918.8,
                ρw = 1028.0,
                C = 1.0e6,
                nbox = 10,
                S0 = 34.6,
                T0 = 0.0,
                ice_front_mask = nothing,
                use_box_mean_depth = false)

    #check that we've input an ice front
    ~(ice_front_mask === nothing) || throw(ArgumentError("Must input an ice front matrix"))

    #map the ice front mask to Boolean
    try
        ice_front_mask = convert(Array{Bool,2}, ice_front_mask)
    catch
        throw(ArgumentError("Ice front mask must be able to be cast as a Boolean matrix"))
    end
    return PICO(γT, λ1, λ2, λ3, L, c, βs, βt, g, ρi, ρw, C, nbox, S0, T0, ice_front_mask, use_box_mean_depth)
end

"""
    function update_melt_rate!(pico::PICO, fields, grid)

Update the ice model melt rate for a PICO melt rate parametrization
"""
function update_melt_rate!(pico::PICO, fields, grid)
    @unpack basal_melt, h, b, grounded_fraction = fields.gh #get the ice thickness and grounded fraction

    #check that the ice front mask has the correct dimensions
    (size(pico.ice_front_mask)) == (grid.nx, grid.ny) || throw(ArgumentError("Ice front mask ($(size(ice_front_mask))) size must be same as model grid ($nx x $ny)"))

    #compute the ice draft
    zb = fields.gh.b .* (grounded_fraction .== 1) + - pico.ρi / pico.ρw .* fields.gh.h .* (grounded_fraction .< 1)

    set_pico_melt_rate!(basal_melt,
                        grounded_fraction,
                        pico,
                        grid,
                        zb)
    return nothing
end

"""
    function set_pico_melt_rate!(melt, grounded_fraction, pico, grid, zb)

Update the melt rate (stored in melt) using the specified pico model.
"""
function set_pico_melt_rate!(melt,grounded_fraction, pico, grid, zb)
    #determine which box certain parts are in (i.e. are you floating, and what is the distance to the ice front?)
    grounding_line_mask    = get_gl_mask(grounded_fraction) #boolean matrix specifying where the grounding line is
    boxes = get_pico_boxes(grounding_line_mask, pico.ice_front_mask, grounded_fraction, pico.nbox)

    #set the melt rate in boxes sequentially
    T_out, S_out,q = set_box_one_melt_rate!(melt, boxes,pico, grid.dx,grid.dy, zb)
    for k = 2:pico.nbox
        T_out, S_out = set_box_k_melt_rate!(melt,k, T_out, S_out,q, boxes,pico, grid.dx,grid.dy, zb)
    end
end

"""
    function get_pico_boxes()

Return an array of boxes for the PICO melt rate. Returns 0 for grounded points.
"""
function get_pico_boxes(grounding_line_mask, ice_front_mask, grounded_fraction, nbox)
    nx, ny = size(grounding_line_mask)
    boxes = zeros(nx,ny)
    r = zeros(nx,ny)

    inds_gl = Tuple.(findall(!iszero, grounding_line_mask))
    inds_if = Tuple.(findall(!iszero, ice_front_mask))     #indices of ice front and grounding line
    for i = 1:nx
        for j = 1:ny
            if grounded_fraction[i,j] == 0 & ~(grounding_line_mask[i,j]) #if we're in the shelf and not at a grounding line point
                #compute distance to the gl
                dgl = findmin(((i .- first.(inds_gl)).^2 + (j .- last.(inds_gl)).^2).^(1/2))[1]
                #compute distance to the ice front
                dif = findmin(((i .- first.(inds_if)).^2 + (j .- last.(inds_if)).^2).^(1/2))[1]
                #assign box
                r[i,j] = dgl / (dgl + dif) #equation 10 in Reese et al 2018
                boxes[i,j] = ceil(r[i,j] * nbox)
            end
         end
    end
    return boxes
end

"""
    function get_gl_mask(grounded_fraction)

Return a matrix with ones corresponding to locations of grounding line points
Grounding line points determines as those points with grounded fraction < 1, and adjacted to a point with grounded fraction = 1
"""
function get_gl_mask(grounded_fraction)
    nx, ny = size(grounded_fraction)
    gl_mask = falses(nx,ny)
    for i = 1:nx
        for j = 1:ny
            if (grounded_fraction[i,j] < 1)
                #compute the box of adjacent points
                adj_box_grounded = get_adj_box(i,j,grounded_fraction)
                #if any adjacent point has grounded_fraction == 1, tag this point as being in gl
                if any(adj_box_grounded .== 1)
                    gl_mask[i,j] = true
                end
            end
        end
    end
    return gl_mask
end

"""
    function get_adj_box(i,j,m)

Returns the entries of m which are adjacent to grid point indexed by i, j, taking into account edges of matrix
"""
function get_adj_box(i,j, m)
    nx, ny = size(m)
    lxn = min(1,i-1) #length of box negative x
    lxp = min(1,nx - i) #length of box positive x
    lyn = min(1,j-1)
    lyp = min(1,ny - j)
    box = m[((i-lxn):(i+lxp)), ((j - lyn):(j+lyp))]
    return box
end

"""
    function set_box_one_melt_rate!(melt, boxes, pico, dx, dy, zb)

Set the melt rate in pico box one (adjacent to the grounding line). Box location specified by boxes array.
Returns the mean temperature T1bar and salinity S1bar in box one (used as boundary conditions on box two), and the overturning flux q (used in all subsequent boxes)
"""
function set_box_one_melt_rate!(melt,boxes,pico, dx, dy, zb)
    m = zeros(size(zb))

    #scalar quantities
    A1 = sum(sum(boxes .== 1))*dx *dy
    g1 = A1 * pico.γT
    ν = pico.ρi / pico.ρw
    λ = pico.L / pico.c
    s  = pico.S0 / ν / λ
    if ~(pico.use_box_mean_depth) #don't use mean depth, below are vector quantities
        Tstar = pico.λ1 .* pico.S0 .+ pico.λ2 .+ pico.λ3 .* zb .- pico.T0
        c =  g1 .* Tstar / pico.C / pico.ρw / (pico.βs * s - pico.βt) #ρw = ρ* in Reese et al 2018
        b =  g1  ./2 ./ pico.C / pico.ρw ./ (pico.βs * s - pico.βt)
        x = -b .+ sqrt.(b.^2 .- c)
        T1 = pico.T0 .- x
        y = pico.S0 .* x ./ ν ./ λ
        S1 = pico.S0 .- y
        m  = -pico.γT ./ ν ./ λ .* (pico.λ1 .* S1 .+ pico.λ2 .+ pico.λ3 .* zb .- T1)*365.25*24*60^2 #melt rate if all domain in box one

        #set melt rate at box one points
        melt[boxes .== 1] .= m[boxes .== 1]

        #compute the mean temperature and salinity
        T1bar = sum(T1[boxes .== 1])/length(T1[boxes .== 1])
        S1bar = sum(S1[boxes .== 1])/length(S1[boxes .== 1])

    else #use the mean depth, then below are scalar quantities
        zbf = zb
        zb = sum(zbf[boxes .== 1])/length(zbf[boxes .== 1]) #use mean depth in the box
        Tstar = pico.λ1 .* pico.S0 .+ pico.λ2 .+ pico.λ3 .* zb .- pico.T0
        c =  g1 .* Tstar / pico.C / pico.ρw / (pico.βs * s - pico.βt) #ρw = ρ* in Reese et al 2018
        b =  g1  ./2 ./ pico.C / pico.ρw ./ (pico.βs * s - pico.βt)
        x = -b .+ sqrt.(b.^2 .- c)
        T1 = pico.T0 .- x
        y = pico.S0 .* x ./ ν ./ λ
        S1 = pico.S0 .- y

        #set melt rate at box one points
        m  = -pico.γT ./ ν ./ λ .* (pico.λ1 .* S1 .+ pico.λ2 .+ pico.λ3 .* zb .- T1)*365.25*24*60^2 #melt rate if all domain in box one
        melt[boxes .== 1] .= m

        #compute the mean
        T1bar = T1
        S1bar = S1
    end

    #compute the overturning
    q = pico.C * pico.ρw * (pico.βs *(pico.S0 - S1bar) - pico.βt *(pico.T0 - T1bar))

    return T1bar, S1bar, q
end

"""
    function set_box_k_melt_rate!(melt,k, T_prev, S_prev,q, boxes, pico, dx, dy, zb)

Set the melt rate in pico box k > 1.  (Box location specified by boxes array.)
Returns the mean temperature and salinity in this box (used as boundary conditions on the next box)
"""
function set_box_k_melt_rate!(melt,k, T_prev, S_prev,q,boxes,pico, dx, dy, zb)
    m = zeros(size(zb))

    #scalar quantities
    Ak = sum(sum(boxes .== k))*dx *dy
    g1 = Ak * pico.γT
    ν = pico.ρi / pico.ρw
    λ = pico.L / pico.c

    #compute salinity and temp in this box
    if ~(pico.use_box_mean_depth) #don't use mean depth, below are vector quantities
        Tstar = pico.λ1 .* S_prev .+ pico.λ2 .+ pico.λ3 .* zb .- T_prev
        x = -g1 .* Tstar ./ (q + g1 - g1*pico.λ1 * S_prev / ν / λ)
        y = S_prev .* x ./ ν / λ
        Tk = T_prev .- x
        Sk = S_prev .- y

        # compute and set melt rate in this box
        m  = -pico.γT ./ ν ./ λ .* (pico.λ1 .* Sk .+ pico.λ2 .+ pico.λ3 .* zb .- Tk)*365.25*24*60^2 #melt rate if all domain in box one
        melt[boxes .== k] .= m[boxes .== k]

        #update the mean temperature and salinity for the next box
        T_out = sum(Tk[boxes .== k])/length(Tk[boxes .== k])
        S_out = sum(Sk[boxes .== k])/length(Sk[boxes .== k])

    else #use the mean depth, then below are scalar quantities
        zbf = zb
        zb = sum(zbf[boxes .== k])/length(zbf[boxes .== k]) #use mean depth in the box
        Tstar = pico.λ1 .* S_prev .+ pico.λ2 .+ pico.λ3 .* zb .- T_prev
        x = -g1 .* Tstar ./ (q + g1 - g1*pico.λ1 * S_prev / ν / λ)
        y = S_prev .* x ./ ν / λ
        Tk = T_prev .- x
        Sk = S_prev .- y

        # compute and set melt rate in this box
        m  = -pico.γT ./ ν ./ λ .* (pico.λ1 .* Sk .+ pico.λ2 .+ pico.λ3 .* zb .- Tk)*365.25*24*60^2 #melt rate if all domain in box one
        melt[boxes .== k] .= m

        T_out = Tk
        S_out = Sk
    end
    return T_out, S_out
end


text:
while true == true
    print('>')
    readline()
end

max_stars_repo_name: fuzzynoise/ono

max_issues_repo_name: fuzzynoise/ono

max_forks_repo_name: fuzzynoise/ono


text:
const tagfield = collect(fieldnames(tag));

function readtags(atag::tag, xfield::Symbol, xval; collectall::Bool = false)

  if(!collectall && !any(x->x==xfield, tagfield))
    println("$xfield not in $(tagfield[:])");
    return nothing;
    end;

  taglist = tag[];
  tagstack = [atag];
  taglim = [size(atag.attach, 1)];
  tagptr = [0];
  tp_ptr = 1;

  while tp_ptr>0

    local ctag = tagptr[tp_ptr];
    local xtag;
    local zval;
    local subtag =0;

    if(ctag == 0)
      xtag = tagstack[tp_ptr];
      zval = getfield(xtag, xfield);
      (zval == xval || collectall) && push!(taglist, xtag);
      ctag += 1;
      end;

    if(ctag <= taglim[tp_ptr])
      xtag = tagstack[tp_ptr].attach[ctag];
      ctag += 1;
      tagptr[tp_ptr] = ctag;
      zval = getfield(xtag, xfield);
      (zval == xval || collectall) && push!(taglist, xtag);
      subtag = size(xtag.attach,1);

      if(ctag > taglim[tp_ptr])
        pop!(tagstack);#, xtag);
        pop!(taglim); #, subtag);
        pop!(tagptr); #, 1);
        tp_ptr -= 1;
        end;

      if(subtag > 0)
        push!(tagstack, xtag);
        push!(taglim, subtag);
        push!(tagptr, 1);
        tp_ptr += 1;
        end;
      end;
    end;
  return taglist;
  end;


text:
"""
The base abstract type for all optimization problems.
`FS` is a type of a problem's `FitnessScheme`.
"""
abstract type OptimizationProblem{FS<:FitnessScheme} end

# common definitions for `OptimizationProblem`
# (enforce field names of subtypes)
name(p::OptimizationProblem) = p.name
fitness_scheme_type(::Type{<:OptimizationProblem{FS}}) where FS = FS
fitness_type(::Type{P}) where P <: OptimizationProblem = fitness_type(fitness_scheme_type(P))
fitness_scheme(p::OptimizationProblem) = p.fitness_scheme
fitness_type(p::P) where P <: OptimizationProblem = fitness_type(P)
numobjectives(p::OptimizationProblem) = numobjectives(fitness_scheme(p))
search_space(p::OptimizationProblem) = p.ss
numdims(p::OptimizationProblem) = numdims(search_space(p))

"""
Checks if the current best fitness is within the given range
of a known best fitness.
By default the best fitness is unknown.
"""
fitness_is_within_ftol(p::OptimizationProblem, fitness, atol::Number) = false

"""
`OptimizationProblem` with the objective function
defined by some Julia `Function` and search space.

Optionally, a known optimal value could be provided to terminate
the optimization once it is reached.
"""
mutable struct FunctionBasedProblem{FS<:FitnessScheme,SS<:SearchSpace,FO} <: OptimizationProblem{FS}
    objfunc::Function     # Objective function
    name::String
    fitness_scheme::FS
    ss::SS                # search space
    opt_value::FO         # known optimal value or nothing

    function FunctionBasedProblem(objfunc::Function, name::String,
                                  fitness_scheme::FS, ss::SS,
                                  opt_value::FO = nothing) where {FS<:FitnessScheme,SS<:SearchSpace,FO}
        if FO <: Number
            fitness_type(fitness_scheme) == FO ||
                throw(ArgumentError("Fitness type ($(fitness_type(fitness_scheme))) and opt_value type ($(FO)) do not match"))
            if isnafitness(opt_value, fitness_scheme) # if NA fitness is given, the problem is unbouned
                opt_value = nothing
                #throw(ArgumentError("Known optimal value cannot be NA"))
            end
        end
        new{FS,SS,typeof(opt_value)}(objfunc, name, fitness_scheme, ss, opt_value)
    end
end

objfunc(p::FunctionBasedProblem) = p.objfunc

"""
    fitness(x, p::OptimizationProblem)

Evaluates the fitness of a candidate.
"""
fitness(x, p::FunctionBasedProblem) = objfunc(p)(x)

Base.copy(p::FunctionBasedProblem) =
    FunctionBasedProblem(deepcopy(p.objfunc), deepcopy(p.name),
                         p.fitness_scheme, p.ss, p.opt_value)

opt_value(p::FunctionBasedProblem) = p.opt_value

# checks if the optimal fitness is reached for bounded scalar problem
fitness_is_within_ftol(p::FunctionBasedProblem{<:FitnessScheme, <:SearchSpace, <:Number},
                       fitness, atol::Number) = norm(opt_value(p) - fitness) < atol

# no check for unbounded problem (by default, see the definition for OptimizationProblem above)
# fitness_is_within_ftol{FS,SS}(p::FunctionBasedProblem{FS,SS,Void}, fitness::Number, atol::Number) = true


text:
function codegen(::Type{Val{:iterate}}, ::Type{AMSMMALA}, job::BasicMCJob)
  local result::Expr
  burninbody = []
  ambody = []
  smmalapasttensorbody = []
  smmalabody = []
  body = []

  vform = variate_form(job.pstate)
  if vform != Multivariate
    error("Only multivariate parameter states allowed in AMSMMALA code generation")
  end

  push!(body, :(_job.sstate.count += 1))

  if (
    job.tuner.totaltuner.verbose ||
    job.tuner.smmalatuner.verbose ||
    job.tuner.malatuner.verbose ||
    isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
  )
    push!(body, :(_job.sstate.tune.totaltune.proposed += 1))
  end

  push!(
    body,
    :(
      if _job.sstate.count > _job.sampler.t0
        _job.sampler.update!(_job.sstate, _job.pstate, _job.sstate.count, _job.range.nsteps)
      end
    )
  )

  push!(smmalabody, :(_job.sstate.updatetensorcount += 1))

  if job.tuner.smmalatuner.verbose
    push!(smmalabody, :(_job.sstate.tune.smmalatune.proposed += 1))
  end

  push!(smmalapasttensorbody, :(_job.parameter.uptotensorlogtarget!(_job.pstate)))

  if job.sampler.transform != nothing
    push!(smmalapasttensorbody, :(_job.pstate.tensorlogtarget = _job.sampler.transform(_job.pstate.tensorlogtarget)))
  end

  push!(smmalapasttensorbody, :(_job.sstate.oldinvtensor = inv(_job.pstate.tensorlogtarget)))

  push!(smmalapasttensorbody, :(_job.sstate.oldfirstterm = _job.sstate.oldinvtensor*_job.pstate.gradlogtarget))

  push!(smmalapasttensorbody, :(_job.sstate.cholinvtensor = ctranspose(chol(Hermitian(_job.sstate.oldinvtensor)))))

  push!(smmalabody, Expr(:if, :(!_job.sstate.pastupdatetensor), Expr(:block, smmalapasttensorbody...)))

  push!(smmalabody, :(_job.sstate.μ = _job.pstate.value+0.5*_job.sstate.tune.totaltune.step*_job.sstate.oldfirstterm))

  push!(
    smmalabody,
    :(
      _job.sstate.pstate.value =
        _job.sstate.μ+_job.sstate.sqrttunestep*_job.sstate.cholinvtensor*randn(_job.pstate.size)
    )
  )

  push!(smmalabody, :(_job.parameter.uptotensorlogtarget!(_job.sstate.pstate)))

  if job.sampler.transform != nothing
    push!(smmalabody, :(_job.sstate.pstate.tensorlogtarget = _job.sampler.transform(_job.sstate.pstate.tensorlogtarget)))
  end

  push!(smmalabody, :(_job.sstate.newinvtensor = inv(_job.sstate.pstate.tensorlogtarget)))

  push!(smmalabody, :(_job.sstate.ratio = _job.sstate.pstate.logtarget-_job.pstate.logtarget))

  push!(
    smmalabody,
    :(
      _job.sstate.ratio += (
        0.5*(
          logdet(_job.sstate.tune.totaltune.step*_job.sstate.oldinvtensor)
          +dot(
            _job.sstate.pstate.value-_job.sstate.μ,
            _job.pstate.tensorlogtarget*(_job.sstate.pstate.value-_job.sstate.μ)
          )/_job.sstate.tune.totaltune.step
        )
      )
    )
  )

  push!(smmalabody, :(_job.sstate.newfirstterm = _job.sstate.newinvtensor*_job.sstate.pstate.gradlogtarget))

  push!(
    smmalabody,
    :(_job.sstate.μ = _job.sstate.pstate.value+0.5*_job.sstate.tune.totaltune.step*_job.sstate.newfirstterm)
  )

  push!(smmalabody, :(
      _job.sstate.ratio -= (
        0.5*(
          logdet(_job.sstate.tune.totaltune.step*_job.sstate.newinvtensor)
          +dot(
            _job.pstate.value-_job.sstate.μ,
            _job.sstate.pstate.tensorlogtarget*(_job.pstate.value-_job.sstate.μ)
          )/_job.sstate.tune.totaltune.step
        )
      )
    )
  )

  update = []
  noupdate = []

  push!(update, :(_job.pstate.value = copy(_job.sstate.pstate.value)))

  push!(update, :(_job.pstate.gradlogtarget = copy(_job.sstate.pstate.gradlogtarget)))

  if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
    push!(update, :(_job.pstate.gradloglikelihood = copy(_job.sstate.pstate.gradloglikelihood)))
  end

  if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
    push!(update, :(_job.pstate.gradlogprior = copy(_job.sstate.pstate.gradlogprior)))
  end

  push!(update, :(_job.pstate.tensorlogtarget = copy(_job.sstate.pstate.tensorlogtarget)))

  if in(:tensorloglikelihood, job.outopts[:monitor]) && job.parameter.tensorloglikelihood! != nothing
    push!(update, :(_job.pstate.tensorloglikelihood = copy(_job.sstate.pstate.tensorloglikelihood)))
  end

  if in(:tensorlogprior, job.outopts[:monitor]) && job.parameter.tensorlogprior! != nothing
    push!(update, :(_job.pstate.tensorlogprior = copy(_job.sstate.pstate.tensorlogprior)))
  end

  push!(update, :(_job.sstate.oldinvtensor = copy(_job.sstate.newinvtensor)))

  push!(update, :(_job.sstate.cholinvtensor = ctranspose(chol(Hermitian(_job.sstate.newinvtensor)))))

  push!(update, :(_job.sstate.oldfirstterm = copy(_job.sstate.newfirstterm)))

  push!(update, :(_job.pstate.logtarget = _job.sstate.pstate.logtarget))

  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(update, :(_job.pstate.loglikelihood = _job.sstate.pstate.loglikelihood))
  end

  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(update, :(_job.pstate.logprior = _job.sstate.pstate.logprior))
  end

  if in(:accept, job.outopts[:diagnostics])
    push!(update, :(_job.pstate.diagnosticvalues[1] = true))
    push!(noupdate, :(_job.pstate.diagnosticvalues[1] = false))
  end

  if job.tuner.totaltuner.verbose || isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
    push!(update, :(_job.sstate.tune.totaltune.accepted += 1))
  end

  if job.tuner.smmalatuner.verbose
    push!(update, :(_job.sstate.tune.smmalatune.accepted += 1))
  end

  push!(
    smmalabody,
    Expr(
      :if,
      :(_job.sstate.ratio > 0 || (_job.sstate.ratio > log(rand()))),
      Expr(:block, update...),
      Expr(:block, noupdate...)
    )
  )

  if job.tuner.malatuner.verbose
    push!(ambody, :(_job.sstate.tune.malatune.proposed += 1))
  end

  push!(
    ambody,
    :(
      if _job.sstate.pastupdatetensor
        # Once fully migrated to Julia 0.5 or higher, use LinAlg.lowrankupdate instead of chol in order to reduce complexity
        _job.sstate.cholinvtensor = ctranspose(chol(Hermitian(_job.sstate.oldinvtensor)))
      else
        covariance!(
          _job.sstate.newinvtensor,
          _job.sstate.oldinvtensor,
          _job.sstate.count-2,
          _job.pstate.value,
          _job.sstate.lastmean,
          _job.sstate.secondlastmean
        )

        _job.sstate.cholinvtensor = ctranspose(chol(Hermitian(_job.sstate.newinvtensor)))
      end
    )
  )

  push!(
    ambody,
    :(
      _job.sstate.pstate.value =
        _job.pstate.value+_job.sstate.sqrttunestep*_job.sstate.cholinvtensor*randn(_job.pstate.size)
    )
  )

  push!(ambody, :(_job.parameter.logtarget!(_job.sstate.pstate)))

  push!(ambody, :(_job.sstate.ratio = _job.sstate.pstate.logtarget-_job.pstate.logtarget))

  update = []
  noupdate = []

  push!(update, :(_job.pstate.value = copy(_job.sstate.pstate.value)))

  push!(update, :(_job.pstate.logtarget = _job.sstate.pstate.logtarget))

  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(update, :(_job.pstate.loglikelihood = _job.sstate.pstate.loglikelihood))
  end

  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(update, :(_job.pstate.logprior = _job.sstate.pstate.logprior))
  end

  push!(update, :(_job.sstate.oldinvtensor = copy(_job.sstate.newinvtensor)))

  if in(:accept, job.outopts[:diagnostics])
    push!(update, :(_job.pstate.diagnosticvalues[1] = true))
    push!(noupdate, :(_job.pstate.diagnosticvalues[1] = false))
  end

  if job.tuner.totaltuner.verbose || isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
    push!(update, :(_job.sstate.tune.totaltune.accepted += 1))
  end

  if job.tuner.malatuner.verbose
    push!(update, :(_job.sstate.tune.malatune.accepted += 1))
  end

  push!(
    ambody,
    Expr(
      :if,
      :(_job.sstate.ratio > 0 || (_job.sstate.ratio > log(rand()))),
      Expr(:block, update...),
      Expr(:block, noupdate...)
    )
  )

  push!(body, Expr(:if, :(_job.sstate.presentupdatetensor), Expr(:block, smmalabody...), Expr(:block, ambody...)))

  push!(body, :(_job.sstate.secondlastmean = copy(_job.sstate.lastmean)))

  push!(body, :(mean!(_job.sstate.lastmean, _job.sstate.count, _job.pstate.value)))

  push!(body, :(_job.sstate.pastupdatetensor = _job.sstate.presentupdatetensor))

  if job.tuner.totaltuner.verbose || isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
    push!(burninbody, :(rate!(_job.sstate.tune.totaltune)))
  end

  if job.tuner.smmalatuner.verbose
    push!(
      burninbody,
      :(_job.sstate.tune.smmalafrequency = _job.sstate.tune.smmalatune.proposed/_job.sstate.tune.totaltune.proposed)
    )
    push!(burninbody, :(rate!(_job.sstate.tune.smmalatune)))
  end

  if job.tuner.malatuner.verbose
    push!(
      burninbody,
      :(_job.sstate.tune.malafrequency = _job.sstate.tune.malatune.proposed/_job.sstate.tune.totaltune.proposed)
    )
    push!(burninbody, :(rate!(_job.sstate.tune.malatune)))
  end

  if job.tuner.totaltuner.verbose || job.tuner.smmalatuner.verbose || job.tuner.malatuner.verbose
    fmt_tot_iter = format_iteration(ndigits(job.range.burnin))
    fmt_burnin_iter = format_iteration(ndigits(job.tuner.totaltuner.period))
    fmt_perc = format_percentage()

    push!(burninbody, :(println(
      "Burnin iteration ",
      $(fmt_tot_iter)(_job.sstate.count),
      " out of ",
      _job.range.burnin,
      "..."
    )))

    if job.tuner.totaltuner.verbose
      push!(burninbody, :(println(
        "  Total : ",
        $(fmt_burnin_iter)(_job.sstate.tune.totaltune.accepted),
        "/",
        $(fmt_burnin_iter)(_job.sstate.tune.totaltune.proposed),
        " (",
        $(fmt_perc)(100*_job.sstate.tune.totaltune.rate),
        "%) acceptance rate"
      )))
    end

    if job.tuner.smmalatuner.verbose
      push!(burninbody, :(println(
        "  SMMALA: ",
        $(fmt_burnin_iter)(_job.sstate.tune.smmalatune.accepted),
        "/",
        $(fmt_burnin_iter)(_job.sstate.tune.smmalatune.proposed),
        " (",
        $(fmt_perc)(100*_job.sstate.tune.smmalatune.rate),
        "%) acceptance rate, ",
        $(fmt_burnin_iter)(_job.sstate.tune.smmalatune.proposed),
        "/",
        $(fmt_burnin_iter)(_job.sstate.tune.totaltune.proposed),
        " (",
        $(fmt_perc)(100*_job.sstate.tune.smmalafrequency),
        "%) running frequency"
      )))
    end

    if job.tuner.malatuner.verbose
      push!(burninbody, :(println(
        "  MALA  : ",
        $(fmt_burnin_iter)(_job.sstate.tune.malatune.accepted),
        "/",
        $(fmt_burnin_iter)(_job.sstate.tune.malatune.proposed),
        " (",
        $(fmt_perc)(100*_job.sstate.tune.malatune.rate),
        "%) acceptance rate, ",
        $(fmt_burnin_iter)(_job.sstate.tune.malatune.proposed),
        "/",
        $(fmt_burnin_iter)(_job.sstate.tune.totaltune.proposed),
        " (",
        $(fmt_perc)(100*_job.sstate.tune.malafrequency),
        "%) running frequency"
      )))
    end
  end

  if isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
    push!(burninbody, :(tune!(_job.sstate.tune.totaltune, _job.tuner.totaltuner)))
    push!(burninbody, :(_job.sstate.sqrttunestep = sqrt(_job.sstate.tune.totaltune.step)))
  end

  if job.tuner.smmalatuner.verbose
    push!(burninbody, :(reset_burnin!(_job.sstate.tune.smmalatune)))
  end

  if job.tuner.malatuner.verbose
    push!(burninbody, :(reset_burnin!(_job.sstate.tune.malatune)))
  end

  if (
    job.tuner.totaltuner.verbose ||
    job.tuner.smmalatuner.verbose ||
    job.tuner.malatuner.verbose ||
    isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
  )
    push!(burninbody, :(_job.sstate.tune.totaltune.totproposed += _job.sstate.tune.totaltune.proposed))

    push!(burninbody, :(_job.sstate.tune.totaltune.proposed = 0))

    if job.tuner.totaltuner.verbose || isa(job.tuner.totaltuner, AcceptanceRateMCTuner)
      push!(burninbody, :(_job.sstate.tune.totaltune.accepted = 0))
      push!(burninbody, :(_job.sstate.tune.totaltune.rate = NaN))
    end

    push!(
      body,
      Expr(
        :if,
        :(
          _job.sstate.tune.totaltune.totproposed <= _job.range.burnin &&
          mod(_job.sstate.tune.totaltune.proposed, _job.tuner.totaltuner.period) == 0
        ),
        Expr(:block, burninbody...)
      )
    )
  end

  if !job.plain
    push!(body, :(produce()))
  end

  @gensym _iterate

  result = quote
    function $_iterate(_job::BasicMCJob)
      $(body...)
    end
  end

  result
end


text:
using Distributions

function multiple_NB_dgp(RNG, N)
    log_sq_foot = Vector{Float64}(undef, N)  # Use zeros?
    live_in_super = Vector{Float64}(undef, N)  # Use zeros?
    traps = Vector{Int64}(undef, N)  # Use zeros?
    complaints = Vector{Int64}(undef, N)  # Use zeros?

    α = rand(RNG, Normal(log(4), .1))
    β = rand(RNG, Normal(-0.25, .1))
    β_super = rand(RNG, Normal(-0.5, 1))
    inv_ϕ = abs(rand(RNG, Normal(0, 1))) #TODO: fabs == abs?

    for i in 1:N
        log_sq_foot[i] = rand(RNG, Normal(1.5, 0.1))
        live_in_super[i] = rand(RNG, Bernoulli(0.5))
        traps[i] = rand(RNG, Poisson(8))
        r = α + log_sq_foot[i] + β * traps[i] + β_super * live_in_super[i]
        p = inv(inv_ϕ)

        # TODO: Log parametrization
        # complaints[i] = rand(RNG, NegativeBinomial(r, p))
    end

    return (α = α,
            β = β,
            β_super = β_super,
            log_sq_foot = log_sq_foot,
            live_in_super = live_in_super,
            traps = traps,
            complaints = complaints)
end


text:
@testset "DataArray" begin
    v = [1, 2, 3, 4]
    dv = DataArray(v, falses(size(v)))

    m = [1 2; 3 4]
    dm = DataArray(m, falses(size(m)))

    t = Array{Int}(uninitialized, 2, 2, 2)
    t[1:2, 1:2, 1:2] = 1
    dt = DataArray(t, falses(size(t)))

    dv = DataArray(v)
    dv = DataArray(v, [false, false, false, false])

    dv = DataArray(Int, 2)
    dm = DataArray(Int, 2, 2)
    dt = DataArray(Int, 2, 2, 2)

    similar(dv)
    similar(dm)
    similar(dt)

    similar(dv, 2)
    similar(dm, 2, 2)
    similar(dt, 2, 2, 2)

    @test isequal(DataArray([missing, missing], [true, true]), DataArray([missing, missing], [false, false]))
    @test isequal(DataArray(Any[1, missing], [false, true]), DataArray(Any[1, missing], [false, false]))
    @test isequal(DataArray(Any[1, 2], [false, true]), DataArray(Any[1, missing], [false, false]))

    x = DataArray([9, 9, 8])
    y = DataArray([1, 9, 3, 2, 2])
    @test append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = DataArray([9, 9, 8])
    y = [1, 9, 3, 2, 2]
    @test append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = @data [1, 2, missing]
    y = @data [3, missing, 5]
    @test isequal(copy(x), x)
    @test isequal(copy!(y, x), x)

    x = @data [1, missing, -2, 1, missing, 4]
    @test isequal(unique(x), @data [1, missing, -2, 4])
    @test isequal(unique(reverse(x)), @data [4, missing, 1, -2])
    @test isequal(unique(skipmissing(x)), @data [1, -2, 4])
    @test isequal(unique(reverse(collect(skipmissing(x)))), @data [4, 1, -2])
    @test isequal(levels(x), @data [1, -2, 4])
    @test isequal(levels(reverse(x)), @data [4, 1, -2])

    # check case where only missing occurs in final position
    @test isequal(unique(@data [1, 2, 1, missing]), @data [1, 2, missing])

    # Test copy!
    function nonbits(dv)
        ret = similar(dv, Integer)
        for i = 1:length(dv)
            # if !ismissing(dv, i)
                ret[i] = dv[i]
            # end
        end
        ret
    end
    set1 = Any[@data([1, missing, 3]),
               @data([missing, 5]), @data([1, 2, 3, 4, 5]), data(Int[]),
               @data([missing, 5, 3]), @data([1, 5, 3])]
    set2 = map(nonbits, set1)
    set3 = map(pdata, set1)

    for (dest, src, bigsrc, emptysrc, res1, res2) in Any[set1, set2, set3]
        @test isequal(copy!(copy(dest), src), res1)
        @test isequal(copy!(copy(dest), 1, src), res1)

        @test isequal(copy!(copy(dest), 2, src, 2), res2)
        @test isequal(copy!(copy(dest), 2, src, 2, 1), res2)

        @test isequal(copy!(copy(dest), 99, src, 99, 0), dest)

        @test isequal(copy!(copy(dest), 1, emptysrc), dest)
        @test_throws BoundsError copy!(dest, 1, emptysrc, 1)

        for idx in [0, 4]
            @test_throws BoundsError copy!(dest, idx, src)
            @test_throws BoundsError copy!(dest, idx, src, 1)
            @test_throws BoundsError copy!(dest, idx, src, 1, 1)
            @test_throws BoundsError copy!(dest, 1, src, idx)
            @test_throws BoundsError copy!(dest, 1, src, idx, 1)
        end

        @test_throws ArgumentError copy!(dest, 1, src, 1, -1)

        @test_throws BoundsError copy!(dest, bigsrc)

        @test_throws BoundsError copy!(dest, 3, src)
        @test_throws BoundsError copy!(dest, 3, src, 1)
        @test_throws BoundsError copy!(dest, 3, src, 1, 2)
        @test_throws BoundsError copy!(dest, 1, src, 2, 2)
    end

    # Inferrability of map (#276)
    @test eltype(map(x -> x > 1, @data [1, 2])) == Union{Bool,Missing}

    @testset "Issue #278" begin
        x = @data ones(4)
        @test parent(view(x, :)).data === x.data
    end
end


text:
using NewtonsMethod
using Test
using LinearAlgebra, Statistics, Compat, ForwardDiff

@testset "NewtonsMethod.jl" begin
    #Several @test for the root of a known functions
    #1
    f1(x) = (x-1.0)^3
    f1′(x) = 3.0*(x-1.0)^2
    x1₀ = 0.1
    @test newtonroot(f1, f1′, x₀ = x1₀).root ≈ 0.9999998168636032
    @test newtonroot(f1,      x₀ = x1₀).root ≈ 0.9999998168636032
    @test newtonroot(f1, f1′, x₀ = BigFloat(x1₀), tolerance = 1E-80).root ≈ BigFloat(1.0)
    @test newtonroot(f1,      x₀ = BigFloat(x1₀), tolerance = 1E-80).root ≈ BigFloat(1.0)
    #2
    f2(x) = log(x)
    f2′(x) = 1.0/x
    x2₀ = 2.0
    @test newtonroot(f2, f2′, x₀ = x2₀).root ≈ 1.0
    @test newtonroot(f2,      x₀ = x2₀).root ≈ 1.0
    @test newtonroot(f2, f2′, x₀ = BigFloat(x2₀), tolerance = 1E-50).root ≈ BigFloat(1.0)
    @test newtonroot(f2,      x₀ = BigFloat(x2₀), tolerance = 1E-50).root ≈ BigFloat(1.0)
    #3
    f3(x) = (sin(x)/x)
    f3′(x) = (x*cos(x) - sin(x))/(x^2)
    x3₀ = 1.0
    @test newtonroot(f3, f3′, x₀ = x3₀).root ≈ 3.141592653589793
    @test newtonroot(f3,      x₀ = x3₀).root ≈ 3.141592653589793
    @test newtonroot(f3, f3′, x₀ = BigFloat(x3₀), tolerance = 1E-50).root ≈ pi
    @test newtonroot(f3,      x₀ = BigFloat(x3₀), tolerance = 1E-50).root ≈ pi
    #4
    f4(x) = (x+10.0)*(x^2+2.0)
    f4′(x) = 3x^2 + 20x + 2.0
    x4₀ = 10.0
    @test newtonroot(f4, f4′, x₀ = x4₀).root ≈ -10.0
    @test newtonroot(f4,      x₀ = x4₀).root ≈ -10.0
    @test newtonroot(f4, f4′, x₀ = BigFloat(x4₀), tolerance = 1E-100).root ≈ BigFloat(-10.0)
    @test newtonroot(f4,      x₀ = BigFloat(x4₀), tolerance = 1E-100).root ≈ BigFloat(-10.0)

    #Some @test of non-convergence
    #5
    f5(x) = x^2 + 2.0
    f5′(x) = 2x
    x5₀ = 10.0
    @test newtonroot(f5, f5′, x₀ = x5₀).root == nothing
    @test newtonroot(f5,      x₀ = x5₀).root == nothing
    @test newtonroot(f5, f5′, x₀ = BigFloat(x5₀), tolerance = 1E-100).root == nothing
    @test newtonroot(f5,      x₀ = BigFloat(x5₀), tolerance = 1E-100).root == nothing
    #6
    f6(x) = 2.0
    f6′(x) = 0.0
    x6₀ = 10.0
    @test newtonroot(f6, f6′, x₀ = x6₀).root == nothing
    @test newtonroot(f6,      x₀ = x6₀).root == nothing
    @test newtonroot(f6, f6′, x₀ = BigFloat(x6₀), tolerance = 1E-100).root == nothing
    @test newtonroot(f6,      x₀ = BigFloat(x6₀), tolerance = 1E-100).root == nothing

    #Some @test for maxiter
    #1
    @test newtonroot(f1, f1′, x₀ = x1₀, maxiter = 10).root == nothing
    @test newtonroot(f1,      x₀ = x1₀, maxiter = 10).root == nothing
    @test newtonroot(f1, f1′, x₀ = BigFloat(x1₀), tolerance = 1E-80, maxiter = 100).root == nothing
    @test newtonroot(f1,      x₀ = BigFloat(x1₀), tolerance = 1E-80, maxiter = 100).root == nothing
    #4
    @test newtonroot(f4, f4′, x₀ = x4₀, maxiter = 10).root == nothing
    @test newtonroot(f4,      x₀ = x4₀, maxiter = 10).root == nothing
    @test newtonroot(f4, f4′, x₀ = BigFloat(x4₀), tolerance = 1E-100, maxiter = 100).root ≈ BigFloat(-10.0)
    @test newtonroot(f4,      x₀ = BigFloat(x4₀), tolerance = 1E-100, maxiter = 100).root ≈ BigFloat(-10.0)

    #Some @test for tolerance and BigFloat
    #2
    @test !(newtonroot(f2, f2′, x₀ = x2₀, tolerance = 1E-2).root ≈ 1.0)
    @test !(newtonroot(f2,      x₀ = x2₀, tolerance = 1E-2).root ≈ 1.0)
    #3
    @test !(BigFloat(0.0) + newtonroot(f3, f3′, x₀ = x3₀, tolerance = 1E-50).root ≈ pi)
    @test !(BigFloat(0.0) + newtonroot(f3,      x₀ = x3₀, tolerance = 1E-50).root ≈ pi)
    @test !(newtonroot(f3, f3′, x₀ = BigFloat(x3₀), tolerance = 1E-10).root ≈ pi)
    @test !(newtonroot(f3,      x₀ = BigFloat(x3₀), tolerance = 1E-10).root ≈ pi)
end


text:
using LinearAlgebra    # 使用线性代数库
f(n)=rand(n)⋅rand(n)   # 定义函数 f
for i in 1:3           # 执行 3 次
  @time @show f(100)   # 显示结果与时间
end

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes