text:
using FFMPEG

tosecond(x::Hour) = Float64(Dates.value(Second(x)))
tosecond(x::Minute) = Float64(Dates.value(Second(x)))
tosecond(x::T) where {T<:TimePeriod} = x/convert(T, Second(1))

# extract(_, __::Missing, ___) = nothing

ffmpeg(start, sourcefile, targetfile) = ffmpeg_exe(` -y -nostats -loglevel 8 -ss $start -i $sourcefile -vf 'yadif=1,format=gray,scale=sar*iw:ih' -pix_fmt gray -vframes 1 $targetfile`)

extract(targetpathfile, t::Temporal{WholeVideo, I}, coffeesource) where {I <: Instantaneous} = ffmpeg(tosecond(start(t.time)), joinpath(coffeesource, t.video.file.name), targetpathfile)

function extract(targetpathfile, t::Temporal{<:FragmentedVideo, I}, coffeesource) where {I <: Instantaneous}
    x = start(t.time)
    for vf in files(t.video)
        d = duration(vf)
        if x ≤ d
            ffmpeg(tosecond(x), joinpath(coffeesource, vf.name), targetpathfile)
            break
        end
        x -= d
    end
end

ffmpeg(start, duration, sourcefile, targetpathfile) = ffmpeg_exe(` -y -nostats -loglevel 8 -ss $start -i $sourcefile -t $duration -r 2 -vf 'yadif=1,format=gray,scale=sar*iw:ih' -pix_fmt gray $targetpathfile`)

getsourcefile(coffeesource, t) = joinpath(coffeesource, t)
function getsourcetarget(coffeesource, t, targetpath)
    sourcefile = getsourcefile(coffeesource, t)
    sourcefile, joinpath(targetpath, "img%03d.png")
end
extract(targetpath, t::Temporal{WholeVideo, P}, coffeesource) where {P <: Prolonged} = ffmpeg(tosecond(start(t.time)), tosecond(duration(t.time)), getsourcetarget(coffeesource, t.video.file.name, targetpath)...)

function extract(targetpath, t::Temporal{<:FragmentedVideo, P}, coffeesource) where {P <: Prolonged}
    t1 = start(t.time)
    t2 = stop(t.time)
    files = Iterators.Stateful(t.video.files)
    for vf in files
        d = duration(vf)
        if t1 ≤ d
            if t2 ≤ d
                ffmpeg(tosecond(t1), tosecond(t2 - t1), getsourcetarget(coffeesource, vf.name, targetpath)...)
                break
            else
                ffmpeg(tosecond(t1), tosecond(d), getsourcetarget(coffeesource, vf.name, targetpath)...)
                t2 -= d
                for _vf in files
                    _d = duration(_vf)
                    if t2 ≤ _d
                        ffmpeg(tosecond(Millisecond(0)), tosecond(t2), getsourcetarget(coffeesource, _vf.name, targetpath)...)
                        break
                    end
                    ffmpeg(tosecond(Millisecond(0)), tosecond(_d), getsourcetarget(coffeesource, _vf.name, targetpath)...)
                    t2 -= _d
                end
            end
        end
        t1 -= d
        t2 -= d
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

#=
using JuliaFEM.API
using JuliaFEM.Interfaces
=#

#= TODO: Fix test
@testset "test basic workflow" begin

    # basic workflow, copied from test_solver.jl
    model = Model("Piston Calculation")

    add_node!(model, 1, [0.0, 0.0])
    model.nodes[2] = [1.0, 0.0]
    add_node!(model, 3, [1.0, 1.0])
    model.nodes[4] = [0.0, 1.0]

    # create elements
    e1 = Element(1, [1, 2, 3, 4], :Quad4)
    e2 = Element(2, [1, 2], :Seg2)
    e3 = Element(3, [3, 4], :Seg2)
    model.elements[1] = e1
    model.elements[2] = e2
    add_element!(model, 3, :Seg2, [3, 4])

    # element set
    elset = ElementSet("body", [e1, e2])
    elset4 = ElementSet("set_material", [e1])

    model.elsets["body"] = elset
    add_element_set!(model, "heat_flux", [2])
    add_element_set!(model, "constant_temp", [e3])
    add_element_set!(model, elset4)

    # material properties
    material = Material("myMaterial")
    material["temperature thermal conductivity"] = 6.0
    material["density"] = 36.0
    add_material!(model, "set_material", material)

    # Create problem
    field_problem = Simulation(:HeatProblem)
    add_element_set!(field_problem, "body") # is this necessary?

    # boundary conditions
    bc = DirichletBC("constant_temp", "temperature" => 0.0)
    ne = NeumannBC("heat_flux", "temperature flux" => ((0.0 =>   0.0),
                                                  (1.0 => 600.0)))

    # LoadCase
    add_boundary_condition!(field_problem, bc)
    add_boundary_condition!(field_problem, ne)

    # Get solver
    add_solver!(field_problem, :LinearSolver)

    # add case
    add_simulation!(model, "Heat problem", field_problem)
    #model.load_cases["Heat problem"] = field_problem

    # Solve problem
    solve!(model, "Heat problem", 1.0)
    xi = [0.0, -1.0]
    T = model.elements[1].results("temperature", xi, 1.0)
    X = model.elements[2].results("geometry", xi, 1.0)
    info("Temperature at point X = $X is T = $T")
    #@test isapprox(T, 200.0)
end
=#

#= TODO: Fix test
@testset "test reading piston model using API" begin
    abaqus_input = open(parse_abaqus, "./geometry/piston/piston_8789_P1.inp")
    model = Model("Piston Calculation", abaqus_input)
    @test length(keys(model.elsets)) == 4
    @test length(keys(model.nsets)) == 1
    @test length(keys(model.elements)) == 37331
    @test length(keys(model.nodes)) == 8789
end
=#

#= TODO: Fix test
function test_piston_107168()
    abaqus_input = open(parse_abaqus, "./geometry/piston/piston_107168_P2.inp")

    model = Model("Piston Calculation", abaqus_input)
    @test length(keys(model.elsets)) == 3
    @test length(keys(model.nsets)) == 1
    @test length(keys(model.elements)) == 65948
    @test length(keys(model.nodes)) == 107168
end
=#

function slow_test_something_that_takes_long_time()
    info("This test is SLOW.")
    test_piston_170168()
end

ext: jl


text:
#
# Types supporting structural definition of models and their components
#

# Objects with a `name` attribute
@class NamedObj <: MimiClass begin
    name::Symbol
end

"""
    nameof(obj::NamedDef) = obj.name

Return the name of `def`.  `NamedDef`s include `DatumDef`, `ComponentDef`,
`CompositeComponentDef`, and `VariableDefReference` and `ParameterDefReference`.
"""
Base.nameof(obj::AbstractNamedObj) = obj.name

# Deprecate old definition in favor of standard name
@deprecate name(obj::AbstractNamedObj) nameof(obj)

# Similar structure is used for variables and parameters (parameters merely adds `default`)
@class mutable DatumDef <: NamedObj begin
    comp_path::Union{Nothing, ComponentPath}
    datatype::DataType
    dim_names::Vector{Symbol}
    description::String
    unit::String
end

@class mutable VariableDef <: DatumDef

@class mutable ParameterDef <: DatumDef begin
    # ParameterDef adds a default value, which can be specified in @defcomp
    default::Any
end

@class mutable ComponentDef <: NamedObj begin
    comp_id::Union{Nothing, ComponentId}    # allow anonynous top-level (composite) ComponentDefs (must be referenced by a ModelDef)
    comp_path::Union{Nothing, ComponentPath}
    dim_dict::OrderedDict{Symbol, Union{Nothing, Dimension}}
    namespace::OrderedDict{Symbol, Any}
    first::Union{Nothing, Int}
    last::Union{Nothing, Int}
    is_uniform::Bool

    # Store a reference to the AbstractCompositeComponent that contains this comp def.
    # That type is defined later, so we declare Any here. Parent is `nothing` for
    # detached (i.e., "template") components and is set when added to a composite.
    parent::Any

    function ComponentDef(self::ComponentDef, comp_id::Nothing)
        error("Leaf ComponentDef objects must have a valid ComponentId name (not nothing)")
    end

    # ComponentDefs are created "empty". Elements are subsequently added.
    function ComponentDef(self::AbstractComponentDef, comp_id::Union{Nothing, ComponentId}=nothing;
                          name::Union{Nothing, Symbol}=nothing)
        if comp_id === nothing
            # ModelDefs are anonymous, but since they're gensym'd, they can claim the Mimi package
            comp_id = ComponentId(Mimi, @or(name, gensym(nameof(typeof(self)))))
        end

        name = @or(name, comp_id.comp_name)
        NamedObj(self, name)

        self.comp_id = comp_id
        self.comp_path = nothing    # this is set in add_comp!() and ModelDef()

        self.dim_dict   = OrderedDict{Symbol, Union{Nothing, Dimension}}()
        self.namespace = OrderedDict{Symbol, Any}()
        self.first = self.last = nothing
        self.is_uniform = true
        self.parent = nothing
        return self
    end

    function ComponentDef(comp_id::Union{Nothing, ComponentId};
                          name::Union{Nothing, Symbol}=nothing)
        self = new()
        return ComponentDef(self, comp_id; name=name)
    end
end

ns(obj::AbstractComponentDef) = obj.namespace
comp_id(obj::AbstractComponentDef) = obj.comp_id
pathof(obj::AbstractComponentDef) = obj.comp_path
dim_dict(obj::AbstractComponentDef) = obj.dim_dict
first_period(obj::AbstractComponentDef) = obj.first
last_period(obj::AbstractComponentDef) = obj.last
isuniform(obj::AbstractComponentDef) = obj.is_uniform

Base.parent(obj::AbstractComponentDef) = obj.parent

# Used by @defcomposite to communicate subcomponent information
struct SubComponent <: MimiStruct
    module_obj::Union{Nothing, Module}
    comp_name::Symbol
    alias::Union{Nothing, Symbol}
    bindings::Vector{Pair{Symbol, Any}}
end

# Stores references to the name of a component variable or parameter
# and the ComponentPath of the component in which it is defined
@class DatumReference <: NamedObj begin
    # name::Symbol is inherited from NamedObj
    root::AbstractComponentDef
    comp_path::ComponentPath
end

@class ParameterDefReference <: DatumReference

@class VariableDefReference  <: DatumReference

function _dereference(ref::AbstractDatumReference)
    comp = find_comp(ref)
    return comp[ref.name]
end

# Might not be useful
# convert(::Type{VariableDef},  ref::VariableDefReference)  = _dereference(ref)
# convert(::Type{ParameterDef}, ref::ParameterDefReference) = _dereference(ref)

# Define type aliases to avoid repeating these in several places
global const Binding = Pair{AbstractDatumReference, Union{Int, Float64, AbstractDatumReference}}

# Define which types can appear in the namespace dict for leaf and composite compdefs
global const LeafNamespaceElement      = AbstractDatumDef
global const CompositeNamespaceElement = Union{AbstractComponentDef, AbstractDatumReference}
global const NamespaceElement          = Union{LeafNamespaceElement, CompositeNamespaceElement}

@class mutable CompositeComponentDef <: ComponentDef begin
    bindings::Vector{Binding}

    internal_param_conns::Vector{InternalParameterConnection}
    external_param_conns::Vector{ExternalParameterConnection}
    external_params::Dict{Symbol, ModelParameter}               # TBD: make key (ComponentPath, Symbol)?

    # Names of external params that the ConnectorComps will use as their :input2 parameters.
    backups::Vector{Symbol}

    sorted_comps::Union{Nothing, Vector{Symbol}}

    function CompositeComponentDef(comp_id::Union{Nothing, ComponentId}=nothing)
        self = new()
        CompositeComponentDef(self, comp_id)
        return self
    end

    function CompositeComponentDef(self::AbstractCompositeComponentDef, comp_id::Union{Nothing, ComponentId}=nothing)
        ComponentDef(self, comp_id) # call superclass' initializer

        self.comp_path = ComponentPath(self.name)
        self.bindings = Vector{Binding}()
        self.internal_param_conns = Vector{InternalParameterConnection}()
        self.external_param_conns = Vector{ExternalParameterConnection}()
        self.external_params = Dict{Symbol, ModelParameter}()
        self.backups = Vector{Symbol}()
        self.sorted_comps = nothing
    end
end

# Used by @defcomposite
function CompositeComponentDef(comp_id::ComponentId, alias::Symbol, subcomps::Vector{SubComponent},
                               calling_module::Module)
    # @info "CompositeComponentDef($comp_id, $alias, $subcomps)"
    composite = CompositeComponentDef(comp_id)

    for c in subcomps
        # @info "subcomp $c: module: $(printable(c.module_obj)), calling module: $(nameof(calling_module))"
        comp_module = @or(c.module_obj, calling_module)
        subcomp_id = ComponentId(comp_module, c.comp_name)
        subcomp = compdef(subcomp_id)
        add_comp!(composite, subcomp, @or(c.alias, c.comp_name))
    end
    return composite
end

add_backup!(obj::AbstractCompositeComponentDef, backup) = push!(obj.backups, backup)

internal_param_conns(obj::AbstractCompositeComponentDef) = obj.internal_param_conns
external_param_conns(obj::AbstractCompositeComponentDef) = obj.external_param_conns

external_params(obj::AbstractCompositeComponentDef) = obj.external_params

is_leaf(c::AbstractComponentDef) = true
is_leaf(c::AbstractCompositeComponentDef) = false
is_composite(c::AbstractComponentDef) = !is_leaf(c)

ComponentPath(obj::AbstractCompositeComponentDef, name::Symbol) = ComponentPath(obj.comp_path, name)

ComponentPath(obj::AbstractCompositeComponentDef, path::AbstractString) = comp_path(obj, path)

ComponentPath(obj::AbstractCompositeComponentDef, names::Symbol...) = ComponentPath(obj.comp_path.names..., names...)

@class mutable ModelDef <: CompositeComponentDef begin
    number_type::DataType
    dirty::Bool

    function ModelDef(number_type::DataType=Float64)
        self = new()
        CompositeComponentDef(self)  # call super's initializer
        return ModelDef(self, number_type, false)       # call @class-generated method
    end
end

#
# Reference types offer a more convenient syntax for interrogating Components.
#

# A container for a component, for interacting with it within a model.
@class ComponentReference <: MimiClass begin
    parent::AbstractComponentDef
    comp_path::ComponentPath
end

function ComponentReference(parent::AbstractComponentDef, name::Symbol)
    return ComponentReference(parent, ComponentPath(parent.comp_path, name))
end

# A container for a variable within a component, to improve connect_param! aesthetics,
# by supporting subscripting notation via getindex & setindex .
@class VariableReference <: ComponentReference begin
    var_name::Symbol
end


text:
using Test

using LinearAlgebra
using ProximalAlgorithms
using ProximalAlgorithms: NesterovExtrapolation, initialize, SimpleNesterovSequence, FixedNesterovSequence, AdaptiveNesterovSequence

for sequence_type in [SimpleNesterovSequence, FixedNesterovSequence]
    for R in [Float32, Float64]
        @testset "Nesterov accel. ($sequence_type, $R)" begin
            H = R[
                0.63287    0.330934   -0.156908   -0.294776    0.10761;
                0.330934   0.673201    0.0459778   0.231011   -0.235265;
               -0.156908   0.0459778   0.635812   -0.232261   -0.388775;
               -0.294776   0.231011   -0.232261    0.726854   -0.0691783;
                0.10761   -0.235265   -0.388775   -0.0691783   0.336262;
            ]
            l = R[1.0, 2.0, 3.0, 4.0, 5.0]

            n = length(l)

            f(x) = dot(x, H * x) / 2 + dot(x, l)
            grad_f(x) = H * x + l

            x_star = -H \ l
            f_star = f(x_star)

            Lip = opnorm(H)
            gamma = 1 / Lip
            x = zeros(R, n)
            y = x

            norm_initial_error = norm(x_star - x)

            @inferred initialize(NesterovExtrapolation(sequence_type), x)

            seq = sequence_type{R}()
            @test eltype(seq) == R

            @inferred Iterators.peel(seq)

            for (it, coeff) in Iterators.take(enumerate(seq), 100)
                if it == 1
                    @test iszero(coeff)
                end

                grad_f_y = grad_f(y)
                x_prev = x
                x = y - gamma * grad_f_y
                y = x + coeff * (x - x_prev)

                # test that iterates satisfy Thm 4.4 from Beck, Teboulle (2009)
                @test f(x) - f_star <= 2 / (gamma * (it + 1)^2) * norm_initial_error^2
            end
        end
    end
end

@testset "Nesterov adaptive accel." begin
    R = Float64
    fixed_gamma = R(1.7)

    adaptive_seq = AdaptiveNesterovSequence(R(0))
    fixed_seq = FixedNesterovSequence{R}()
    for el in Iterators.take(fixed_seq, 20)
        @test isapprox(el, ProximalAlgorithms.next!(adaptive_seq, fixed_gamma))
    end

    m = R(1)
    adaptive_seq = AdaptiveNesterovSequence(m)
    for _ in 1:20
        @test isapprox(
            (1 - sqrt(m * fixed_gamma)) / (1 + sqrt(m * fixed_gamma)),
            ProximalAlgorithms.next!(adaptive_seq, fixed_gamma)
        )
    end
end


text:
# This example demonstrated how to calculate the bubble diagram of free electrons using the Monte Carlo module

using QuantumStatistics, LinearAlgebra, Random, Printf, BenchmarkTools, InteractiveUtils, Parameters, StaticArrays
# using ProfileView
using PyCall

special = pyimport("scipy.special")
const l = 7
coeff = SVector{l+1,Float64}(special.legendre(l).c)
exp = zeros(l+1)
for i=1:l+1
    exp[i]=l+1-i
end
exp = SVector{l+1,Float64}(exp)

function legendre(x)
    return dot((x .^ exp) ,coeff)
end

function test()
    x = 1.3
    y1 = legendre(x)
    y2 = special.legendre(l)(x)
    @printf("%10.6f , %10.6f\n", y1, y2)
end

test()
const Steps = 1e8

const e0 = sqrt(2)  # electric charge
const me = 0.5  # electron mass
const dim = 3    # dimension (D=2 or 3, doesn't work for other D!!!)
const spin = 2  # number of spins

const rs = 2.0
const kF =1.0# (dim == 3) ? (9π / (2spin))^(1 / 3) / rs : sqrt(4 / spin) / rs
const EF = kF^2 / (2me)
const β = 1000.0 / kF^2
const mass2 = 0.01

const e2 = rs * 2.0/(9π/4.0)^(1.0/3)

# function P_(z,u)
#     return 1.0+(1.0-z^2+u^2)/4/z*log(((1+z)^2+u^2)/((1-z)^2+u^2)) - u *atan(2*u/(u^2+z^2-1))
# end

function P_(z::Float64)
    return 1.0+(1.0-z^2)/2.0/z*log(((1+z))/(abs(1-z)))
end

# function Π(q,ω)
#     return me * kF/2/π^2 * P_(q/2.0/kF, me*ω/q/kF)
# end

function Π0(q::Float64)
    return me * kF/2/π^2 * P_(q/2.0/kF)
end

function integrand(config)
    if config.curr != 1
        error("impossible")
    end

    T = config.var[1]
    z = 2kF * T[1]

    k, p = kF, kF

    W = 1.0/( z^2/(4π*e2) + Π0(z))

    return 2kF * z/(k*p) * legendre((k^2+p^2-z^2)/2/k/p) * W
end

function measure(config)
    obs = config.observable
    factor = 1.0 / config.reweight[config.curr]
    weight = integrand(config)
    obs[1] += weight / abs(weight) * factor
end

function run(steps)

    T = MonteCarlo.Tau(1.0, 1.0 / 2.0)

    dof = [[1,],] # degrees of freedom of the normalization diagram and the bubble
    obs = zeros(Float64, 1) # observable for the normalization diagram and the bubble

    config = MonteCarlo.Configuration(steps, (T,), dof, obs)
    avg, std = MonteCarlo.sample(config, integrand, measure; print=0, Nblock=16)
    # @profview MonteCarlo.sample(config, integrand, measure; print=0, Nblock=1)
    # sleep(100)

    if isnothing(avg) == false
        @printf("%10.6f ± %10.6f\n", avg[1], std[1])
    end
end

run(Steps)
#@time run(Steps)


text:
module UtilTest

using Parametron
using Test

using Parametron: sort_and_combine!
using Base.Sort: QuickSort

combinepair(x, y) = first(x) => last(x) + last(y)

@testset "sort_and_combine!" begin
    n = 100
    @test sort_and_combine!([2 => 1.0, 2 => 2.0]; combine=combinepair, by=first) == [2 => 3.0]
    @test sort_and_combine!([3 => 1.0, 2 => 4.0, 3 => 2.0, 1 => 2.0], combine=combinepair, by=first) == [1 => 2.0, 2 => 4.0, 3 => 3.0]
    @test length(sort_and_combine!(Pair.(1 : n, rand(n)), combine=combinepair, by=first)) == n
    @test length(sort_and_combine!(Pair.(fill(3, n), rand(n)), combine=combinepair, by=first)) == 1

    v = [rand(1 : round(Int, n / 2)) => rand() for i = 1 : n]
    global allocs
    for i = 1 : 2
        vcopy = copy(v)
        allocs = @allocated sort_and_combine!(vcopy; combine=combinepair, by=first, alg=QuickSort)
    end
    @test allocs == 0
end

end # module


text:
module MFUtilities

    export comp_inv_pr_pick, get_n_elem, BatchIterator

    function get_n_elem(block, nusers, nitems)
        n_elem_row = zeros(Int64,nusers)
        n_elem_col = zeros(Int64,nitems)
        n = size(block)[1]
        for j in 1:n
            u,i,r = block[j,:]
            u = round(Int64,u)
            i = round(Int64,i)
            n_elem_row[u] += 1
            n_elem_col[i] += 1
        end
        return n_elem_row, n_elem_col
    end

    # compute local inverse pick prob
    function comp_inv_pr_pick(n_elem_row,n_elem_col,batchsize)
        N = sum(n_elem_row)
        try
            @assert N == sum(n_elem_col)
        catch
            println("$(myid()) assertion failed")
        end

        Pr_us = 1 - (1 - n_elem_row./ N ) .^ batchsize
        Pr_is = 1 - (1 - n_elem_col./ N ) .^ batchsize

        return 1 ./ Pr_us, 1 ./ Pr_is
    end

    """
    BatchIterator type to get minibatches
    """
    type BatchIterator
      nitems::Int
      batchsize::Int
      niters::Float64 # float to allow infinite iterations.
      randomize::Bool

      function BatchIterator(nitems::Int,batchsize::Int,niters::Float64;randomize::Bool=false)
        @assert batchsize > 0
        @assert nitems >= batchsize
        @assert niters >= 0
        new(nitems,batchsize,niters,randomize)
      end
    end

    function Base.start(B::BatchIterator)
      (0,0,B.randomize ? randperm(B.nitems) : (1:B.nitems))
    end
    function Base.next(B::BatchIterator,state)
      iter = state[1]
      id = state[2]
      order = state[3]

      # +1 is necessary since the Base.start uses a 0 index
      # The remainder means that the batches wrap around if the batch size doesn't evenly divide the data size
      batch = order[rem(Int[i for i = id:id+B.batchsize-1],B.nitems)+1]

      # get new permutation for next epoch
      if id+B.batchsize >= B.nitems
        if B.randomize
          order = randperm(B.nitems)
        end
      end

      id = rem(id+B.batchsize,B.nitems)
      iter += 1

      # Return the indices of the next minibatch, and update the state (which actually points to the next minibatch)
      (batch,(iter,id,order))
    end

    # Iterating is done when the iteration in the state matches the preset iteration limit in the iterator object
    function Base.done(B::BatchIterator,state)
      state[1] >= B.niters
    end

end


text:
#UTILITY PROGRESSMETER, REPORTS WORKER NUMBER AND CURRENT ITERATE
mutable struct ProgressHMM{T<:Real} <: ProgressMeter.AbstractProgress
    thresh::T
    dt::AbstractFloat
    val::T
    counter::Integer
    triggered::Bool
    tfirst::AbstractFloat
    tlast::AbstractFloat
    printed::Bool        # true if we have issued at least one status update
    desc::AbstractString # prefix to the percentage, e.g.  "Computing..."
    color::Symbol        # default to green
    output::IO           # output stream into which the progress is written
    numprintedvalues::Integer   # num values printed below progress in last iteration
    offset::Integer             # position offset of progress bar (default is 0)
    steptime::AbstractFloat
    function ProgressHMM{T}(thresh;
                               dt::Real=0.1,
                               desc::AbstractString="Progress: ",
                               color::Symbol=:green,
                               output::IO=stderr,
                               offset::Integer=0,
                               start_it::Integer=1) where T
        tfirst = tlast = time()
        printed = false
        new{T}(thresh, dt, typemax(T), start_it, false, tfirst, tlast, printed, desc, color, output, 0, offset, 0.0)
    end
end

ProgressHMM(thresh::Real, dt::Real=0.1, desc::AbstractString="Progress: ",
         color::Symbol=:green, output::IO=stderr;
         offset::Integer=0, start_it::Integer=1) = ProgressHMM{typeof(thresh)}(thresh, dt=dt, desc=desc, color=color, output=output, offset=offset, start_it=start_it)

ProgressHMM(thresh::Real, desc::AbstractString, offset::Integer=0, start_it::Integer=1) = ProgressHMM{typeof(thresh)}(thresh, desc=desc, offset=offset, start_it=start_it)

function update!(p::ProgressHMM, val, steptime; options...)
    p.val = val
    p.counter += 1
    p.steptime = steptime
    updateProgress!(p; options...)
end

function updateProgress!(p::ProgressHMM; showvalues = Any[], valuecolor = :blue, offset::Integer = p.offset, keep = (offset == 0))
    p.offset = offset
    t = time()
    if p.val <= p.thresh && !p.triggered
        p.triggered = true
        if p.printed
            p.triggered = true
            dur = ProgressMeter.durationstring(t-p.tfirst)
            msg = @sprintf "%s Time: %s (%d iterations)" p.desc dur p.counter
            print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
            ProgressMeter.move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
            ProgressMeter.printover(p.output, msg, p.color)
            ProgressMeter.printvalues!(p, showvalues; color = valuecolor)
            if keep
                println(p.output)
            else
                print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))
            end
        end
        return
    end

    if t > p.tlast+p.dt && !p.triggered
        elapsed_time = t - p.tfirst
        msg = @sprintf "%s (convergence: %g => %g, iterate: %g, step time: %s)" p.desc p.val p.thresh p.counter hmss(p.steptime)
        print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
        ProgressMeter.move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
        ProgressMeter.printover(p.output, msg, p.color)
        ProgressMeter.printvalues!(p, showvalues; color = valuecolor)
        print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))
        # Compensate for any overhead of printing. This can be
        # especially important if you're running over a slow network
        # connection.
        p.tlast = t + 2*(time()-t)
        p.printed = true
    end
end

            function hmss(dt)
                isnan(dt) && return "NaN"
                (h,r) = divrem(dt,60*60)
                (m,r) = divrem(r, 60)
                (isnan(h)||isnan(m)||isnan(r)) && return "NaN"
                string(Int(h),":",Int(m),":",Int(ceil(r)))
            end


text:
mutable struct GenerateContext{T <: Trace, K <: ConstrainedSelection, P <: Parameters} <: ExecutionContext
    tr::T
    select::K
    weight::Float64
    score::Float64
    visited::Visitor
    params::P
    GenerateContext(tr::T, select::K) where {T <: Trace, K <: ConstrainedSelection} = new{T, K, EmptyParameters}(tr, select, 0.0, 0.0, Visitor(), Parameters())
    GenerateContext(tr::T, select::K, params::P) where {T <: Trace, K <: ConstrainedSelection, P <: Parameters} = new{T, K, P}(tr, select, 0.0, 0.0, Visitor(), params)
end
Generate(select::ConstrainedSelection) = GenerateContext(Trace(), select, Parameters())
Generate(select::ConstrainedSelection, params) = GenerateContext(Trace(), select, params)
Generate(tr::Trace, select::ConstrainedSelection) = GenerateContext(tr, select)

# ------------ Convenience ------------ #

function generate(sel::L, fn::Function, args...) where L <: ConstrainedSelection
    ctx = Generate(sel)
    ret = ctx(fn, args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.weight
end

function generate(sel::L, params, fn::Function, args...) where L <: ConstrainedSelection
    ctx = Generate(sel, params)
    ret = ctx(fn, args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.weight
end

function generate(sel::L, fn::typeof(rand), d::Distribution{K}) where {L <: ConstrainedSelection, K}
    ctx = Generate(sel)
    addr = gensym()
    ret = ctx(fn, addr, d)
    return ret, get_top(ctx.tr, addr), ctx.weight
end

function generate(sel::L, fn::typeof(markov), call::Function, len::Int, args...) where L <: ConstrainedSelection
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = Generate(v_sel)
    ret = ctx(fn, addr, call, len, args...)
    return ret, get_sub(ctx.tr, addr), ctx.weight
end

function generate(sel::L, params, fn::typeof(markov), call::Function, len::Int, args...) where L <: ConstrainedSelection
    addr = gensym()
    v_sel = selection(addr => sel)
    v_ps = parameters(addr => params)
    ctx = Generate(v_sel, v_ps)
    ret = ctx(fn, addr, call, len, args...)
    return ret, get_sub(ctx.tr, addr), ctx.weight
end

function generate(sel::L, fn::typeof(plate), call::Function, args::Vector) where L <: ConstrainedSelection
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = Generate(v_sel)
    ret = ctx(fn, addr, call, args)
    return ret, get_sub(ctx.tr, addr), ctx.weight
end

function generate(sel::L, params, fn::typeof(plate), call::Function, args::Vector) where L <: ConstrainedSelection
    addr = gensym()
    v_sel = selection(addr => sel)
    v_ps = parameters(addr => params)
    ctx = Generate(v_sel, v_ps)
    ret = ctx(fn, addr, call, args)
    return ret, get_sub(ctx.tr, addr), ctx.weight
end

function generate(sel::L, fn::typeof(plate), d::Distribution{K}, len::Int) where {L <: ConstrainedSelection, K}
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = Generate(v_sel)
    ret = ctx(fn, addr, d, len)
    return ret, get_sub(ctx.tr, addr), ctx.weight
end

# ------------ includes ------------ #

include("hierarchical/generate.jl")
include("plate/generate.jl")
include("markov/generate.jl")
include("conditional/generate.jl")
include("factor/generate.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct GenerateContext{T <: Trace, K <: ConstrainedSelection, P <: Parameters} <: ExecutionContext
     tr::T
     select::K
     weight::Float64
     score::Float64
     visited::Visitor
     params::P
end
```
`GenerateContext` is used to generate traces, as well as record and accumulate likelihood weights given observations at addressed randomness.

Inner constructors:
```julia
GenerateContext(tr::T, select::K) where {T <: Trace, K <: ConstrainedSelection} = new{T, K}(tr, select, 0.0, Visitor(), Parameters())
GenerateContext(tr::T, select::K, params::P) where {T <: Trace, K <: ConstrainedSelection, P <: Parameters} = new{T, K, P}(tr, select, 0.0, Visitor(), params)
```
Outer constructors:
```julia
Generate(select::ConstrainedSelection) = GenerateContext(Trace(), select)
Generate(select::ConstrainedSelection, params) = GenerateContext(Trace(), select, params)
Generate(tr::Trace, select::ConstrainedSelection) = GenerateContext(tr, select)
```
""", GenerateContext)

@doc(
"""
```julia
ret, cl, w = generate(sel::L, fn::Function, args...; params = Parameters()) where L <: ConstrainedSelection
ret, cs, w = generate(sel::L, fn::typeof(rand), d::Distribution{K}; params = Parameters()) where {L <: ConstrainedSelection, K}
ret, v_cl, w = generate(sel::L, fn::typeof(markov), call::Function, len::Int, args...; params = Parameters()) where L <: ConstrainedSelection
ret, v_cl, w = generate(sel::L, fn::typeof(plate), call::Function, args::Vector; params = Parameters()) where L <: ConstrainedSelection
ret, v_cl, w = generate(sel::L, fn::typeof(plate), d::Distribution{K}, len::Int; params = Parameters()) where {L <: ConstrainedSelection, K}
```

`generate` provides an API to the `GenerateContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, a `RecordSite` instance specialized to the call, and the score/weight `w` computed with respect to the constraints `sel`.
""", generate)


text:
#
# The functions in this file initialize the iterations
#
"""
Newtoninit: set up Newton's method
"""
function Newtoninit(
    x0,
    dx,
    F!,
    J!,
    solver,
    sham,
    armmax,
    armfix,
    resdec,
    maxit,
    printerr,
    pdata,
    jfact,
    keepsolhist,
)
    #
    #   Initialize the iteration.
    #
    n = length(x0)
    x = copy(x0)
    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    ItRules = (
        dx = dx,
        f = F!,
        fp = J!,
        solver = solver,
        sham = sham,
        armmax = armmax,
        armfix = armfix,
        resdec = resdec,
        maxit = maxit,
        printerr = printerr,
        pdata = pdata,
        fact = jfact,
    )
    return (ItRules, x, n, solhist)
end

"""
Secantinit(x0, dx, f, solver,
         armmax, armfix, maxit, printerr, pdata, jfact)
"""
function Secantinit(x0, dx, f, solver, armmax, armfix, maxit, printerr, pdata, jfact)
    n = length(x0)
    x = copy(x0)
    ItRules = (
        f = f,
        solver = solver,
        armmax = armmax,
        armfix = armfix,
        maxit = maxit,
        printerr = printerr,
        pdata = pdata,
        fact = jfact,
    )
    return (ItRules, x, n)
end

"""
PTCinit(x0, dx, F!, J!, delta0, maxit, pdata, jfact, keepsolhist)

PTCinit: get organized for PTC
"""
function PTCinit(x0, dx, F!, J!, delta0, maxit, pdata, jfact, keepsolhist, jknowsdt = false)
    #
    #   Initialize the iteration.
    #
    n = length(x0)
    x = copy(x0)
    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    ItRules = (
        dx = dx,
        f = F!,
        fp = J!,
        delta0 = delta0,
        maxit = maxit,
        pdata = pdata,
        fact = jfact,
        jknowsdt = jknowsdt,
    )
    return (ItRules, x, n, solhist)
end

"""
Newton_Krylov_Init( x0, dx, F!, Jvec, Pvec, pside, lsolver, eta,
    fixedeta, armmax, armfix, maxit, lmaxit, printerr, pdata, keepsolhist)

Newton_Krylov_Init: set up nsoli
"""
function Newton_Krylov_Init(
    x0,
    dx,
    F!,
    Jvec,
    Pvec,
    pside,
    lsolver,
    eta,
    fixedeta,
    armmax,
    armfix,
    maxit,
    lmaxit,
    printerr,
    pdata,
    keepsolhist,
)
    #
    #   Initialize the iteration.
    #
    eta > 0 || error("eta must be positive")
    n = length(x0)
    x = copy(x0)
    tmp1 = zeros(n)
    tmp2 = zeros(n)
    tmp3 = zeros(n)
    tmp4 = zeros(n)
    kl_store = (tmp1, tmp2, tmp3, tmp4)
    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    ((lmaxit == -1) && (lsolver=="bicgstab") ) && (lmaxit = 5)
    ItRules = (
        dx = dx,
        f = F!,
        Jvec = Jvec,
        Pvec = Pvec,
        pside = pside,
        lsolver = lsolver,
        kl_store = kl_store,
        eta = eta,
        fixedeta = fixedeta,
        lmaxit = lmaxit,
        armmax = armmax,
        armfix = armfix,
        maxit = maxit,
        printerr = printerr,
        pdata = pdata,
    )
    return (ItRules, x, n, solhist)
end

"""
solhistinit(n, maxit, x)

Am I keeping the solution history? If so, allocate the space.
"""
function solhistinit(n, maxit, x)
    #
    # If you are keeping a solution history, make some room for it.
    #
    solhist = zeros(n, maxit + 1)
    @views solhist[:, 1] .= x
    return solhist
end


text:
module Requires

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
    @eval Base.Experimental.@compiler_options compile=min optimize=0 infer=false
end

using UUIDs

function _include_path(relpath::String)
    # Reproduces include()'s runtime relative path logic
    # See Base._include_dependency()
    prev = Base.source_path(nothing)
    if prev === nothing
        path = abspath(relpath)
    else
        path = normpath(joinpath(dirname(prev), relpath))
    end
end

"""
    @include("somefile.jl")

Behaves like `include`, but caches the target file content at macro expansion
time, and uses this as a fallback when the file doesn't exist at runtime. This
is useful when compiling a sysimg. The argument `"somefile.jl"` must be a
string literal, not an expression.

`@require` blocks insert this automatically when you use `include`.
"""
macro include(relpath::String)
    compiletime_path = joinpath(dirname(String(__source__.file)), relpath)
    s = String(read(compiletime_path))
    quote
        # NB: Runtime include path may differ from the compile-time macro
        # expansion path if the source has been relocated.
        runtime_path = _include_path($relpath)
        if isfile(runtime_path)
            # NB: For Revise compatibility, include($relpath) needs to be
            # emitted where $relpath is a string *literal*.
            $(esc(:(include($relpath))))
        else
            include_string($__module__, $s, $relpath)
        end
    end
end

include("init.jl")
include("require.jl")

function __init__()
    push!(package_callbacks, loadpkg)
end

if isprecompiling()
    precompile(loadpkg, (Base.PkgId,)) || @warn "Requires failed to precompile `loadpkg`"
    precompile(withpath, (Any, String)) || @warn "Requires failed to precompile `withpath`"
    precompile(err, (Any, Module, String)) || @warn "Requires failed to precompile `err`"
    precompile(parsepkg, (Expr,)) || @warn "Requires failed to precompile `parsepkg`"
    precompile(listenpkg, (Any, Base.PkgId)) || @warn "Requires failed to precompile `listenpkg`"
    precompile(callbacks, (Base.PkgId,)) || @warn "Requires failed to precompile `callbacks`"
    precompile(withnotifications, (String, Module, String, String, Expr)) || @warn "Requires failed to precompile `withnotifications`"
    precompile(replace_include, (Expr, LineNumberNode)) || @warn "Requires failed to precompile `replace_include`"
    precompile(getfield(Requires, Symbol("@require")), (LineNumberNode, Module, Expr, Any)) || @warn "Requires failed to precompile `@require`"

    precompile(_include_path, (String,)) || @warn "Requires failed to precompile `_include_path`"
    precompile(getfield(Requires, Symbol("@include")), (LineNumberNode, Module, String)) || @warn "Requires failed to precompile `@include`"

    precompile(__init__, ()) || @warn "Requires failed to precompile `__init__`"
end

end # module


text:
using JuliaUtils
using Documenter

makedocs(;
    modules=[JuliaUtils],
    authors="Frank Urbach",
    repo="https://GitHub.com/FrankUrbach/JuliaUtils.jl/blob/{commit}{path}#L{line}",
    sitename="JuliaUtils.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://FrankUrbach.github.io/JuliaUtils.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="GitHub.com/FrankUrbach/JuliaUtils.jl",
    devbranch="master",
)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

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

    function LinearKernel(; kwargs...)
        obj = new(meta(LinearKernel), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct LinearKernel
const __meta_LinearKernel = Ref{ProtoMeta}()
function meta(::Type{LinearKernel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_LinearKernel)
            __meta_LinearKernel[] = target = ProtoMeta(LinearKernel)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, LinearKernel, 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_LinearKernel[]
    end
end

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

    function RBFKernel(; kwargs...)
        obj = new(meta(RBFKernel), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct RBFKernel
const __meta_RBFKernel = Ref{ProtoMeta}()
function meta(::Type{RBFKernel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RBFKernel)
            __meta_RBFKernel[] = target = ProtoMeta(RBFKernel)
            allflds = Pair{Symbol,Union{Type,String}}[:gamma => Float64]
            meta(target, RBFKernel, 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_RBFKernel[]
    end
end
function Base.getproperty(obj::RBFKernel, name::Symbol)
    if name === :gamma
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

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

    function PolyKernel(; kwargs...)
        obj = new(meta(PolyKernel), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PolyKernel
const __meta_PolyKernel = Ref{ProtoMeta}()
function meta(::Type{PolyKernel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PolyKernel)
            __meta_PolyKernel[] = target = ProtoMeta(PolyKernel)
            allflds = Pair{Symbol,Union{Type,String}}[:degree => Int32, :c => Float64, :gamma => Float64]
            meta(target, PolyKernel, 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_PolyKernel[]
    end
end
function Base.getproperty(obj::PolyKernel, name::Symbol)
    if name === :degree
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :c
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :gamma
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

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

    function SigmoidKernel(; kwargs...)
        obj = new(meta(SigmoidKernel), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SigmoidKernel
const __meta_SigmoidKernel = Ref{ProtoMeta}()
function meta(::Type{SigmoidKernel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SigmoidKernel)
            __meta_SigmoidKernel[] = target = ProtoMeta(SigmoidKernel)
            allflds = Pair{Symbol,Union{Type,String}}[:gamma => Float64, :c => Float64]
            meta(target, SigmoidKernel, 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_SigmoidKernel[]
    end
end
function Base.getproperty(obj::SigmoidKernel, name::Symbol)
    if name === :gamma
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :c
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

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

    function Kernel(; kwargs...)
        obj = new(meta(Kernel), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Kernel
const __meta_Kernel = Ref{ProtoMeta}()
function meta(::Type{Kernel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Kernel)
            __meta_Kernel[] = target = ProtoMeta(Kernel)
            allflds = Pair{Symbol,Union{Type,String}}[:linearKernel => LinearKernel, :rbfKernel => RBFKernel, :polyKernel => PolyKernel, :sigmoidKernel => SigmoidKernel]
            oneofs = Int[1,1,1,1]
            oneof_names = Symbol[Symbol("kernel")]
            meta(target, Kernel, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Kernel[]
    end
end
function Base.getproperty(obj::Kernel, name::Symbol)
    if name === :linearKernel
        return (obj.__protobuf_jl_internal_values[name])::LinearKernel
    elseif name === :rbfKernel
        return (obj.__protobuf_jl_internal_values[name])::RBFKernel
    elseif name === :polyKernel
        return (obj.__protobuf_jl_internal_values[name])::PolyKernel
    elseif name === :sigmoidKernel
        return (obj.__protobuf_jl_internal_values[name])::SigmoidKernel
    else
        getfield(obj, name)
    end
end

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

    function SparseNode(; kwargs...)
        obj = new(meta(SparseNode), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SparseNode
const __meta_SparseNode = Ref{ProtoMeta}()
function meta(::Type{SparseNode})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SparseNode)
            __meta_SparseNode[] = target = ProtoMeta(SparseNode)
            allflds = Pair{Symbol,Union{Type,String}}[:index => Int32, :value => Float64]
            meta(target, SparseNode, 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_SparseNode[]
    end
end
function Base.getproperty(obj::SparseNode, name::Symbol)
    if name === :index
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

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

    function SparseVector(; kwargs...)
        obj = new(meta(SparseVector), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SparseVector
const __meta_SparseVector = Ref{ProtoMeta}()
function meta(::Type{SparseVector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SparseVector)
            __meta_SparseVector[] = target = ProtoMeta(SparseVector)
            allflds = Pair{Symbol,Union{Type,String}}[:nodes => Base.Vector{SparseNode}]
            meta(target, SparseVector, 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_SparseVector[]
    end
end
function Base.getproperty(obj::SparseVector, name::Symbol)
    if name === :nodes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{SparseNode}
    else
        getfield(obj, name)
    end
end

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

    function SparseSupportVectors(; kwargs...)
        obj = new(meta(SparseSupportVectors), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SparseSupportVectors
const __meta_SparseSupportVectors = Ref{ProtoMeta}()
function meta(::Type{SparseSupportVectors})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SparseSupportVectors)
            __meta_SparseSupportVectors[] = target = ProtoMeta(SparseSupportVectors)
            allflds = Pair{Symbol,Union{Type,String}}[:vectors => Base.Vector{SparseVector}]
            meta(target, SparseSupportVectors, 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_SparseSupportVectors[]
    end
end
function Base.getproperty(obj::SparseSupportVectors, name::Symbol)
    if name === :vectors
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{SparseVector}
    else
        getfield(obj, name)
    end
end

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

    function DenseVector(; kwargs...)
        obj = new(meta(DenseVector), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct DenseVector
const __meta_DenseVector = Ref{ProtoMeta}()
function meta(::Type{DenseVector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DenseVector)
            __meta_DenseVector[] = target = ProtoMeta(DenseVector)
            pack = Symbol[:values]
            allflds = Pair{Symbol,Union{Type,String}}[:values => Base.Vector{Float64}]
            meta(target, DenseVector, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_DenseVector[]
    end
end
function Base.getproperty(obj::DenseVector, name::Symbol)
    if name === :values
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

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

    function DenseSupportVectors(; kwargs...)
        obj = new(meta(DenseSupportVectors), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct DenseSupportVectors
const __meta_DenseSupportVectors = Ref{ProtoMeta}()
function meta(::Type{DenseSupportVectors})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DenseSupportVectors)
            __meta_DenseSupportVectors[] = target = ProtoMeta(DenseSupportVectors)
            allflds = Pair{Symbol,Union{Type,String}}[:vectors => Base.Vector{DenseVector}]
            meta(target, DenseSupportVectors, 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_DenseSupportVectors[]
    end
end
function Base.getproperty(obj::DenseSupportVectors, name::Symbol)
    if name === :vectors
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{DenseVector}
    else
        getfield(obj, name)
    end
end

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

    function Coefficients(; kwargs...)
        obj = new(meta(Coefficients), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Coefficients
const __meta_Coefficients = Ref{ProtoMeta}()
function meta(::Type{Coefficients})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Coefficients)
            __meta_Coefficients[] = target = ProtoMeta(Coefficients)
            pack = Symbol[:alpha]
            allflds = Pair{Symbol,Union{Type,String}}[:alpha => Base.Vector{Float64}]
            meta(target, Coefficients, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Coefficients[]
    end
end
function Base.getproperty(obj::Coefficients, name::Symbol)
    if name === :alpha
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

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

    function SupportVectorRegressor(; kwargs...)
        obj = new(meta(SupportVectorRegressor), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SupportVectorRegressor
const __meta_SupportVectorRegressor = Ref{ProtoMeta}()
function meta(::Type{SupportVectorRegressor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SupportVectorRegressor)
            __meta_SupportVectorRegressor[] = target = ProtoMeta(SupportVectorRegressor)
            allflds = Pair{Symbol,Union{Type,String}}[:kernel => Kernel, :sparseSupportVectors => SparseSupportVectors, :denseSupportVectors => DenseSupportVectors, :coefficients => Coefficients, :rho => Float64]
            oneofs = Int[0,1,1,0,0]
            oneof_names = Symbol[Symbol("supportVectors")]
            meta(target, SupportVectorRegressor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_SupportVectorRegressor[]
    end
end
function Base.getproperty(obj::SupportVectorRegressor, name::Symbol)
    if name === :kernel
        return (obj.__protobuf_jl_internal_values[name])::Kernel
    elseif name === :sparseSupportVectors
        return (obj.__protobuf_jl_internal_values[name])::SparseSupportVectors
    elseif name === :denseSupportVectors
        return (obj.__protobuf_jl_internal_values[name])::DenseSupportVectors
    elseif name === :coefficients
        return (obj.__protobuf_jl_internal_values[name])::Coefficients
    elseif name === :rho
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

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

    function SupportVectorClassifier(; kwargs...)
        obj = new(meta(SupportVectorClassifier), 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))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SupportVectorClassifier
const __meta_SupportVectorClassifier = Ref{ProtoMeta}()
function meta(::Type{SupportVectorClassifier})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SupportVectorClassifier)
            __meta_SupportVectorClassifier[] = target = ProtoMeta(SupportVectorClassifier)
            fnum = Int[1,2,3,4,5,6,7,8,100,101]
            pack = Symbol[:numberOfSupportVectorsPerClass,:rho,:probA,:probB]
            allflds = Pair{Symbol,Union{Type,String}}[:kernel => Kernel, :numberOfSupportVectorsPerClass => Base.Vector{Int32}, :sparseSupportVectors => SparseSupportVectors, :denseSupportVectors => DenseSupportVectors, :coefficients => Base.Vector{Coefficients}, :rho => Base.Vector{Float64}, :probA => Base.Vector{Float64}, :probB => Base.Vector{Float64}, :stringClassLabels => StringVector, :int64ClassLabels => Int64Vector]
            oneofs = Int[0,0,1,1,0,0,0,0,2,2]
            oneof_names = Symbol[Symbol("supportVectors"),Symbol("ClassLabels")]
            meta(target, SupportVectorClassifier, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_SupportVectorClassifier[]
    end
end
function Base.getproperty(obj::SupportVectorClassifier, name::Symbol)
    if name === :kernel
        return (obj.__protobuf_jl_internal_values[name])::Kernel
    elseif name === :numberOfSupportVectorsPerClass
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    elseif name === :sparseSupportVectors
        return (obj.__protobuf_jl_internal_values[name])::SparseSupportVectors
    elseif name === :denseSupportVectors
        return (obj.__protobuf_jl_internal_values[name])::DenseSupportVectors
    elseif name === :coefficients
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Coefficients}
    elseif name === :rho
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :probA
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :probB
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :stringClassLabels
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :int64ClassLabels
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    else
        getfield(obj, name)
    end
end

export LinearKernel, RBFKernel, PolyKernel, SigmoidKernel, Kernel, SparseNode, SparseVector, SparseSupportVectors, DenseVector, DenseSupportVectors, Coefficients, SupportVectorRegressor, SupportVectorClassifier


text:
function price2string(price, precision; round_up=false)
    if price <= 0.0001
        if round_up
            return "0."*lpad(Int(ceil(price * 10^precision)), precision, "0")
        else
            return "0."*lpad(Int(floor(price * 10^precision)), precision, "0")
        end
    else
        if round_up
            return string(ceil(price, digits=precision))
        else
            return string(floor(price, digits=precision))
        end
    end
end
export price2string

function get_precision(client, symbol)
    info = client.get_symbol_info(symbol)

    min_qty = parse(Float64, info["filters"][4]["minNotional"])
    stepSize  = parse(Float64, info["filters"][3]["stepSize"])
    buy_precision = 0
    while true
        if stepSize * 10^buy_precision == 1
            break
        end
        buy_precision += 1
    end

    tickSize = parse(Float64, info["filters"][1]["tickSize"])
    sell_precision = 0
    while tickSize * 10^sell_precision < 1
        sell_precision += 1
    end
    return buy_precision, sell_precision, min_qty, stepSize, tickSize
end
export get_precision


text:
const featurenames = [
    :DN_HistogramMode_5 # These shouldn't change between versions of catch22
    :DN_HistogramMode_10
    :CO_Embed2_Dist_tau_d_expfit_meandiff
    :CO_f1ecac
    :CO_FirstMin_ac
    :CO_HistogramAMI_even_2_5
    :CO_trev_1_num
    :DN_OutlierInclude_p_001_mdrmd
    :DN_OutlierInclude_n_001_mdrmd
    :FC_LocalSimple_mean1_tauresrat
    :FC_LocalSimple_mean3_stderr
    :IN_AutoMutualInfoStats_40_gaussian_fmmi
    :MD_hrv_classic_pnn40
    :SB_BinaryStats_diff_longstretch0
    :SB_BinaryStats_mean_longstretch1
    :SB_MotifThree_quantile_hh
    :SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1
    :SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1
    :SP_Summaries_welch_rect_area_5_1
    :SP_Summaries_welch_rect_centroid
    :SB_TransitionMatrix_3ac_sumdiagcov
    :PD_PeriodicityWang_th0_01]

const featuretypes = Dict(featurenames .=> [
                                    #DN_HistogramMode_5
                                    Cdouble
                                    #DN_HistogramMode_10
                                    Cdouble
                                    #CO_Embed2_Dist_tau_d_expfit_meandiff
                                    Cdouble
                                    #CO_f1ecac
                                    Cint
                                    #CO_FirstMin_ac
                                    Cint
                                    #CO_HistogramAMI_even_2_5
                                    Cdouble
                                    #CO_trev_1_num
                                    Cdouble
                                    #DN_OutlierInclude_p_001_mdrmd
                                    Cdouble
                                    #DN_OutlierInclude_n_001_mdrmd
                                    Cdouble
                                    #FC_LocalSimple_mean1_tauresrat
                                    Cdouble
                                    #FC_LocalSimple_mean3_stderr
                                    Cdouble
                                    #IN_AutoMutualInfoStats_40_gaussian_fmmi
                                    Cdouble
                                    #MD_hrv_classic_pnn40
                                    Cdouble
                                    #SB_BinaryStats_diff_longstretch0
                                    Cdouble
                                    #SB_BinaryStats_mean_longstretch1
                                    Cdouble
                                    #SB_MotifThree_quantile_hh
                                    Cdouble
                                    #SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1
                                    Cdouble
                                    #SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1
                                    Cdouble
                                    #SP_Summaries_welch_rect_area_5_1
                                    Cdouble
                                    #SP_Summaries_welch_rect_centroid
                                    Cdouble
                                    #SB_TransitionMatrix_3ac_sumdiagcov
                                    Cdouble
                                    #PD_PeriodicityWang_th0_01
                                    Cint])

"""
    Catch22.featurekeywords
A vector listing keywords of features as vectors of strings.
"""
const featurekeywords = [ # See hctsa
                                    #DN_HistogramMode_5
                                    ["distribution", "location"],
                                    #DN_HistogramMode_10
                                    ["distribution", "location"],
                                    #CO_Embed2_Dist_tau_d_expfit_meandiff
                                    ["correlation", "embedding"],
                                    #CO_f1ecac
                                    ["correlation", "timescale"],
                                    #CO_FirstMin_ac
                                    ["correlation", "timescale"],
                                    #CO_HistogramAMI_even_2_5
                                    ["information", "correlation", "AMI"],
                                    #CO_trev_1_num
                                    ["correlation", "nonlinear"],
                                    #DN_OutlierInclude_p_001_mdrmd
                                    ["distribution", "outliers"],
                                    #DN_OutlierInclude_n_001_mdrmd
                                    ["distribution", "outliers"],
                                    #FC_LocalSimple_mean1_tauresrat
                                    ["forecasting"],
                                    #FC_LocalSimple_mean3_stderr
                                    ["forecasting"],
                                    #IN_AutoMutualInfoStats_40_gaussian_fmmi
                                    ["information", "correlation", "AMI"],
                                    #MD_hrv_classic_pnn40
                                    ["medical"],
                                    #SB_BinaryStats_diff_longstretch0
                                    ["distribution", "stationarity"],
                                    #SB_BinaryStats_mean_longstretch1
                                    ["distribution", "stationarity"],
                                    #SB_MotifThree_quantile_hh
                                    ["symbolic", "motifs"],
                                    #SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1
                                    ["scaling"],
                                    #SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1
                                    ["scaling"],
                                    #SP_Summaries_welch_rect_area_5_1
                                    ["FourierSpectrum"],
                                    #SP_Summaries_welch_rect_centroid
                                    ["FourierSpectrum"],
                                    #SB_TransitionMatrix_3ac_sumdiagcov
                                    ["symbolic", "transitionmat"],
                                    #PD_PeriodicityWang_th0_01
                                    ["periodicity", "spline"]]

"""
    Catch22.featuredescriptions
A vector listing short descriptions of each feature, as strings.
"""
const featuredescriptions = [  # See catch22 paper
                                    #DN_HistogramMode_5
                                    "Mode of z-scored distribution (5-bin histogram)"
                                    #DN_HistogramMode_10
                                    "Mode of z-scored distribution (10-bin histogram)"
                                    #CO_Embed2_Dist_tau_d_expfit_meandiff
                                    "Exponential fit to successive distances in 2-d embedding space"
                                    #CO_f1ecac
                                    "First 1/𝑒 crossing of autocorrelation function"
                                    #CO_FirstMin_ac
                                    "First minimum of autocorrelation function"
                                    #CO_HistogramAMI_even_2_5
                                    "Automutual information, 𝑚=2,𝜏=5"
                                    #CO_trev_1_num
                                    "Time-reversibility statistic, ⟨(𝑥ₜ₊₁−𝑥ₜ)³⟩ₜ"
                                    #DN_OutlierInclude_p_001_mdrmd
                                    "Time intervals between successive extreme events above the mean"
                                    #DN_OutlierInclude_n_001_mdrmd
                                    "Time intervals between successive extreme events below the mean"
                                    #FC_LocalSimple_mean1_tauresrat
                                    "Change in correlation length after iterative differencing"
                                    #FC_LocalSimple_mean3_stderr
                                    "Mean error from a rolling 3-sample mean forecasting"
                                    #IN_AutoMutualInfoStats_40_gaussian_fmmi
                                    "First minimum of the automutual information function"
                                    #MD_hrv_classic_pnn40
                                    "Proportion of successive differences exceeding 0.04𝜎 (Mietus 2002)"
                                    #SB_BinaryStats_diff_longstretch0
                                    "Longest period of successive incremental decreases"
                                    #SB_BinaryStats_mean_longstretch1
                                    "Longest period of consecutive values above the mean"
                                    #SB_MotifThree_quantile_hh
                                    "Shannon entropy of two successive letters in equiprobable 3-letter symbolization"
                                    #SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1
                                    "Proportion of slower timescale fluctuations that scale with linearly rescaled range fits"
                                    #SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1
                                    "Proportion of slower timescale fluctuations that scale with DFA (50% sampling)"
                                    #SP_Summaries_welch_rect_area_5_1
                                    "Total power in lowest fifth of frequencies in the Fourier power spectrum"
                                    #SP_Summaries_welch_rect_centroid
                                    "Centroid of the Fourier power spectrum"
                                    #SB_TransitionMatrix_3ac_sumdiagcov
                                    "Trace of covariance of transition matrix between symbols in 3-letter alphabet"
                                    #PD_PeriodicityWang_th0_01
                                    "Periodicity measure of (Wang et al. 2007)"]


text:
export endEncoding!

abstract type MtlCommandEncoder end

Base.:(==)(a::T, b::T) where {T <: MtlCommandEncoder} = a.handle == b.handle
Base.hash(q::MtlCommandEncoder, h::UInt) = hash(q.handle, h)

function unsafe_destroy!(cce::MtlCommandEncoder)
    if cce.handle !== C_NULL
        mtRelease(cce)
    end
end

endEncoding!(ce::MtlCommandEncoder) = mtCommandEncoderEndEncoding(ce);
Base.close(ce::MtlCommandEncoder) = endEncoding!(ce)


text:
module ComputerVision

using StaticArrays, GeometryTypes, LinearAlgebra, Optim

include("util.jl")
include("model/image.jl")
include("model/camera.jl")
include("model/plane.jl")
include("model/world.jl")

export AbstractAnalgoueImage,
       AbstractDigitalImage,
       AbstractImage,
       AnalogueImage

#export  get_data

export  AbstractCamera,
        AbstractCameraModel,
        AbstractIntrinsicParameters,
        AbstractExtrinsicParameters,
        IntrinsicParameters,
        ExtrinsicParameters,
        Camera

export Plane,
       AbstractWorld,
       PrimitiveWorld

export  get_xaxis,
        set_xaxis!,
        get_yaxis,
        set_yaxis!,
        get_zaxis,
        set_zaxis!,
        get_points,
        set_points!,
        get_planes,
        set_planes!,
        get_normal,
        get_distance

export  hom⁻¹,
        hom

export  get_focal_length,
        set_focal_length!,
        get_width,
        set_width!,
        get_height!,
        set_height!,
        get_centroid,
        set_centroid,
        get_intrinsics,
        set_intrinsics,
        get_extrinsics,
        set_extrinsics,
        get_e₁,
        set_e₁!,
        get_e₂,
        set_e₂!,
        get_e₃,
        set_e₃!,
        get_origin,
        set_origin,
        get_model,
        set_model!,
        get_image_type,
        set_image_type!

end # module


text:
using Test
using Cambrian

cfg = get_config("test.yaml")

function test_individual(ind::Individual)
    ind2 = mutate(ind, 0.0)
    @test all(ind2.fitness[1] .== -Inf)
    @test all(ind.genes .== ind2.genes)

    ind2 = mutate(ind, 1.0)
    @test all(ind2.fitness[1] == -Inf)
    @test all(ind.genes .!= ind2.genes)

    ind3 = crossover(ind, ind2)
    @test all(ind3.fitness[1] == -Inf)
    @test all((ind3.genes .== ind.genes) .| (ind3.genes .== ind2.genes))
end

@testset "bool individual" begin
    b_ind = BoolIndividual(cfg)
    fval = rand()
    b_ind.fitness[1] = fval

    b_str = string(b_ind)
    @test typeof(b_str) == String

    b_ind2 = BoolIndividual(b_str)
    @test typeof(b_ind2) == BoolIndividual
    @test b_ind2.fitness[1] == fval
    @test all(b_ind.genes .== b_ind2.genes)

    test_individual(b_ind)
end

@testset "float individual" begin
    f_ind = FloatIndividual(cfg)
    @test f_ind.fitness[1] == -Inf

    f_str = string(f_ind)
    @test typeof(f_str) == String

    f_ind2 = FloatIndividual(f_str)
    @test typeof(f_ind2) == FloatIndividual
    @test f_ind2.fitness[1] == -Inf
    @test all(f_ind.genes .== f_ind2.genes)

    test_individual(f_ind)
end


text:
export S²
export ComplexLine
export Cartesian
export Spherical
export Geographic

"""
    Represents a point in 2-sphere.
"""
abstract type S² end

"""
    Represents a point in the complex line.

field: z.
"""
struct ComplexLine <: S²
    z::Complex
end

"""
    Represents a point in the Cartesian coordinate system.

fields: x, y and z.
"""
struct Cartesian <: S²
    x::Float64
    y::Float64
    z::Float64
    Cartesian(x::Real, y::Real, z::Real) = new(float(x), float(y), float(z))
    Cartesian(a::Array) = begin
        @assert(length(a) == 3, "The length of the input vector must be exactly 3.")
        Cartesian(a...)
    end
    Cartesian(r::ℝ³) = Cartesian(vec(r))
end

"""
    Represents a point in the spherical coordinate system.

fields: r, ϕ and θ.
"""
struct Spherical <: S²
    r::Float64
    ϕ::Float64
    θ::Float64
    Spherical(r::Real, ϕ::Real, θ::Real) = begin
        @assert(r ≥ 0, "r must be non-negative.")
        @assert(0 ≤ ϕ ≤ 2π, "ϕ must be in [0, 2π].")
        @assert(0 ≤ θ ≤ π, "θ must be in [0, π].")
        new(float(r), float(ϕ), float(θ))
    end
    Spherical(a::Array) = begin
        @assert(length(a) == 3, "The length of the input vector must be exactly 3.")
        Spherical(a...)
    end
end

"""
    Represents a point in the geographic coordinate system.

fields: ϕ and θ.
"""
struct Geographic <: S²
    ϕ::Float64
    θ::Float64
    Geographic(ϕ::Real, θ::Real) = begin
        @assert(-π ≤ ϕ ≤ π, "ϕ must be in [-π, π].")
        @assert(-π/2 ≤ θ ≤ π/2, "θ must be in [-π/2, π/2].")
        new(float(ϕ), float(θ))
    end
    Geographic(a::Array) = begin
        @assert(length(a) == 2, "The length of the input vector must be exactly 2.")
        Geographic(a...)
    end
end

Base.real(cl::ComplexLine) = real(cl.z)
Base.imag(cl::ComplexLine) = imag(cl.z)
Base.abs(cl::ComplexLine) = abs(cl.z)
Base.angle(cl::ComplexLine) = angle(cl.z)
Base.conj(cl::ComplexLine) = conj(cl.z)

Base.vec(cl::ComplexLine) = [cl.z; conj(cl)]
Base.vec(c::Cartesian) = [c.x; c.y; c.z]
Base.vec(s::Spherical) = [s.r; s.ϕ; s.θ]
Base.vec(g::Geographic) = [g.ϕ; g.θ]

Base.isapprox(c1::Cartesian, c2::Cartesian) = isapprox(vec(c1), vec(c2))
Base.isapprox(s1::S², s2::S²) = isapprox(Cartesian(s1), Cartesian(s2))
Base.isapprox(a::Array{<:S²,N} where N,
              b::Array{<:S²,N} where N;
              atol::Float64 = TOLERANCE) = begin
    for (elementa, elementb) in zip(a, b)
        if isapprox(vec(elementa), vec(elementb), atol = atol) == false
            return false
        end
    end
    return true
end

ℝ³(c::Cartesian) = ℝ³(vec(c))
ℝ³(s::S²) = ℝ³(Cartesian(s))

norm(c::Cartesian) = norm(ℝ³(c))

ComplexLine(r3::ℝ³) = ComplexLine(Cartesian(r3))
ComplexLine(s::Spherical) = ComplexLine(cot(s.θ / 2) * exp(im * s.ϕ))
ComplexLine(c::Cartesian) = ComplexLine((c.x + im * c.y) / (1 - c.z))
ComplexLine(g::Geographic) = ComplexLine(Spherical(g))

Cartesian(cl::ComplexLine) = begin
    x, y = real(cl), imag(cl)
    d = x^2 + y^2 + 1
    Cartesian(2x / d, 2y / d, (d - 2) / d)
end
Cartesian(s::Spherical) = begin
    r, ϕ, θ = s.r, s.ϕ, s.θ
    Cartesian(r * sin(θ) * cos(ϕ), r * sin(θ) * sin(ϕ), r * cos(θ))
end
Cartesian(g::Geographic) = Cartesian(Spherical(g))

Spherical(cl::ComplexLine) = Spherical([1; angle(cl); 2acot(abs(cl))])
Spherical(c::Cartesian) = begin
    x, y, z = vec(c)
    if x > 0
        ϕ = atan(y / x)
    elseif y > 0
        ϕ = atan(y / x) + pi
    else
        ϕ = atan(y / x) - pi
    end
    r = norm(c)
    θ = asin(z / r)
    Spherical(r, ϕ, pi/2 - θ)
end
Spherical(g::Geographic) = Spherical(1, g.ϕ + pi, pi/2 - g.θ)
Spherical(r::ℝ³) = Spherical(Cartesian(r))

Geographic(g::Geographic) = g
Geographic(s::Spherical) = Geographic(s.ϕ - pi, pi/2 - s.θ)
Geographic(r::ComplexLine) = Geographic(Spherical(r))
Geographic(c::Cartesian) = Geographic(ComplexLine(c))
Geographic(r::ℝ³) = Geographic(Cartesian(r))


text:
using StandardMarketData
using HDF5Utils
using PandasLite
using Statistics
using Dates
using Random
using Test

cd(mktempdir())

Random.seed!(1234)
F, N, T = 2, 5, 100

特征名 = idxmap(string.("f", 1:F))
特征 = randn(Float32, F, N, T)
涨幅 = dropdims(mean(特征, dims = 1), dims = 1)
特征[:, 1:1, 1:3] .= NaN
特征[:, 2:2, 1:3] .= Inf
ti, Δt = DateTime(2019, 1, 1), Hour(1)
时间戳 = range(ti, step = Δt, length = T ÷ 2)
时间戳 = datetime2unix.(repeat(reshape(时间戳, 1, :), N, 2))
代码 = MLString{8}[string(2t <= T ? n : N + n) for n in 1:N, t in 1:T]
最新价 = 买1价 = 卖1价 = cumsum(涨幅, dims = 2)
手续费率 = fill(1f-4, N, T)
涨停 = 跌停 = zeros(Float32, N, T)
交易池 = ones(Float32, N, T)
data = Data(特征名, 特征, 涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 涨停, 跌停, 交易池)

savedata("test.h5", data)
@test data == reloaddata(data)
@test nticks(downsample(data, "5H")) == T ÷ 5
Sys.iswindows() || @test size(pivot(data)[1]) == (2N, T ÷ 2)
@test isapprox(pct_change(data, 5)[:, 10], sum(涨幅[:, 11:15], dims = 2))

@test nfeats(data) == F
@test ncodes(data) == N
@test nticks(data) == T
@test period(data) == Second(Δt).value
@test length(featnames(data)) == F

@test ncodes(data["1"]) == 1
data["20190101":"20190103"]
split(data, "20190110")

@uncol f1, f2 = data
@test isapprox(getfeat(data, r"f1$"), 特征[1, :, :], nans = true)
@test isequal(getfeats(data, r"f1$"), [特征[1, :, :]])
@test nfeats(dropfeats(data, r"f1$")) == F - 1
@test nfeats(keepfeats(data, r"f1$")) == 1
@test nfeats(keepcats(data, r"f1$")) == F
@test getcats(data) == featnames(data)

df = to_df(data)
@test isequal(df["f1"].values, vec(f1))
to_data(df, "tmp.h5")
@test all(isone, setcomm(data, 1).手续费率)
@test isempty(setpool(data, "f1"))

@test firstdate(data) == minimum(unix2date, 时间戳)
@test lastdate(data) == maximum(unix2date, 时间戳)
@test datespan(data) == (firstdate(data), lastdate(data))

@test ncodes(concat([data, data], dims = 1)) == 2N
@test nticks(concat([data, data], dims = 2)) == 2T

@test epochsof(data) == unique(时间戳)
@test datetimesof(data) == unix2datetime.(epochsof(data))
@test datesof(data) == unix2date.(epochsof(data))
@test Set(codesof(data)) == Set(unique(代码))

@test parsefreq("1T") == 60
@test str2date("20100101") == Date(2010, 01, 01)
@test int2unix(20100101) == datetime2unix(DateTime(2010, 01, 01))

@test nunique([1, 2, 2, 3]) == 3
@test sortednunique([1, 2, 3, 3, 4]) == 4
@test isempty(rolldata(data, "6M", "6M"))

if v"0.25" <= PandasLite.version() < v"1"
    extract_tsfresh_feats(to_df(data), shifts = ["10H"], horizon = "3H")
end

if !Sys.iswindows()
    df = DataFrame("code" => vec(代码), "close" => vec(最新价))
    df["high"] = df["low"] = df["open"] = df["close"]
    df["volume"] = rand(N * T)
    extract_talib_feats(df, "code")
end

if !isnothing(Sys.which("mpiexec"))
    x = data.特征
    x8, bin_edges = discretize(x, mpiarg = `-n 1`)
    x′ = undiscretize(x8, bin_edges)
    @test mean(x .- x′) do z
        abs(!isnan(z) * !isinf(z) * z)
    end < 0.02
end

# foo(x) = unsqueeze(x; dims = 2)
# x = rand(10, 10)
# @code_warntype foo(x)

# xs = [rand(10, 10) for i in 1:2]
# concatenate(xs, dims = 1)
# foo(xs) = concatenate(xs, dims = 1)
# foo(xs); @code_warntype foo(xs)


text:
class Test1 {
    field x
    def new() {}
}

def test(a,b,c) {
    a+1, b+2, c+3
}

def main {
    test1 = Test1.new()

    {test1.x, y, z} = test(1, 2, 3)
    ;print(test1.x)
    ;print(y)
    ;print(z)

    arr = [1, 2, 3, 4, 5]

    print((test1.x, y, z)[1])

}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
module ML_Hopf

# Add core packages to the enviroment
# add DifferentialEquations, DiffEqFlux, PGFPlotsX, MAT, Plots, NLsolve, StaticArrays, ForwardDiff, LinearAlgebra, Statistics

using DifferentialEquations
using LinearAlgebra
using DiffEqFlux
using PGFPlotsX
using LaTeXStrings
using MAT
using Plots
using NLsolve
using Statistics
using StaticArrays
using ForwardDiff

export generate_data, Array_chain

function flutter_eq_CBC(du, u, p, t) #Flutter equation of motion
    μ, Kp, Kd, Fa, Fω = p[1:5]
    c = p[6:end]
    ind1 = 1
    ind2 = 3
    theta = atan(u[ind1], u[ind2]) #theta is computed from two state variables
    r = c[1]
    h = c[2:end]
    nh = Int(length(h) / 2)
    for i in 1:nh
        r += h[i] * cos(i * theta) + h[i + nh] * sin(theta * i)
    end
    b = 0.15
    a = -0.5
    rho = 1.204
    x__alpha = 0.2340
    c__0 = 1
    c__1 = 0.1650
    c__2 = 0.0455
    c__3 = 0.335
    c__4 = 0.3
    c__alpha = 0.562766779889303
    c__h = 15.4430
    I__alpha = 0.1726
    k__alpha = 54.116182926744390
    k__h = 3.5294e+03
    m = 5.3
    m__T = 16.9
    U = μ

    MM = [
        b^2 * pi * rho+m__T -a * b^3 * pi * rho+b * m * x__alpha 0
        -a * b^3 * pi * rho+b * m * x__alpha I__alpha+pi * (0.1e1 / 0.8e1 + a^2) * rho * b^4 0
        0 0 1
    ]
    DD = [
        c__h+2 * pi * rho * b * U * (c__0 - c__1 - c__3) (1+(c__0 - c__1 - c__3) * (1 - 2 * a))*pi*rho*b^2*U 2*pi*rho*U^2*b*(c__1 * c__2+c__3 * c__4)
        -0.2e1*pi*(a+0.1e1 / 0.2e1)*rho*(b^2)*(c__0 - c__1-c__3)*U c__alpha+(0.1e1 / 0.2e1 - a) * (1 - (c__0 - c__1 - c__3) * (1 + 2 * a)) * pi * rho * (b^3) * U -0.2e1*pi*rho*(U^2)*(b^2)*(a+0.1e1 / 0.2e1)*(c__1 * c__2+c__3 * c__4)
        -1/b a-0.1e1 / 0.2e1 (c__2 + c__4) * U/b
    ]
    KK = [
        k__h 2*pi*rho*b*U^2*(c__0 - c__1-c__3) 2*pi*rho*U^3*c__2*c__4*(c__1+c__3)
        0 k__alpha-0.2e1 * pi * (a + 0.1e1 / 0.2e1) * rho * (c__0 - c__1 - c__3) * (b^2) * (U^2) -0.2e1*pi*rho*b*(U^3)*(a+0.1e1 / 0.2e1)*c__2*c__4*(c__1+c__3)
        0 -U/b c__2 * c__4 * U^2/b^2
    ]

    K1 = -inv(MM) * KK
    D1 = -inv(MM) * DD

    J1 = [0 1 0 0 0 0]
    J2 = [K1[1, 1] D1[1, 1] K1[1, 2] D1[1, 2] K1[1, 3] D1[1, 3]]
    J3 = [0 0 0 1 0 0]
    J4 = [K1[2, 1] D1[2, 1] K1[2, 2] D1[2, 2] K1[2, 3] D1[2, 3]]
    J5 = [0 0 0 0 0 1]
    J6 = [K1[3, 1] D1[3, 1] K1[3, 2] D1[3, 2] K1[3, 3] D1[3, 3]]

    J = [J1; J2; J3; J4; J5; J6]
    du_ = J * u
    #Control added on heave
    M2 = inv(MM)
    control = Kp * (r * cos(theta) - u[ind1]) + Kd * (r * sin(theta) - u[ind2])
    ka2 = 751.6
    ka3 = 5006
    nonlinear_h = -M2[1, 2] * (ka2 * u[3]^2 + ka3 * u[3]^3)
    nonlinear_theta = -M2[2, 2] * (ka2 * u[3]^2 + ka3 * u[3]^3)

    du[1] = du_[1]
    du[2] = du_[2] + nonlinear_h + M2[1, 1] * control + M2[1, 1] * Fa * sin(Fω * t)
    du[3] = du_[3]
    du[4] = du_[4] + nonlinear_theta + M2[2, 1] * control + M2[2, 1] * Fa * sin(Fω * t)
    du[5] = du_[5]
    return du[6] = du_[6]
end

function generate_data(vel_l, Vel, nh) # Training data
    #Generate training data
    u0 = Float32[2e-1, 0, 2e-1, 0, 0, 0]
    tol = 1e-7
    stol = 1e-8
    eq = flutter_eq_CBC
    rp = 5
    ind1 = 1
    ind2 = 3
    pp = [zeros(10) for i in 1:vel_l]
    AA = zeros(vel_l, Int(nh * 2 + 1))

    p_ = zeros(9)
    p_ = vcat(Vel[1], p_)
    tl = 10.0
    g = get_stable_LCO(p_, u0, tl, tol, eq, stol, rp, ind1, ind2, 0.0, 0.0, st)

    u₀ = mean(g.u[:, 1])
    v₀ = mean(g.u[:, 3])
    for i in 1:vel_l
        pp[i][1] = Vel[i]
        g = get_stable_LCO(pp[i], u0, tl, tol, eq, stol, rp, ind1, ind2, u₀, v₀, st)
        r = g.r
        t = g.t
        c = LS_harmonics(r, t, 1, nh).coeff
        AA[i, :] = c
    end
    cc = [
        LS_harmonics(
            get_sol(U[s_ind[i]], N, 1, 3).r, get_sol(U[s_ind[i]], N, 1, 3).t, 1, nh
        ).coeff for i in 1:length(s_ind)
    ]
    cc = hcat(cc)

    AA = transpose(AA)
    Al = AA
    AA = hcat(AA, cc)
    t_series = [
        Transpose(
            get_stable_LCO([Vel[i] transpose(zeros(9))], u0, tl, tol, eq, stol, rp, ind1, ind2, u₀, v₀, st).u[
                :, [1, 3]
            ],
        ) for i in 1:vel_l
    ]
    θ_series = [
        get_stable_LCO(
            [Vel[i] transpose(zeros(9))], u0, tl, tol, eq, stol, rp, ind1, ind2, u₀, v₀, st
        ).t for i in 1:vel_l
    ]
    θ = range(0; stop=2π, length=θ_l)
    coθ = cos.(θ)
    siθ = sin.(θ)
    return (data=AA, ts=t_series, d0=[u₀, v₀], data2=Al, coθ=coθ, siθ=siθ, theta_s=θ_series)
end

## Numerical model

function nf_dis(U₀, s, Vel, Vel2, coθ, siθ)
    del = Vel - U₀ * ones(length(Vel))
    del2 = Vel2 - U₀ * ones(length(Vel2))
    va2 = s * ones(length(Vel))
    va2_2 = s * ones(length(Vel2))
    s_amp = sqrt.(va2 / 2 + sqrt.(va2 .^ 2 + 4 * del) / 2)
    u_amp = sqrt.(va2_2 / 2 - sqrt.(va2_2 .^ 2 + 4 * del2) / 2)

    vl = [s_amp[i] * [coθ'; siθ'] for i in 1:length(Vel)]
    vl2 = [u_amp[i] * [coθ'; siθ'] for i in 1:length(Vel2)]
    return (v=vl, v2=vl2)
end

function f_coeff(vlT, Vel, u₀, v₀, nh, θ_l)
    Pr = zeros(2 * nh + 1, 0)
    for k in 1:length(Vel)
        z1 = vlT[k][1, :] - u₀ * ones(θ_l)
        z2 = vlT[k][2, :] - v₀ * ones(θ_l)
        theta = atan.(z2, z1)
        r = sqrt.(z1 .^ 2 + z2 .^ 2)
        tM = Array{Float64}(undef, 0, 2 * nh + 1)
        rr = Array{Float64}(undef, θ_l)
        for j in 1:θ_l
            tM1 = Array{Float64}(undef, 0, nh + 1)
            tM2 = Array{Float64}(undef, 0, nh)
            tM1_ = [cos(theta[j] * i) for i in 1:nh]
            tM2_ = [sin(theta[j] * i) for i in 1:nh]
            tM1_ = vcat(1, tM1_)
            tM1 = vcat(tM1, Transpose(tM1_))
            tM2 = vcat(tM2, Transpose(tM2_))
            tM_ = hcat(tM1, tM2)
            tM = vcat(tM, tM_)
        end
        MM = Transpose(tM) * tM
        rN = Transpose(tM) * r
        c = inv(MM) * rN
        Pr = hcat(Pr, c)
        Pr
    end
    return Pr
end

function predict_lt(θ_t, U₀, s, Vel, Vel2, uu0, scale_f_l, nh, θ_l, coθ, siθ) #predict the linear transformation
    nf = nf_dis(U₀, s, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    p1, p2, p3, p4, p5, p6 = θ_t[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    T = [p1 p3; p2 p4]
    vlT = [dis * norm(vl[i][:, 1])^2 + T * (vl[i]) for i in 1:length(Vel)]
    vlT2 = [dis * norm(vl2[i][:, 1])^2 + T * (vl2[i]) for i in 1:length(Vel2)]
    u₀ = uu0[1]
    v₀ = uu0[2]
    Pr = f_coeff(vlT, Vel, u₀, v₀, nh, θ_l)
    Pr2 = f_coeff(vlT2, Vel2, 0, 0, nh, θ_l)
    return hcat(Pr, Pr2)
end

function lt_pp(θ_t, U₀, s, Vel, Vel2, scale_f_l, θ_l, coθ, siθ) # This function gives phase portrait of the transformed system from the normal form
    nf = nf_dis(U₀, s, Vel, Vel2, coθ, siθ)
    vl = nf.v
    p1, p2, p3, p4, p5, p6 = θ_t[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    vlT = [dis * norm(vl[i][:, 1])^2 + T * (vl[i]) for i in 1:length(Vel)]
    return vlT
end

function Array_chain(gu, ann, p) # vectorized input-> vectorized neural net
    al = length(gu[1, :])
    AC = zeros(2, 0)
    for i in 1:al
        AC = hcat(AC, ann(gu[:, i], p))
    end
    return AC
end

function loss_lt(θ_t, U₀, s, Vel, Vel2, uu0, AA, scale_f_l, nh, θ_l, coθ, siθ)
    pred = predict_lt(θ_t, U₀, s, Vel, Vel2, uu0, scale_f_l, nh, θ_l, coθ, siθ)
    return sum(abs2, AA .- pred) # + 1e-5*sum(sum.(abs, params(ann)))
end

function predict_nt(θ_t, θ_, θ_l, scale_f_l, Vel, Vel2, scale_f, ann, nh, uu0, coθ, siθ)
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    np1, np2 = θ_t[1:2]

    pn = θ_t[3:end]
    T = [p1 p3; p2 p4]

    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2

    vlT = [
        dis * norm(vl[i][:, 1])^2 +
        T * (vl[i]) +
        Array_chain([vl[i]; (Vel[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel)
    ]
    vlT2 = [
        dis * norm(vl2[i][:, 1])^2 +
        T * (vl2[i]) +
        Array_chain([vl2[i]; (Vel2[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel2)
    ]

    u₀ = uu0[1]
    v₀ = uu0[2]

    Pr = f_coeff(vlT, Vel, u₀, v₀, nh, θ_l)
    Pr2 = f_coeff(vlT2, Vel2, 0, 0, nh, θ_l)
    return hcat(Pr, Pr2)
end

function lt_pp_n(θ_t, θ_, θ_l, scale_f_l, Vel, Vel2, scale_f, ann, coθ, siθ) # This function gives phase portrait of the transformed system from the normal form (stable LCO)
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    np1, np2 = θ_t[1:2]
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v

    pn = θ_t[3:end]
    vlT = [
        dis * norm(vl[i][:, 1])^2 +
        T * (vl[i]) +
        Array_chain([vl[i]; (Vel[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel)
    ]
    return vlT
end

function lt_pp_n_u(θ_t, θ_, θ_l, scale_f_l, Vel, Vel2, scale_f, ann, coθ, siθ) # This function gives phase portrait of the transformed system from the normal form (unstable LCO)
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    np1, np2 = θ_t[1:2]
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl2 = nf.v2
    pn = θ_t[3:end]

    vlT = [
        dis * norm(vl2[i][:, 1])^2 +
        T * (vl2[i]) +
        Array_chain([vl2[i]; (Vel2[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel2)
    ]
    return vlT
end

function loss_nt(θ_t, θ_, θ_l, scale_f_l, Vel, Vel2, scale_f, AA, ann, nh, uu0, coθ, siθ)
    pred = predict_nt(θ_t, θ_, θ_l, scale_f_l, Vel, Vel2, scale_f, ann, nh, uu0, coθ, siθ)
    return sum(abs2, AA .- pred)
end

function lt_b_dia(θ_t, ind, vel_l, ann, θ_l, coθ, siθ, scale_f_l, scale_f)
    p1, p2, p3, p4, p5, p6 = θ_t[1:6]
    T = [p1 p3; p2 p4]
    np1, np2 = θ_t[7:8]
    Vel = range(np1 - np2^2 / 4 + 1e-7; stop=np1, length=vel_l)
    nf = nf_dis(np1, np2, Vel, Vel, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    pn = θ_t[9:end]

    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    vlT = [
        dis * norm(vl[i][:, 1])^2 +
        T * (vl[i]) +
        Array_chain([vl[i]; (Vel[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel)
    ]
    vlT2 = [
        dis * norm(vl[i][:, 1])^2 +
        T * (vl2[i]) +
        Array_chain([vl2[i]; (Vel[i] - np1) * ones(1, θ_l)], ann, pn) / scale_f for
        i in 1:length(Vel)
    ]
    vlTas = [maximum(vlT[i][ind, :]) - minimum(vlT[i][ind, :]) for i in 1:length(Vel)]
    vlTau = [maximum(vlT2[i][ind, :]) - minimum(vlT2[i][ind, :]) for i in 1:length(Vel)]
    return (s=vlTas, u=vlTau, v=Vel)
end

## Experimental model

function LT(vl, T, np1, dis, Vel, pn1, ann_l, scale_f2)
    TT =
        dis * norm(vl[:, 1])^2 +
        (T + reshape(ann_l([norm(vl[:, 1]), Vel - np1], pn1), 2, 2) / scale_f2) * (vl)
    return TT
end

function predict_lt_nn(θ_t, ann_l, scale_f_l, scale_f2, θ_l, nh, Vel, Vel2, coθ, siθ)
    np1 = θ_t[end - 1]
    np2 = θ_t[end]
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    p1, p2, p3, p4, p5, p6 = θ_t[1:6]
    T = [p1 p3; p2 p4]

    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    pn1 = θ_t[7:(end - 2)]
    vlT = [LT(vl[i], T, np1, dis, Vel[i], pn1, ann_l, scale_f2) for i in 1:length(Vel)]
    vlT2 = [LT(vl2[i], T, np1, dis, Vel2[i], pn1, ann_l, scale_f2) for i in 1:length(Vel2)]

    Pr = f_coeff(vlT, Vel, 0, 0, nh, θ_l)
    Pr2 = f_coeff(vlT2, Vel2, 0, 0, nh, θ_l)
    PP = hcat(Pr, Pr2)
    return PP
end

function predict_lt_nn2(
    θ_t, ann_l, scale_f_l, scale_f2, U₀, s_, θ_l, nh, Vel, Vel2, coθ, siθ
) #predict the linear transformation
    np1 = U₀
    np2 = s_
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    p1, p2, p3, p4, p5, p6 = θ_t[1:6]
    T = [p1 p3; p2 p4]

    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    pn1 = θ_t[7:(end - 2)]
    vlT = [LT(vl[i], T, np1, dis, Vel[i], pn1, ann_l, scale_f2) for i in 1:length(Vel)]
    vlT2 = [LT(vl2[i], T, np1, dis, Vel2[i], pn1, ann_l, scale_f2) for i in 1:length(Vel2)]

    Pr = f_coeff(vlT, Vel, 0, 0, nh, θ_l)
    Pr2 = f_coeff(vlT2, Vel2, 0, 0, nh, θ_l)
    PP = hcat(Pr, Pr2)
    return PP
end

function loss_lt_nn(θ_t, ann_l, scale_f_l, scale_f2, AA, θ_l, nh, Vel, Vel2, coθ, siθ)
    pred = predict_lt_nn(θ_t, ann_l, scale_f_l, scale_f2, θ_l, nh, Vel, Vel2, coθ, siθ)
    return sum(abs2, AA .- pred)
end

function loss_lt_nn2(
    θ_t, ann_l, scale_f_l, scale_f2, U₀, s_, AA, θ_l, nh, Vel, Vel2, coθ, siθ
)
    pred = predict_lt_nn2(
        θ_t, ann_l, scale_f_l, scale_f2, U₀, s_, θ_l, nh, Vel, Vel2, coθ, siθ
    )
    return sum(abs2, AA .- pred)
end

function Nt(vl, T, Vel, dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann)
    vlT =
        LT(vl, T, np1, dis, Vel, pn1, ann_l, scale_f2) +
        Array_chain([vl; (Vel - np1) * ones(1, length(vl[1, :]))], ann, pn2) / scale_f
    return vlT
end

function predict_nt2(
    θ_t, θ_, ann_l, scale_f, scale_f2, scale_f_l, ann, nh, Vel, Vel2, θ_l, coθ, siθ
)
    np1 = θ_t[end - 1]
    np2 = θ_t[end]
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    pn1 = θ_[7:(end - 2)]
    pn2 = θ_t[1:(end - 2)]

    vlT = [
        Nt(vl[i], T, Vel[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel)
    ]
    vlT2 = [
        Nt(vl2[i], T, Vel2[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel2)
    ]

    Pr = f_coeff(vlT, Vel, 0, 0, nh, θ_l)
    Pr2 = f_coeff(vlT2, Vel2, 0, 0, nh, θ_l)
    return hcat(Pr, Pr2)
end

function lt_pp_n2(
    θ_t, θ_, ann_l, scale_f, scale_f2, scale_f_l, ann, Vel, Vel2, θ_l, coθ, siθ
) # This function gives phase portrait of the transformed system from the normal form (stable LCO)
    np1 = θ_t[end - 1]
    np2 = θ_t[end]
    nf = nf_dis(np1, np2, Vel, Vel2, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    dis = transpose([p5 * ones(θ_l) p6 * ones(θ_l)]) / scale_f_l
    pn1 = θ_[7:(end - 2)]
    pn2 = θ_t[1:(end - 1)]

    vlT = [
        Nt(vl[i], T, Vel[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel)
    ]
    vlT2 = [
        Nt(vl2[i], T, Vel2[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel2)
    ]

    return vcat(vlT, vlT2)
end

function loss_nt2(
    θ_t, θ_, ann_l, scale_f, scale_f2, scale_f_l, ann, AA, nh, Vel, Vel2, θ_l, coθ, siθ
) # Loss function
    pred = predict_nt2(
        θ_t, θ_, ann_l, scale_f, scale_f2, scale_f_l, ann, nh, Vel, Vel2, θ_l, coθ, siθ
    )
    return sum(abs2, AA .- pred)
end

function lt_b_dia2(ind, θ_n, θ_, ann_l, ann, scale_f_l, scale_f, scale_f2, coθ, siθ)
    vel_l = 300
    np1 = θ_n[end - 1]
    np2 = θ_n[end]
    Vel = range(np1 - np2^2 / 4 + 1e-7; stop=np1, length=vel_l)
    p1, p2, p3, p4, p5, p6 = θ_[1:6]
    T = [p1 p3; p2 p4]
    pn1 = θ_[7:(end - 2)]
    pn2 = θ_n[1:(end - 2)]

    nf = nf_dis(np1, np2, Vel, Vel, coθ, siθ)
    vl = nf.v
    vl2 = nf.v2
    dis = transpose([p5 * ones(length(vl)) p6 * ones(length(vl))]) / scale_f_l

    vlT = [
        Nt(vl[i], T, Vel[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel)
    ]
    vlT2 = [
        Nt(vl2[i], T, Vel[i], dis, pn1, np1, pn2, ann_l, scale_f, scale_f2, ann) for
        i in 1:length(Vel)
    ]

    vlTas = [maximum(vlT[i][ind, :]) - minimum(vlT[i][ind, :]) for i in 1:length(Vel)]
    vlTau = [maximum(vlT2[i][ind, :]) - minimum(vlT2[i][ind, :]) for i in 1:length(Vel)]
    return (s=vlTas, u=vlTau, v=Vel)
end

include("Numerical_Cont_Hopf_CBC.jl")

end