text:
export MinorBody, pluto, Pluto, pluto_barycenter, PlutoBarycenter

abstract type MinorBody <: CelestialBody end

const MINOR_BODY_NAMES = (
    "Ceres",
    "Pallas",
    "Vesta",
    "Psyche",
    "Lutetia",
    "Ida",
    "Eros",
    "Davida",
    "Gaspra",
    "Steins",
    "Itokawa",
    "Tempel1",
    "Borrelly",
)

for body in MINOR_BODY_NAMES
    str_sym = Symbol(body)
    con_sym = Symbol(lowercase(body))
    @eval begin
        struct $str_sym <: MinorBody end
        const $con_sym = $str_sym()
        parent(::$str_sym) = ssb
        export $str_sym, $con_sym
    end
end

struct PlutoBarycenter <: Barycenter end
const pluto_barycenter = PlutoBarycenter()
Base.show(io::IO, ::PlutoBarycenter) = print(io, "Pluto Barycenter")
parent(::PlutoBarycenter) = ssb
naifid(::PlutoBarycenter) = 9
from_naifid(::Val{9}) = pluto_barycenter

struct Pluto <: MinorBody end
const pluto = Pluto()
parent(::Pluto) = pluto_barycenter
naifid(::Pluto) = 999
from_naifid(::Val{999}) = pluto

naifid(::Ceres) = 2000001
from_naifid(::Val{2000001}) = ceres

naifid(::Pallas) = 2000002
from_naifid(::Val{2000002}) = pallas

naifid(::Vesta) = 2000004
from_naifid(::Val{2000004}) = vesta

naifid(::Psyche) = 2000016
from_naifid(::Val{2000016}) = psyche

naifid(::Lutetia) = 2000021
from_naifid(::Val{2000021}) = lutetia

naifid(::Ida) = 2431010
from_naifid(::Val{2431010}) = ida

naifid(::Eros) = 2000433
from_naifid(::Val{2000433}) = eros

naifid(::Davida) = 2000511
from_naifid(::Val{2000511}) = davida

naifid(::Gaspra) = 9511010
from_naifid(::Val{9511010}) = gaspra

naifid(::Steins) = 2002867
from_naifid(::Val{2002867}) = steins

naifid(::Itokawa) = 2025143
from_naifid(::Val{2025143}) = itokawa

naifid(::Tempel1) = 1000093
Base.show(io::IO, ::Tempel1) = print(io, "Tempel 1")
from_naifid(::Val{1000093}) = tempel1

naifid(::Borrelly) = 1000005
from_naifid(::Val{1000005}) = borrelly


text:
# Copied from https://github.com/yuyichao/FunctionWrappers.jl.
# Unfortunately, version 1.0.0 of FunctionWrappers doesn't work on Julia 1.0 due to
# https://github.com/yuyichao/FunctionWrappers.jl/issues/8 (caused by a bug in Julia Base)
# The code below is a copy of https://github.com/tkoolen/FunctionWrappers.jl/tree/tk/julia-0.7-quickfix,
# which works on 1.0 but uses deprecated Julia functionality.

# > Copyright (c) 2016: Yichao Yu
# > and other contributors:
# >
# > https://github.com/yuyichao/FunctionWrappers.jl/contributors
# >
# > Permission is hereby granted, free of charge, to any person obtaining
# > a copy of this software and associated documentation files (the
# > "Software"), to deal in the Software without restriction, including
# > without limitation the rights to use, copy, modify, merge, publish,
# > distribute, sublicense, and/or sell copies of the Software, and to
# > permit persons to whom the Software is furnished to do so, subject to
# > the following conditions:
# >
# > The above copyright notice and this permission notice shall be
# > included in all copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# > NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# > LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# > OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# > WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
module FunctionWrappersQuickFix

# Used to bypass NULL check
@inline function assume(v::Bool)
    Base.llvmcall(("declare void @llvm.assume(i1)",
                    """
                    %v = trunc i8 %0 to i1
                    call void @llvm.assume(i1 %v)
                    ret void
                    """), Cvoid, Tuple{Bool}, v)
end

is_singleton(@nospecialize(T)) = isdefined(T, :instance)

# Convert return type and generates cfunction signatures
Base.@pure map_rettype(T) =
    (isbitstype(T) || T === Any || is_singleton(T)) ? T : Ref{T}
Base.@pure function map_cfunc_argtype(T)
    if is_singleton(T)
        return Ref{T}
    end
    return (isbitstype(T) || T === Any) ? T : Ref{T}
end
Base.@pure function map_argtype(T)
    if is_singleton(T)
        return Any
    end
    return (isbitstype(T) || T === Any) ? T : Any
end
Base.@pure get_cfunc_argtype(Obj, Args) =
    Tuple{Ref{Obj}, (map_cfunc_argtype(Arg) for Arg in Args.parameters)...}

# Call wrapper since `cfunction` does not support non-function
# or closures
struct CallWrapper{Ret} <: Function end
(::CallWrapper{Ret})(f, args...) where {Ret} = convert(Ret, f(args...))

# Specialized wrapper for
for nargs in 0:128
    @eval function (::CallWrapper{Ret})(f, $((Symbol("arg", i) for i in 1:nargs)...)) where Ret
        convert(Ret, f($((Symbol("arg", i) for i in 1:nargs)...)))
    end
end

mutable struct FunctionWrapper{Ret,Args<:Tuple}
    ptr::Ptr{Cvoid}
    objptr::Ptr{Cvoid}
    obj
    objT

    function FunctionWrapper{Ret,Args}(obj::objT) where {Ret,Args,objT}
        objref = Base.cconvert(Ref{objT}, obj)
        # ptr = cfunction(
        #     CallWrapper{Ret}(), map_rettype(Ret),
        #     get_cfunc_argtype(objT, Args))
        # FIXME: use @cfunction (problem: it expects a literal tuple for the argument types)
        ptr = ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), CallWrapper{Ret}(), map_rettype(Ret), get_cfunc_argtype(objT, Args))
        new{Ret,Args}(ptr, Base.unsafe_convert(Ref{objT}, objref), objref, objT)
    end

    FunctionWrapper{Ret,Args}(obj::FunctionWrapper{Ret,Args}) where {Ret, Args} = obj
end

Base.convert(::Type{T}, obj) where {T<:FunctionWrapper} = T(obj)
Base.convert(::Type{T}, obj::T) where {T<:FunctionWrapper} = obj

@noinline function reinit_wrapper(f::FunctionWrapper{Ret,Args}) where {Ret,Args}
    objref = f.obj
    objT = f.objT
    # ptr = cfunction(CallWrapper{Ret}(), map_rettype(Ret),
    #                 get_cfunc_argtype(objT, Args))
    # FIXME: use @cfunction (problem: it expects a literal tuple for the argument types)
    ptr = ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), CallWrapper{Ret}(), map_rettype(Ret), get_cfunc_argtype(objT, Args))
    f.ptr = ptr
    f.objptr = Base.unsafe_convert(Ref{objT}, objref)
    return ptr
end

@generated function do_ccall(f::FunctionWrapper{Ret,Args}, args::Args) where {Ret,Args}
    # Has to be generated since the arguments type of `ccall` does not allow
    # anything other than tuple (i.e. `@pure` function doesn't work).
    quote
        Base.@_inline_meta
        ptr = f.ptr
        if ptr == C_NULL
            # For precompile support
            ptr = reinit_wrapper(f)
        end
        assume(ptr != C_NULL)
        objptr = f.objptr
        ccall(ptr, $(map_rettype(Ret)),
                (Ptr{Cvoid}, $((map_argtype(Arg) for Arg in Args.parameters)...)),
                objptr, $((:(args[$i]) for i in 1:length(Args.parameters))...))
    end
end

@inline (f::FunctionWrapper)(args...) = do_ccall(f, args)

# Testing only
const identityAnyAny = FunctionWrapper{Any,Tuple{Any}}(identity)

end

ext: jl


text:
const array_h = """
template<typename T, int N>
struct Array {
    T *data;
    const int dims[N];
    const int strides[N];
public:
    __device__ int length() {
        int n = dims[0];
        for (int i = 1; i < N; i++) n *= dims[i];
        return n;
    }
    __device__ T& operator[](const int idx) { return data[idx]; }
    __device__ T& operator()(const int idx) {
        int i = idx0 * strides[0];
        return data[i];
    }
    __device__ T& operator()(const int idx0, const int idx1) {
        int i = idx0 * strides[0];
        i += idx1 * strides[1];
        return data[i];
    }
    __device__ T& operator()(const int idx0, const int idx1, const int idx2) {
        int i = idx0 * strides[0];
        i += idx1 * strides[1];
        i += idx2 * strides[2];
        return data[i];
    }
    __device__ T& operator()(const int idx0, const int idx1, const int idx2, const int idx3) {
        int i = idx0 * strides[0];
        i += idx1 * strides[1];
        i += idx2 * strides[2];
        i += idx3 * strides[3];
        return data[i];
    }
    __device__ T& operator()(const int idx0, const int idx1, const int idx2, const int idx3, const int idx4) {
        int i = idx0 * strides[0];
        i += idx1 * strides[1];
        i += idx2 * strides[2];
        i += idx3 * strides[3];
        i += idx4 * strides[4];
        return data[i];
    }
    __device__ T& operator()(const int *idxs) {
        int i = 0;
        for (int k = 0; k < N; k++) {
            i += idxs[k] * strides[k];
        }
        return data[i];
    }
    __device__ void subindex(int *out, const int idx) {
        int cumdims[N];
        cumdims[0] = 1;
        for (int i = 1; i < N; i++) cumdims[i] = cumdims[i-1] * dims[i-1];

        for (int i = N-1; i >= 1; i--) {
            int k = idx / cumdims[i];
            out[i] = k;
            idx -= k * cumdims[i];
        }
        out[0] = idx;
        return;
    }
};
"""

immutable _CuArray{T,N}
    ptr::Ptr{T}
    dims::Ntuple{N,Cint}
    strides::Ntuple{N,Cint}
end

box{T}(x::CuArray{T}) = _CuArray(Ptr{T}(x), map(Cint,size(x)), map(Cint,strides(x)))
box{T}(x::CuSubArray{T}) = _CuArray(Ptr{T}(x), map(Cint,size(x)), map(Cint,strides(x)))


text:
export nlcg

"""
nlcg(f,df,x)

Nonlinear conjugate gradient method for solving min_x f(x)

"""
function nlcg(f::Function,df::Function,x::Vector;maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,
	lineSearch::Function=(f,df,fk,dfk,xk,pk)-> armijo(f,fk,dfk,xk,pk,maxIter=30) )

    his = zeros(maxIter,3)
    X = (storeInterm) ? zeros(length(x),maxIter) : []

	fk    = f(x)
    dfk   = df(x)
    dfOld = copy(dfk)
    pk    = -dfk

    i = 1; flag = -1;
    while i<=maxIter
        his[i,1:2] = [fk norm(dfk)]
        if storeInterm; X[:,i] = x; end;

        if (norm(dfk)<atol)
            flag=0
            his = his[1:i,:]
            break
        end

        # line search
        ak,his[i,3] = lineSearch(f,df,fk,dfk,x,pk)
        if out>0
            @printf "iter=%04d\t|f|=%1.2e\t|df|=%1.2e\tLS=%d\n" i his[i,1] his[i,2] his[i,3]
        end
        if his[i,3]==-1
            flag = -3
            his = his[1:i,:]
            break;
        end

        # update x and H
        x  += ak*pk
        fk  = f(x)
        dfk  = df(x)

        # beta  = dot(df, df-dfOld)/norm(dfOld)^2 # Polak-Ribiere
        beta  = norm(dfk)^2/norm(dfOld)^2 # Fletcher-Reeves
        pk    = -dfk  + beta*pk
        dfOld = copy(dfk)
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("nlcg iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif flag==-3
            println(@sprintf("nlcg stopped of a line search fail at iteration %d.",i))
        elseif out>1
            println(@sprintf("nlcg achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end


text:
"""
This functions returns an AxisArray of project parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_parameter_vector(
    structs::Vector, id::Symbol, parameter::Symbol)
    ids = getproperty.(structs, id)
    vals = getproperty.(structs, parameter)
    return AxisArrays.AxisArray(vals, ids)
end

"""
This functions returns an array of market parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_parameter_vector(
    structs::Vector{<:MarketCollection}, market::Symbol, parameter::Symbol)
    markets = getproperty.(structs, market)
    return getproperty.(markets, parameter)
end

"""
This functions returns a matrix of demand parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_demand_matrix(structs::Vector{<:MarketCollection}, market::Symbol)
    markets = getproperty.(structs, market)
    matrix = hcat(getproperty.(markets, :demand)...)
    return permutedims(matrix)
end

"""
This functions returns capacity market demand curve parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_capacity_demand_vectors(capmarkets::Vector{CapacityMarket})
   break_points =  getproperty.(capmarkets, :break_points)
   price_points =  getproperty.(capmarkets, :price_points)
   num_segments = length.(break_points) .- 1
   segment_size = [zeros(num_segments[i]) for i in 1:length(capmarkets)]
   segment_grad = [zeros(num_segments[i]) for i in 1:length(capmarkets)]
   for p in 1:length(capmarkets)
        for segment in 1:num_segments[p]
           segment_size[p][segment] = break_points[p][segment + 1] - break_points[p][segment]
           if segment_size[p][segment] != 0
                segment_grad[p][segment] = (price_points[p][segment + 1] - price_points[p][segment]) /  segment_size[p][segment]
           end
        end
    end

   return segment_size, segment_grad, price_points, num_segments
end

"""
This functions returns operating reserve demand curve parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_ORDC_vectors(reserveup_markets::Vector{ReserveUpMarket})
    break_points =  getproperty.(reserveup_markets, :break_points)
    price_points_raw =  getproperty.(reserveup_markets, :price_points)

    zones = AxisArrays.axisvalues(break_points[1])[1]

    num_segments = AxisArrays.AxisArray(Array{Int64, 2}(undef, length(zones), length(break_points)), zones, 1:length(break_points))
    segment_size = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(zones), length(break_points)), zones, 1:length(break_points))
    segment_grad = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(zones), length(break_points)), zones, 1:length(break_points))
    price_points = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(zones), length(break_points)), zones, 1:length(break_points))

    for z in zones
        for p in 1:length(break_points)
            num_segments[z, p] = length(break_points[p][z]) - 1
            segment_size[z, p] = zeros(num_segments[z, p])
            segment_grad[z, p] = zeros(num_segments[z, p])
            price_points[z, p] = price_points_raw[p][z]

            for segment in 1:num_segments[z, p]
                segment_size[z, p][segment] = break_points[p][z][segment + 1] - break_points[p][z][segment]
                if segment_size[z, p][segment] != 0
                    segment_grad[z, p][segment] = (price_points[z, p][segment + 1] - price_points[z, p][segment]) /  segment_size[z, p][segment]
                end
            end
        end
    end

    return segment_size, segment_grad, price_points, num_segments
 end

 """
 This functions returns capacity market demand curve parameters
 included in the actual clearing of capacity markets.
 """
function make_capacity_demand(capmarket::CapacityMarket)
    break_points =  getproperty(capmarket, :break_points)
    price_points =  getproperty(capmarket, :price_points)
    num_segments = length(break_points) - 1
    segment_size = zeros(num_segments)
    segment_grad = zeros(num_segments)
         for segment in 1:num_segments
            segment_size[segment] = break_points[segment + 1] - break_points[segment]
            segment_grad[segment] = (price_points[segment + 1] - price_points[segment]) /  segment_size[segment]
         end

    return segment_size, segment_grad, price_points
 end

 """
This functions returns capital cost multiplier curve parameters
included in CEM for price projection and endogeneous Economic Dispatch.
"""
function make_capital_cost_curve(options_by_type::Dict{String, Vector{String}},
                              annualized_cap_cost::AxisArrays.AxisArray{Float64, 2},
                              basecostunits::AxisArrays.AxisArray{Int64, 1},
                              maxnewoptions::AxisArrays.AxisArray{Int64, 1},
                              capital_cost_multiplier::Float64)

    invperiods = size(annualized_cap_cost, 2)
    types = collect(keys(options_by_type))

    num_segments = AxisArrays.AxisArray(round.(Int, 2 * ones(length(types), invperiods)), types, 1:invperiods)

    segment_size = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(types), invperiods), types, 1:invperiods)
    segment_grad = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(types), invperiods), types, 1:invperiods)
    costpoints = AxisArrays.AxisArray(Array{Vector{Float64}, 2}(undef, length(types), invperiods), types, 1:invperiods)

    for t in types
        for p in 1:invperiods
            segment_size[t, p] = zeros(num_segments[t, p])
            segment_grad[t, p] = zeros(num_segments[t, p])
            costpoints[t, p] = zeros(num_segments[t, p])

            if length(options_by_type[t]) >= 1
                total_max_new_options = sum(maxnewoptions[g] for g in options_by_type[t])
                g = options_by_type[t][1]
                segment_size[t, p][1] = basecostunits[g]
                segment_size[t, p][2] = total_max_new_options - basecostunits[g]
                costpoints[t, p] = [annualized_cap_cost[g, p], annualized_cap_cost[g, p]]
                segment_grad[t, p][1] = 0.0
                segment_grad[t, p][2] = costpoints[t, p][2] * capital_cost_multiplier
            end
        end
    end

    return segment_size, segment_grad, costpoints, num_segments
end

"""
This function updates the expected capacity factors of projects which are not in Option phase.
"""
function update_capacity_factors!(project::P,
                                 scenario_name::String,
                                 capacity_factors::Dict{String, Array{Float64, 2}}) where P <: Project{<: BuildPhase}
            for product in get_products(project)
                set_capacity_factors!(product, scenario_name, capacity_factors[get_name(project)])
            end
    return
end

"""
This function updates the expected capacity factors of Option projects.
"""
function update_capacity_factors!(project::P,
                                 scenario_name::String,
                                 capacity_factors::Dict{String, Array{Float64, 2}}) where P <: Project{Option}

            tech = get_tech(project)
            for product in get_products(project)
                set_capacity_factors!(product, scenario_name, capacity_factors["option_$(get_type(tech))_$(get_zone(tech))"])
            end
    return
end

"""
This function updates the expected accepted capacity percentage of projects which are not in Option phase.
"""
function update_capacity_accepted_perc!(project::P,
                                 scenario_name::String,
                                 capacity_accepted_perc::Dict{String, Array{Float64, 1}}) where P <: Project{<: BuildPhase}
            for product in get_products(project)
                set_accepted_perc!(product, scenario_name, capacity_accepted_perc[get_name(project)])
            end
    return
end

"""
This function updates the expected accepted capacity percentage of Option projects.
"""
function update_capacity_accepted_perc!(project::P,
                                 scenario_name::String,
                                 capacity_accepted_perc::Dict{String, Array{Float64, 1}}) where P <: Project{Option}

            tech = get_tech(project)
            for product in get_products(project)
                set_accepted_perc!(product, scenario_name, capacity_accepted_perc["option_$(get_type(tech))_$(get_zone(tech))"])
            end
    return
end

"""
This function does nothing if the product is not of Capacity type.
"""
function update_initial_capacity_revenues!(project::P,
                                 product::T,
                                 initial_capacity_prices::Vector{Float64},
                                 year::Int64) where {P <: Project{<: BuildPhase}, T <: Product}
    return
end

"""
This function updates the capacity market revenues of project at the start of the simulation.
"""
function update_initial_capacity_revenues!(project::P,
                                 product::Capacity,
                                 initial_capacity_prices::Vector{Float64},
                                 year::Int64) where P <: Project{<: BuildPhase}
    capacity_revenue = initial_capacity_prices[year] * get_maxcap(project) * get_project_derating(project)

    finance_data = get_finance_data(project)
    set_realized_profit!(get_finance_data(project),
                    get_name(product),
                    year,
                    capacity_revenue)

    for scenario_name in keys(get_scenario_profit(finance_data))
                update_forward_profit!(product, finance_data, scenario_name, year, capacity_revenue)
            end

    return
end

"""
This function does nothing if the product is not of Capacity or REC type.
"""
function update_bid!(product::T,
                     capacity_market_bid::Float64,
                     rec_market_bid::Float64,
                     energy_production::Float64) where T <: Product
    return
end

"""
This function updates the Capacity market bid of projects.
"""
function update_bid!(product::Capacity,
                     capacity_market_bid::Float64,
                     rec_market_bid::Float64,
                     energy_production::Float64)

    set_capacity_bid!(product, capacity_market_bid)

    return
end

"""
This function updates the REC market bid of projects.
"""
function update_bid!(product::REC,
                     capacity_market_bid::Float64,
                     rec_market_bid::Float64,
                     energy_production::Float64)

    set_rec_bid!(product, rec_market_bid)
    set_rec_certificates!(product, energy_production)
    return
end


text:
using BilliardModels

include("StatisticsAccumulator.jl")

function run_lorentz(num_particles, num_collisions, delta_t=5.0)

    billiard_table = Sinai_billiard(0.354, true, true)  # periodic in x and y

    stats = StatisticsObject[]

    num_stats_objects = 10
    for i in 1:num_stats_objects+1
        push!(stats, StatisticsObject())  # use all times as origin
    end

    # Divide statistics up according to first free path

    skip = 5

    for i in 1:num_particles
        x0, v0 = initial_condition(billiard_table, -.5, .5, -.5, .5)
        l = Vector2D(0, 0)
        p = ParticleOnLattice(x0, v0, l)

        xs, ls, free_paths = billiard_dynamics_on_lattice(p, billiard_table, num_collisions)
        positions, times = continuous_time(xs, ls, free_paths, delta_t)

        first_free_path = free_paths[1]

        if first_free_path < num_stats_objects
            stats_object = stats[ceil(first_free_path)]
        else
            stats_object = stats[end]
        end

        # "all" statistics: use each (or nearly each) time and position as a new origin
        for i in 1:skip:length(times)
            t00 = times[i]
            x00 = positions[i][1]  # x component

            for j in i+1:length(times) #skip:length(times)

                time_displacement = times[j] - t00
                displacement = positions[j][1] - x00

                #r = norm(displacement)
                #r = displacement
                add_data!(stats_object, time_displacement, displacement)
            end
        end

    end

    map(get_statistics, stats)

end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
function sqlite3_errmsg()
    return ccall( (:sqlite3_errmsg, sqlite3_lib),
                  Ptr{Uint8}, () )
end

function sqlite3_errmsg(db::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_errmsg, sqlite3_lib),
                  Ptr{Uint8}, (Ptr{Void},), db )
end

function sqlite3_open(file::AbstractString, handle::Vector{Ptr{Void}})
    return ccall( (:sqlite3_open, sqlite3_lib),
                  Cint, (Ptr{Uint8}, Ptr{Void}),
                  file, handle )
end

function sqlite3_open16(file::UTF16String, handle::Array{Ptr{Void},1})
    return ccall( (:sqlite3_open16, sqlite3_lib),
                  Cint, (Ptr{Uint16}, Ptr{Void}),
                  file, handle )
end

function sqlite3_close(handle::Ptr{Void})
    @NULLCHECK handle
    return ccall( (:sqlite3_close, sqlite3_lib),
                  Cint, (Ptr{Void},), handle )
end

function sqlite3_close_v2(handle::Ptr{Void})
    @NULLCHECK handle
    try
        return ccall( (:sqlite3_close_v2, sqlite3_lib),
                      Cint, (Ptr{Void},), handle )
    catch
        # Older versions of the library don't have this, abort to other close
        warn("sqlite3_close_v2 not available.")
        sqlite3_close(handle)
    end
end

function sqlite3_next_stmt(db::Ptr{Void}, stmt::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_next_stmt, sqlite3_lib),
                  Ptr{Void}, (Ptr{Void}, Ptr{Void}),
                  db, stmt )
end

function sqlite3_prepare_v2(handle::Ptr{Void}, query::AbstractString,
                            stmt::Array{Ptr{Void},1}, unused::Array{Ptr{Void},1})
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare_v2, sqlite3_lib),
                  Cint, (Ptr{Void}, Ptr{Uint8}, Cint, Ptr{Void}, Ptr{Void}),
                  handle, query, sizeof(query), stmt, unused )
end

function sqlite3_prepare16_v2(handle::Ptr{Void}, query::AbstractString,
                              stmt::Array{Ptr{Void},1}, unused::Array{Ptr{Void},1})
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare16_v2, sqlite3_lib),
                  Cint, (Ptr{Void}, Ptr{Uint16}, Cint, Ptr{Void}, Ptr{Void}),
                  handle, query, sizeof(query), stmt, unused )
end

function sqlite3_finalize(stmt::Ptr{Void})
    @NULLCHECK stmt
    return ccall( (:sqlite3_finalize, sqlite3_lib),
                  Cint, (Ptr{Void},), stmt )
end

function sqlite3_step(stmt::Ptr{Void})
    @NULLCHECK stmt
    return ccall( (:sqlite3_step, sqlite3_lib),
                  Cint, (Ptr{Void},), stmt )
end

function sqlite3_column_count(stmt::Ptr{Void})
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_count, sqlite3_lib),
                  Cint, (Ptr{Void},), stmt )
end

function sqlite3_column_type(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_type, sqlite3_lib),
                  Cint, (Ptr{Void}, Cint),
                  stmt, col )
end

function sqlite3_column_name(stmt::Ptr{Void},n::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_name, sqlite3_lib),
                  Ptr{Uint8}, (Ptr{Void}, Cint),
                  stmt, n )
end

function sqlite3_reset(stmt::Ptr{Void})
    @NULLCHECK stmt
    return ccall( (:sqlite3_reset, sqlite3_lib),
                  Cint, (Ptr{Void},), stmt )
end

function sqlite3_column_blob(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_blob, sqlite3_lib),
                  Ptr{Void}, (Ptr{Void}, Cint),
                  stmt, col )
end
function sqlite3_column_bytes(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_bytes, sqlite3_lib),
                  Cint, (Ptr{Void}, Cint),
                  stmt, col )
end
function sqlite3_column_bytes16(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_bytes16, sqlite3_lib),
                  Cint, (Ptr{Void}, Cint),
                  stmt, col )
end
@inline function sqlite3_column_double(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_double, sqlite3_lib),
                  Cdouble, (Ptr{Void}, Cint),
                  stmt, col )
end
function sqlite3_column_int(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_int, sqlite3_lib),
                  Cint, (Ptr{Void}, Cint),
                  stmt, col )
end
@inline function sqlite3_column_int64(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_int64, sqlite3_lib),
                  Clonglong, (Ptr{Void}, Cint),
                  stmt, col )
end
function sqlite3_column_text(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_text, sqlite3_lib),
                  Ptr{Uint8}, (Ptr{Void}, Cint),
                  stmt, col )
end
function sqlite3_column_text16(stmt::Ptr{Void},col::Int)
    @NULLCHECK stmt
    return ccall( (:sqlite3_column_text16, sqlite3_lib),
                  Ptr{Void}, (Ptr{Void}, Cint),
                  stmt, col )
end

# SQLITE_API void sqlite3_result_double(sqlite3_context*, double);
function sqlite3_result_double(context::Ptr{Void},value::Float64)
    return ccall( (:sqlite3_result_double, sqlite3_lib),
                  Void, (Ptr{Void}, Float64),
                  context, value )
end
# SQLITE_API void sqlite3_result_error(sqlite3_context*, const char*, int)
function sqlite3_result_error(context::Ptr{Void},msg::AbstractString)
    return ccall( (:sqlite3_result_error, sqlite3_lib),
                  Void, (Ptr{Void}, Ptr{Uint8}, Cint),
                  context, value, sizeof(msg) + 1 )
end
# SQLITE_API void sqlite3_result_error16(sqlite3_context*, const void*, int)
function sqlite3_result_error(context::Ptr{Void},msg::UTF16String)
    return ccall( (:sqlite3_result_error16, sqlite3_lib),
                  Void, (Ptr{Void}, Ptr{Uint16}, Cint),
                  context, value, sizeof(msg) + 1 )
end
# SQLITE_API void sqlite3_result_int(sqlite3_context*, int);
function sqlite3_result_int(context::Ptr{Void},value::Int32)
    return ccall( (:sqlite3_result_int, sqlite3_lib),
                  Void, (Ptr{Void}, Int32),
                  context, value )
end
# SQLITE_API void sqlite3_result_int64(sqlite3_context*, sqlite3_int64);
function sqlite3_result_int64(context::Ptr{Void},value::Int64)
    return ccall( (:sqlite3_result_int64, sqlite3_lib),
                  Void, (Ptr{Void}, Int64),
                  context, value )
end
# SQLITE_API void sqlite3_result_null(sqlite3_context*);
function sqlite3_result_null(context::Ptr{Void})
    return ccall( (:sqlite3_result_null, sqlite3_lib),
                  Void, (Ptr{Void},),
                  context )
end
# SQLITE_API void sqlite3_result_text(sqlite3_context*, const char*, int n, void(*)(void*));
function sqlite3_result_text(context::Ptr{Void},value::AbstractString)
    return ccall( (:sqlite3_result_text, sqlite3_lib),
                  Void, (Ptr{Void}, Ptr{Uint8}, Cint, Ptr{Void}),
                  context, value, sizeof(value) + 1, SQLITE_TRANSIENT )
end
# SQLITE_API void sqlite3_result_text16(sqlite3_context*, const void*, int, void(*)(void*));
function sqlite3_result_text16(context::Ptr{Void},value::UTF16String)
    return ccall( (:sqlite3_result_text, sqlite3_lib),
                  Void, (Ptr{Void}, Ptr{Uint16}, Cint, Ptr{Void}),
                  context, value, sizeof(value) + 1, SQLITE_TRANSIENT )
end
# SQLITE_API void sqlite3_result_blob(sqlite3_context*, const void*, int n, void(*)(void*));
function sqlite3_result_blob(context::Ptr{Void},value)
    return ccall( (:sqlite3_result_blob, sqlite3_lib),
                  Void, (Ptr{Void}, Ptr{Uint8}, Cint,Ptr{Void}),
                  context, value, sizeof(value), SQLITE_TRANSIENT )
end

function sqlite3_initialize()
    return ccall( (:sqlite3_initialize, sqlite3_lib),
        Cint, (),
        )
end
function sqlite3_shutdown()
    return ccall( (:sqlite3_shutdown, sqlite3_lib),
        Cint, (),
        )
end

sqlite3_types = Dict{Int, DataType}(
    1 => Int,
    2 => Float64,
    3 => UTF8String,
    4 => Any,
    5 => Nullable,
)


text:
#=
language_binding:
- Julia version:
- Author: thinkAmi
- Date: 2019-09-17
=#

using HDF5

h5open("output/h5py/binding.h5", "r") do f
    dataset = read(f, "hello")
    print(dataset[1])  # Juliaのindexは1始まり
    # => ワールド
end


text:
using DataFrames, CSV, Plots

# julia --project cli/main.jl mcmc "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Linear(nutrient) + Noise" --data test/testin/input_pair_109.tsv --output test/testin/mcmc_109.tsv --samples 100 --burnin 50 --thin 2
# julia --project cli/main.jl mcmc "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Noise" --data test/testin/input_pair_109.tsv --output test/testout/mcmc2_109.tsv --samples 100 --burnin 50 --thin 2
# julia --project cli/main.jl select --mcmc test/testin/mcmc_109.tsv --mcmc2 test/testout/mcmc2_109.tsv > test/testout/select_109.txt
# julia --project cli/main.jl predict "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Linear(nutrient) + Noise" --data test/testin/input_pair_109.tsv --mcmc test/testin/mcmc_109.tsv --at "nutrient=-5:0.1:5;PersonID=0;StoolPairs=0" --output test/testout/prediction_109.tsv

pred = CSV.File("test/testout/prediction_109.tsv") |> DataFrame
data = CSV.File("test/testin/input_pair_109.tsv") |> DataFrame

x = pred.nutrient
ymu = pred.ymu
y95 = pred.yQ950
y05 = pred.yQ050

scatter(data.bug, data.nutrient)
plot!(x, ymu, ribbon=(y95 .- ymu, ymu .- y05))
xlims!(extrema(data.bug))

describe(data.bug)

# julia --project cli/main.jl mcmc "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Linear(nutrient) + Noise" --data test/testin/input_pair_3206.tsv --output test/testin/mcmc_3206.tsv --samples 100 --burnin 50 --thin 2
# julia --project cli/main.jl mcmc "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Noise" --data test/testin/input_pair_3206.tsv --output test/testout/mcmc2_3206.tsv --samples 100 --burnin 50 --thin 2
# julia --project cli/main.jl select --mcmc test/testin/mcmc_3206.tsv --mcmc2 test/testout/mcmc2_3206.tsv > test/testout/select_3206.txt
# julia --project cli/main.jl predict "bug :~| Cat(PersonID) * Cat(StoolPairs) + Cat(PersonID) + Linear(nutrient) + Noise" --data test/testin/input_pair_3206.tsv --mcmc test/testin/mcmc_3206.tsv --at "nutrient=-5:0.1:5;PersonID=0;StoolPairs=0" --output test/testout/prediction_3206.tsv

pred = CSV.File("test/testout/prediction_3206.tsv") |> DataFrame
data = CSV.File("test/testin/input_pair_3206.tsv") |> DataFrame

x = pred.nutrient
ymu = pred.ymu
y95 = pred.yQ950
y05 = pred.yQ050

scatter(data.bug, data.nutrient)
plot!(x, ymu, ribbon=(y95 .- ymu, ymu .- y05))
xlims!(extrema(data.bug))

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
"""
    AxesConstraint

Abstract type for define mask's constraint. return a tuple of concrete subtype of AxesConstraint.
"""
abstract type AxesConstraint end

struct All1Constraint <: AxesConstraint
    from::Int
    n::Int
end

struct DimConstraint <: AxesConstraint
    dim::Int
    val::Int
end

struct NDimConstraint <: AxesConstraint
    n::Int
    least::Bool
end
NDimConstraint(n) = NDimConstraint(n, false)

@noinline thrdm(s) = throw(DimensionMismatch("arrays could not be broadcast to a common size; $s"))

axislength(x::Integer) = x
axislength(x) = length(x)

function check_constraint(c::All1Constraint, x)
    foreach(c.from:length(x)) do i
        if axislength(x[i]) != 1
            thrdm("mask require $i-th dimension to be 1, but get $(axislength(x[i]))")
        end
    end
    return x
end

function check_constraint(c::DimConstraint, x)
    if c.dim > 0
        dim = c.dim
        xi = length(x) ≥ dim ? axislength(x[dim]) : 1
        xi != c.val && thrdm("mask require $dim-th dimension to be $(c.val), but get $xi")
    else
        dim = length(x) + c.dim + 1
        xi = axislength(x[dim])
        xi != c.val && thrdm("mask require $dim-th dimension to be $(c.val), but get $xi")
    end
    return x
end

function check_constraint(c::NDimConstraint, x)
    cmp = c.least ? (≥) : (==)
    cmp(length(x), c.n) || thrdm("mask require ndims(A) $(cmp) $(c.n)")
    return x
end

check_constraint(cs::Tuple{}, x) = x
check_constraint(cs::Tuple, x) = check_constraint(Base.tail(cs), check_constraint(cs[1], x))

Base.iterate(c::AxesConstraint) = (c, nothing)
Base.iterate(c::AxesConstraint, ::Nothing) = nothing

@inline merge_constraint(cs1, cs2, cs...) = merge_constraint(merge_constraint(cs1, cs2), cs...)
@inline merge_constraint(::Tuple{}, ::Tuple{}) = ()
@inline merge_constraint(cs1::Tuple, ::Tuple{}) = cs1
@inline merge_constraint(::Tuple{}, cs2::Tuple) = cs2
@inline function merge_constraint(cs1::Tuple, cs2::Tuple)
    ndc = _merge_ndim_c(cs1[1], cs2[1])
    dc = _merge_dim_c(Base.tails(cs1, cs2)...)
    return (ndc, dc...)
end

_merge_dim_c(::Tuple{},::Tuple{}) = ()
_merge_dim_c(a::Tuple, ::Tuple{}) = a
_merge_dim_c(::Tuple{}, b::Tuple) = b
_merge_dim_c(a::Tuple{All1Constraint}, b::Tuple{All1Constraint}) = All1Constraint(min(a.from, b.from), max(a.n, b.n))

_merge_post_dim_c(::Tuple{}, ::Tuple{}) = ()
_merge_post_dim_c(::Tuple{}, b::Tuple) = b
_merge_post_dim_c(a::Tuple, ::Tuple{}) = a
function _merge_post_dim_c(a::Tuple, b::Tuple)::Tuple{Vararg{DimConstraint}}
    a1 = a[1]
    b1 = b[1]
    if a1.dim < b1.dim
        return (a1, _merge_post_dim_c(Base.tail(a), b)...)
    elseif a1.dim > b1.dim
        return (b1, _merge_post_dim_c(a, Base.tail(b))...)
    else
        a1.val != b1.val && thrdm("mask require $(a1.dim)-th dimension to be both $(a1.val) and $(b1.val)")
        return (a1, _merge_post_dim_c(Base.tails(a, b)...)...)
    end
end

function _merge_dim_c(a::Tuple{DimConstraint, Vararg{DimConstraint}}, b::Tuple{DimConstraint, Vararg{DimConstraint}})
    a1 = a[1]
    b1 = b[1]
    if (a1.dim > 0) == (b1.dim > 0)
        return _merge_post_dim_c(a, b)
    else
        if a1.dim < 0
            n = b[end].dim
            return _merge_post_dim_c(ntuple(i->DimConstraint(a[i].dim+n+1, a[i].val), length(a)), b)
        else
            n = a[end].dim
            return _merge_post_dim_c(a, ntuple(i->DimConstraint(b[i].dim+n+1, b[i].val), length(b)))
        end
    end
end

_merge_dim_c(a::Tuple{DimConstraint, Vararg{DimConstraint}}, b::Tuple{All1Constraint}) = _merge_dim_c(b, a)
function _merge_dim_c(a::Tuple{All1Constraint}, b::Tuple{DimConstraint, Vararg{DimConstraint}})
    c = a[]
    offset = c.from - 1
    return _merge_dim_c(ntuple(i->DimConstraint(i+offset, 1), b[end].dim - offset), b)
end

@inline function _merge_ndim_c(a, b)
    a.least || b.least || thrdm("mask require both ndims(A) == $(a.n) and ndim(A) == $(b.n)")
    n = max(a.n, b.n)
    a.least || a.n == n || thrdm("mask require both ndims(A) == $(a.n) and ndims(A) ≥ $(b.n)")
    b.least || b.n == n || thrdm("mask require both ndims(A) == $(b.n) and ndims(A) ≥ $(a.n)")
    return NDimConstraint(n, a.least == b.least)
end


text:
export
    fit_lane,
    move_along_curve,
    get_curvePt_cec_per

#####under developing

curveModel(x, p) = p[1]*x.^3+p[2]*x.^2+p[3]*x+p[4]

function fit_lane(model::MultiPtsDriver, scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway, ego_index::Int)
    sampleDistance=25.0;
    sampleInterval=1.0;
    #model(x, p) = p[1]*x^3+p[2]*x^2+p[3]*x+p[4]
    xdata=[];
    ydata=[];
    distance=0.0;
    veh = scene[ego_index]
    roadind=veh.state.posF.roadind
    #frenet=Frenet(roadind, roadway[roadind].s,0,0)
    mode=1
    x=-Inf
    y=-Inf
    while distance<=sampleDistance
        roadind_n=move_along(roadind, roadway, distance)
        frenet=Frenet(roadind_n, roadway[roadind_n].s,0.0,0.0)
        pos=get_posG(frenet,roadway)
        if pos.x==x #repeating on x, change to mode 2
            mode=2
        end
        x=pos.x
        y=pos.y
        push!(xdata,pos.x)
        push!(ydata,pos.y)
        distance=distance+sampleInterval
    end
    #println("xdata : ",xdata)
    #println("ydata : ",ydata)
    if mode==1
        fit = curve_fit(curveModel, xdata, ydata, [0.5, 0.5, 0.5, 0.5])
        return fit,xdata[1],ydata[1],xdata[2]>xdata[1],mode
    else
        #println("change to mode 2!")
        fit = curve_fit(curveModel, ydata, xdata, [0.5, 0.5, 0.5, 0.5])
        return fit,xdata[1],ydata[1],ydata[2]>ydata[1],mode
    end
end

function move_along_curve(fit::LsqFit.LsqFitResult,x::Float64,y::Float64,Δs::Float64,direction::Bool,mode::Int)
    p=fit.param
    x_n=x
    y_n=y
    if mode==1
        for i=1:10
            slope=3*p[1]*x_n^2+2*p[2]*x_n+p[3]
            dx=Δs/sqrt(slope^2+1)/10
            if direction==1
                x_n=x_n+dx
            else
                x_n=x_n-dx
            end
            y_n=p[1]*x_n^3+p[2]*x_n^2+p[3]*x_n+p[4]
        end
    else
        for i=1:10
            slope=3*p[1]*y_n^2+2*p[2]*y_n+p[3]
            dy=Δs/sqrt(slope^2+1)/10
            if direction==1
                y_n=y_n+dy
            else
                y_n=y_n-dy
            end
            x_n=p[1]*y_n^3+p[2]*y_n^2+p[3]*y_n+p[4]
        end
    end
    return x_n,y_n
end

function get_curvePt_vec_per(fit::LsqFit.LsqFitResult,x::Float64,y::Float64,turnDirection::Int,direction::Bool,mode::Int)
    p=fit.param
    if mode==1
        slope=3*p[1]*x^2+2*p[2]*x+p[3]
        curvePtVecTan=[1/sqrt(slope^2+1);slope/sqrt(slope^2+1)];
        if  direction==0
            curvePtVecTan=-curvePtVecTan
        end
        if turnDirection==1 #turn left
            curvePtVecPer=[0 -1;1 0]* curvePtVecTan
        else #turn right
            curvePtVecPer=[0 1;-1 0]* curvePtVecTan
        end
    else
        slope=3*p[1]*y^2+2*p[2]*y+p[3]
        curvePtVecTan=[slope/sqrt(slope^2+1);1/sqrt(slope^2+1)];
        if  direction==0
            curvePtVecTan=-curvePtVecTan
        end
        if turnDirection==1 #turn left
            curvePtVecPer=[0 -1;1 0]* curvePtVecTan
        else #turn right
            curvePtVecPer=[0 1;-1 0]* curvePtVecTan
        end
    end
    return curvePtVecPer
end

function set_lane_changing_pts!(model::MultiPtsDriver, scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway, ego_index::Int)
    #println("enter function")
    #laneFunction=fitLane(model, scene, roadway, ego_index)
    veh = scene[ego_index]
    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    steer=model.siji.steer
    steerdotmax=model.steerdotmax
    steermax=model.steermax
    accmax=model.accmax
    v = veh.state.v
    roadind=veh.state.posF.roadind
    Δt = model.C.Δt

    laneNum=roadind.tag.lane
    transitionStep=0

    frenet=Frenet(roadind, roadway[roadind].s,0.0,0.0)
    pos=get_posG(frenet,roadway)
    fit,xstart,ystart,direction,mode=fit_lane(model::DriverModel, scene::Scene, roadway::Roadway, ego_index::Int)

    transitionMagnitude=0.2
    offSetDistance=sqrt((x-xstart)^2+(y-ystart)^2)
    offSetVec=[x-xstart;y-ystart]

    v_desire=v

    v_desire=model.v_desire

    if model.laneNum_desire==laneNum
        curvePtVecPer=get_curvePt_vec_per(fit,x,y,1,direction,mode) #curvePtVecPer pointing to left
        offSetDistance=sqrt((x-xstart)^2+(y-ystart)^2)
        if (curvePtVecPer'*offSetVec)[1]>0
            turnDirection=1
        else
            turnDirection=-1
        end
        i=1
        x_n=x
        y_n=y
        while offSetDistance-transitionMagnitude>0 && i<=20
            v_desire_n=v_desire+model.acc_desire*(Δt*2)
            v_desire_n=clamp(v_desire_n,model.v_min,model.v_max)
            Δs=(v_desire+v_desire_n)*(Δt*2)/2
            v_desire=v_desire_n

            if i==6
                model.v_desire=v_desire
            end

            offSetDistance=offSetDistance-transitionMagnitude
            x_n,y_n=move_along_curve(fit,x_n,y_n,Δs,direction,mode)
            curvePtVecPer=get_curvePt_vec_per(fit,x_n,y_n,turnDirection,direction,mode)
            pos=[x_n;y_n]+offSetDistance*curvePtVecPer
            model.Pts[1,i]=pos[1]
            model.Pts[2,i]=pos[2]
            i=i+1
        end
        while offSetDistance-transitionMagnitude<=0 && i<=20
            v_desire_n=v_desire+model.acc_desire*(Δt*2)
            v_desire_n=clamp(v_desire_n,model.v_min,model.v_max)
            Δs=(v_desire+v_desire_n)*(Δt*2)/2
            v_desire=v_desire_n

            if i==6
                model.v_desire=v_desire
            end

            x_n,y_n=move_along_curve(fit,x_n,y_n,Δs,direction,mode)
            pos=[x_n;y_n]
            model.Pts[1,i]=pos[1]
            model.Pts[2,i]=pos[2]
            i=i+1
        end
    else
        if model.laneNum_desire>laneNum #turn to left lane
            turnDirection=1
            #println("turn left")
        else #turn to left right
            turnDirection=-1
            #println("right")
        end
        curvePtVecPer=curvePtVecPer=get_curvePt_vec_per(fit,x,y,turnDirection,direction,mode)
        if (curvePtVecPer'*offSetVec)[1]>0
            offSetDistance=sqrt((x-xstart)^2+(y-ystart)^2)
        else
            offSetDistance=-sqrt((x-xstart)^2+(y-ystart)^2)
        end
        i=1
        x_n=x
        y_n=y
        while offSetDistance+transitionMagnitude<1.0*DEFAULT_LANE_WIDTH && i<=20
            v_desire_n=v_desire+model.acc_desire*(Δt*2)
            v_desire_n=clamp(v_desire_n,model.v_min,model.v_max)
            Δs=(v_desire+v_desire_n)*(Δt*2)/2
            v_desire=v_desire_n

            if i==6
                model.v_desire=v_desire
            end

            offSetDistance=offSetDistance+transitionMagnitude
            x_n,y_n=move_along_cure(fit,x_n,y_n,Δs,direction,mode)
            curvePtVecPer=get_curvePt_vec_per(fit,x_n,y_n,turnDirection,direction,mode)
            pos=[x_n;y_n]+offSetDistance*curvePtVecPer
            model.Pts[1,i]=pos[1]
            model.Pts[2,i]=pos[2]
            i=i+1
        end
        while offSetDistance+transitionMagnitude>=1.0*DEFAULT_LANE_WIDTH && i<=20
            v_desire_n=v_desire+model.acc_desire*(Δt*2)
            v_desire_n=clamp(v_desire_n,model.v_min,model.v_max)
            Δs=(v_desire+v_desire_n)*(Δt*2)/2
            v_desire=v_desire_n

            if i==6
                model.v_desire=v_desire
            end

            x_n,y_n=move_along_curve(fit,x_n,y_n,Δs,direction,mode)
            curvePtVecPer=get_curvePt_vec_per(fit,x_n,y_n,turnDirection,direction,mode)
            pos=[x_n;y_n]+1.0*DEFAULT_LANE_WIDTH*curvePtVecPer
            model.Pts[1,i]=pos[1]
            model.Pts[2,i]=pos[2]
            i=i+1
        end
    end

    #println("x waypoint : ",model.Pts[1,:])
    #println("y waypoint : ",model.Pts[2,:])
    model.index=1
    model.subindex=0
    model
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
using Wanderer
using Test

DataSource1 = Dict(
    :a => [1, 2, 3],
    :b => [2, 3, 4],
    :c => [3, 4, 5]
)

DataSource2 = Dict(
    :a => ["a", "b", "c"],
    :b => [2, 3, 4],
    :c => [3, 4, 5]
)

@testset "Wanderer.jl" begin
    @test (DataSource1 |> @select _.a + _.b + _.c => sum)[:sum] == [6, 9, 12]
    @test (DataSource1 |> @where _.a != 2, @select _.b - _.c => sub)[:sub] == [-1, -1]
    df = (DataSource1 |> @select _.(!1), @where _.c === 4)
    @test df[:b][1] == 3 && df[:c][1] == 4
    @test ["aa", "bbb", "cccc"] == (DataSource2 |> @select repeat(_.a, _.b) => a)[:a]
    f(x) = 2x + 1

    df = DataSource2 |>
    @select f(_.b) => b, _.(!b),
    @where _.c > 2, _.b > 2,
    @orderby _.b % 3

    @test df[:b] == [9, 7, 5] && df[:a] == ["c", "b", "a"] && df[:c] == [5, 4, 3]

    agg = DataSource1 |>
    @groupby (_.a + _.b - _.c) % 2 == 0 => is_even,
    @having is_even

    @test agg[:a][1] == [1, 3] && agg[:b][1] == [2, 4] && agg[:c][1] == [3, 5]

    df = DataSource1 |>
    @select _, _.a + _.b => "a + b",
    @limit 2

    @test df[Symbol("a + b")] == [3, 5]
end


text:
"""
Plot the trace of an SCF, i.e. the absolute error of the total energy at
each iteration versus the converged energy in a semilog plot. By default
a new plot canvas is generated, but an existing one can be passed and reused
along with `kwargs` for the call to `plot!`.
"""
function ScfPlotTrace(plt=plot(yaxis=:log); kwargs...)
    energies = Float64[]
    function callback(info)
        if info.stage == :finalize
            minenergy = minimum(energies[max(1, end-5):end])
            error = abs.(energies .- minenergy)
            error[error .== 0] .= NaN
            extra = ifelse(:mark in keys(kwargs), (), (mark=:x, ))
            plot!(plt, error; extra..., kwargs...)
            display(plt)
        else
            push!(energies, info.energies.total)
        end
    end
end

"""
Default callback function for `self_consistent_field`, which prints a convergence table
"""
function ScfDefaultCallback()
    prev_energies = nothing
    function callback(info)
        info.stage == :finalize && return
        if info.n_iter == 1
            E_label = haskey(info.energies, "Entropy") ? "Free energy" : "Energy"
            @printf "n     %-12s      Eₙ-Eₙ₋₁     ρout-ρin   Diag\n" E_label
            @printf "---   ---------------   ---------   --------   ----\n"
        end
        E = isnothing(info.energies) ? Inf : info.energies.total
        Estr  = (@sprintf "%+15.12f" round(E, sigdigits=13))[1:15]
        prev_E = prev_energies === nothing ? Inf : prev_energies.total
        Δρ = norm(info.ρout.fourier - info.ρin.fourier)
        ΔE = prev_E == Inf ? "      NaN" : @sprintf "% 3.2e" E - prev_E
        diagiter = sum(info.diagonalization.iterations) / length(info.diagonalization.iterations)
        @printf "% 3d   %s   %s   %2.2e   % 3.1f \n" info.n_iter Estr ΔE Δρ diagiter
        prev_energies = info.energies
    end
    callback
end

"""
Flag convergence as soon as total energy change drops below tolerance
"""
function ScfConvergenceEnergy(tolerance)
    energy_total = NaN

    function is_converged(info)
        info.energies === nothing && return false # first iteration

        # The ρ change should also be small, otherwise we converge if the SCF is just stuck
        norm(info.ρout.fourier - info.ρin.fourier) > 10sqrt(tolerance) && return false

        etot_old = energy_total
        energy_total = info.energies.total
        abs(energy_total - etot_old) < tolerance
    end
    return is_converged
end

"""
Flag convergence by using the L2Norm of the change between
input density and unpreconditioned output density (ρout)
"""
function ScfConvergenceDensity(tolerance)
    info -> norm(info.ρout.fourier - info.ρin.fourier) < tolerance
end

"""
Determine the tolerance used for the next diagonalization. This function takes
``|ρnext - ρin|`` and multiplies it with `ratio_ρdiff` to get the next `diagtol`,
ensuring additionally that the returned value is between `diagtol_min` and `diagtol_max`
and never increases.
"""
function ScfDiagtol(;ratio_ρdiff=0.2, diagtol_min=nothing, diagtol_max=0.03)
    function determine_diagtol(info)
        isnothing(diagtol_min) && (diagtol_min = 100eps(real(eltype(info.ρin))))
        info.n_iter ≤ 0 && return diagtol_max
        info.n_iter == 1 && (diagtol_max /= 5)  # Enforce more accurate Bloch wave

        diagtol = norm(info.ρnext.fourier - info.ρin.fourier) * ratio_ρdiff
        diagtol = min(diagtol_max, diagtol)  # Don't overshoot
        diagtol = max(diagtol_min, diagtol)  # Don't undershoot
        @assert isfinite(diagtol)

        # Adjust maximum to ensure diagtol may only shrink during an SCF
        diagtol_max = min(diagtol, diagtol_max)
        diagtol
    end
end


text:
# Conservative Estimator
type ConservativeEstimator <: Pi0Estimator
end

const conservative = ConservativeEstimator()

function estimatepi0(pvalues::Array{Float64,1}, pi0estimator::ConservativeEstimator)
  one(Float64)
end


text:
using Dates

export
    substitute_latex

struct Equation
	text::AbstractString
	scale::Float64
	type::String # Either `display` or `inline`.
	extra_packages::String # For example, \usepackage{tikz}.
    filename::String # Output filename without extension.
end

struct DisplayImage
	eq::Equation
	im_dir::AbstractString
	im_name::AbstractString
end

struct InlineImage
	eq::Equation
	im_dir::AbstractString
	im_name::AbstractString
	# Height could allow for scaling depth if necessary.
	myheight::Float64 # Height according to `.sizes` file.
	mydepth::Float64 # Depth according to `.sizes` file.
end

cache_path(im_dir) = joinpath(im_dir, "cache.txt")

# Source: https://stackoverflow.com/questions/14182879.
latex_regex = r"(?<!\\)    # negative look-behind to make sure start is not escaped
    (?:        # start non-capture group for all possible match starts
      # group 1, match dollar signs only
      # single or double dollar sign enforced by look-arounds
      ((?<!\$)\${1,2}(?!\$))|
      # group 2, match escaped parenthesis
      (\\\()|
      # group 3, match escaped bracket
      (\\\[)|
      # group 4, match begin equation
      (\\begin\{equation\})
    )
    # if group 1 was start
    (?(1)
      # non greedy match everything in between
      # group 1 matches do not support recursion
      (.*?)(?<!\\)
      # match ending double or single dollar signs
      (?<!\$)\1(?!\$)|
    # else
    (?:
      # greedily and recursively match everything in between
      # groups 2, 3 and 4 support recursion
      (.*(?R)?.*)(?<!\\)
      (?:
        # if group 2 was start, escaped parenthesis is end
        (?(2)\\\)|
        # if group 3 was start, escaped bracket is end
        (?(3)\\\]|
        # else group 4 was start, match end equation
        \\end\{equation\}
      )
    ))))|
    # Match multiline double dollar symbols.
    \${2}(?!\$)[\s\S]*?\${2}|
    # Match multiline backslash and square bracket
    \\\[(?!\$)[\s\S]*?\\\]
    "x

floatregex = "([0-9]*[.])?[0-9]+"
match2num(m::RegexMatch) = parse(Float64, match(Regex(floatregex), m.match).match)

preamble(extra_packages="") = """
    \\documentclass[12pt]{article}
    \\nonstopmode
    \\usepackage{amsmath}
    \\usepackage{amssymb}
    $(extra_packages)
    \\pagestyle{empty}"""

function wrap_eq(eq::Equation)::String
    if eq.type == "display"
        # The `lrbox` errors on display math (`$$ ... $$`).
        return string(
            preamble(eq.extra_packages), """
            \\begin{document}
            $(eq.text)
            \\end{document}"""
        )
    else
        # Using the `lrbox` to determine the baseline and resulting depth.
        # Source: http://mactextoolbox.sourceforge.net/articles/baseline.html
        return string(
            preamble(eq.extra_packages), raw"""
            \newsavebox{\mybox}
            \newlength{\mywidth}
            \newlength{\myheight}
            \newlength{\mydepth}
            \begin{lrbox}{\mybox}
            """,
            eq.text, '\n',
            raw"""
            \end{lrbox}
            \settowidth {\mywidth}  {\usebox{\mybox}}
            \settoheight{\myheight} {\usebox{\mybox}}
            \settodepth {\mydepth}  {\usebox{\mybox}}
            \newwrite\foo
            \immediate\openout\foo=\jobname.sizes
                \immediate\write\foo{Depth = \the\mydepth}
                \immediate\write\foo{Height = \the\myheight}
                \addtolength{\myheight} {\mydepth}
                \immediate\write\foo{TotalHeight = \the\myheight}
                \immediate\write\foo{Width = \the\mywidth}
            \closeout\foo
            \begin{document}
            \usebox{\mybox}
            \end{document}"""
        )
    end
end

function check_latex()
    try
        run(pipeline(`pdflatex --help`, devnull))
    catch
        throw(ErrorException("Failed to run pdflatex"))
    end
end

function get_sizes(path::AbstractString)::Tuple{Float64,Float64}
    io = open(path, "r")
    sizes = read(io, String)
    close(io)
    height = match(Regex("TotalHeight = $(floatregex)pt"), sizes)
    depth = match(Regex("Depth = $(floatregex)pt"), sizes)
    (match2num(height), match2num(depth))
end

function latex_im(eq::Equation, im_dir::String)::Union{DisplayImage,InlineImage}
    check_latex()
    tmpdir = tempname() * '/'
    mkdir(tmpdir)
    file(extension) = joinpath(tmpdir * "eq.$extension")
    open(file("tex"), "w") do io
        write(io, wrap_eq(eq))
    end
    old_pwd = pwd()
    # pdflatex uses current working directory to store intermediate files.
    cd(tmpdir)
    pdflatex = `pdflatex $(file("tex"))`
    mktemp() do path, file
        try
            run(pipeline(pdflatex, stdout=path, stderr=devnull))
        catch e
            open(path, "r") do io
                println(read(path, String))
            end
            throw(ErrorException("""Failed to run pdflatex on
                $(wrap_eq(eq))
                """))
        end
    end
    crop = `pdfcrop $(file("pdf")) $(file("crop.pdf"))`
    try
        run(pipeline(crop, stdout=devnull, stderr=devnull))
    catch e
        throw(ErrorException("Failed to run pdfcrop. Is pdfcrop installed?"))
    end
    pdf2svg = `dvisvgm --pdf $(file("crop.pdf"))`
    try
        run(pipeline(pdf2svg, stdout=devnull, stderr=devnull))
    catch e
        throw(ErrorException("Failed to run dvisvgm. Is dvisvgm installed?"))
    end
    mv(file("crop"), file("svg"))
    tmpfilename = split(tmpdir, '/')[end-1]
    im_name = eq.filename * ".svg"
    im_path = joinpath(im_dir, im_name)
    mv(file("svg"), im_path, force=true)
    cd(old_pwd)
    if eq.type == "display"
        eq_image = DisplayImage(eq, im_dir, im_name)
    else
        (height, depth) = get_sizes(file("sizes"))
        eq_image =  InlineImage(eq, im_dir, im_name, height, depth)
    end
    rm(tmpdir, recursive=true)
    return eq_image
end

function dimensions(svg_path::AbstractString)::Tuple{Float64,Float64}
    svg = ""
    try
        io = open(svg_path, "r")
        svg = read(io, String)
        close(io)
    catch
        throw(ErrorException("Failed to read file. Is the cache invalid?"))
    end
    w = match(Regex("width='$(floatregex)pt'"), svg)
    h = match(Regex("height='$(floatregex)pt'"), svg)
    return (match2num(w), match2num(h))
end

function determine_param(eq_image)::Array{String,1}
    link = '/' * joinpath("latex", eq_image.im_name)
    im_path = joinpath(eq_image.im_dir, eq_image.im_name)
    (w, h) = dimensions(im_path)
    width = round(eq_image.eq.scale * w; digits=3)
    height = round(eq_image.eq.scale * h; digits=3)
    params = [
        "src=\"$(link)\"",
        "width=\"$(width)\"",
        "height=\"$(height)\""
    ]
    if eq_image.eq.type == "inline"
        valign = round(eq_image.eq.scale * (eq_image.mydepth); digits=3)
        style = "style=\"margin:0;vertical-align:-$(valign)px\""
        push!(params, style)
    end
    return params
end

function _eq!(cache, eq::Equation, im_dir)::String
    eq_image = check_cache(cache, eq)
    if eq_image == nothing
        eq_image = latex_im(eq, im_dir)
        update_cache!(cache, eq_image)
    end
    param = determine_param(eq_image)
    img = """<img $(join(param, ' '))>"""
    s = eq.type == "display" ? "<center>$(img)</center>" : img
end

function eq_replace!(cache, eq::SubString, scale, im_dir, extra_packages)::String
    is_display = (startswith(eq, raw"$$") || startswith(eq, raw"\[") || startswith(eq, raw"\begin{eq"))
    eq_type = is_display ? "display" : "inline"
    _eq!(cache, Equation(eq, scale, eq_type, extra_packages, string(hash(eq))), im_dir)
end

"""
    substitute_latex(md::AbstractString, scale::Number, im_dir; extra_packages="")::String

Substitute LaTeX in Markdown string `md`.
The LaTeX images will be placed at `im_dir/<h>.svg` where `h` is a hash calculated over the math expression.
Tweaking the image size is possible by setting `scale`.

The benefit of using a hash is that the browser downloads only one image per math expression.
"""
function substitute_latex(md::AbstractString, scale::Number, im_dir; extra_packages="")::String
	if !(isdir(im_dir)); mkdir(im_dir) end
	cache = load_cache(scale, im_dir)
	initial_cache_length = length(cache.images)

    md = replace(md, latex_regex => m -> eq_replace!(cache, m, scale, im_dir, extra_packages))

	if length(cache.images) != initial_cache_length
		write_cache!(cache, im_dir)
	end
	md
end


text:
function RoboDojo.indices_z(model::DoublePendulum{T,Nominal}) where T
    nq = model.nq
    nc = model.nc
    q = collect(1:nq)
    γ = collect(1:0)
    sγ = collect(1:0)
    ψ = collect(1:0)
    b = collect(1:0)
    sψ = collect(1:0)
    sb = collect(1:0)
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

RoboDojo.num_var(model::DoublePendulum{T,Nominal}) where T = model.nq

function RoboDojo.indices_optimization(model::DoublePendulum{T,Nominal}) where T
    nq = model.nq
    nz = nq
    IndicesOptimization(
        nz,
        nz,
        [collect(1:0), collect(1:0)],
        [collect(1:0), collect(1:0)],
        Vector{Vector{Vector{Int}}}(),
        Vector{Vector{Vector{Int}}}(),
        collect(1:nq),
        collect(1:0),
        collect(1:0),
        Vector{Vector{Int}}(),
        collect(1:0))
end

function RoboDojo.initialize_z!(z, model::DoublePendulum{T,Nominal}, idx::IndicesZ, q) where T
    z[idx.q] .= q
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9c2efd8a-26b6-11ec-2927-f55ab3486629
begin
	using ConstraintSolver
	using JuMP
	using IterTools: iterated
	using Base.Iterators: take
	using CommonMark
	using PlutoUI
end

# ╔═╡ 0bb13e45-85da-45e4-a17a-bc56baccf28d
const CS = ConstraintSolver

# ╔═╡ fc760fa6-b9b3-4eac-bf7c-4ad3c356caaa
@bind description TextField(default="SEND + MORE = MONEY")

# ╔═╡ 3c149796-1639-4af1-b2b9-f49078db9906
begin
	m = Model(CS.Optimizer)
	parts = strip.(split(description, "="))
	str_result = parts[2]
	str_sums = strip.(split(parts[1], "+"))
	str_sums, str_result
end

# ╔═╡ 7d409d6c-b9cc-4ea4-9932-6d003e0abe0b

# ╔═╡ 25c5e105-ded0-40ca-89e0-85e70eb32a35
var_names = unique(vcat(collect(str_result), collect.(str_sums)...))

# ╔═╡ 4e211e3d-cec8-48ed-a5b1-bc198a8c0fc2
@variable(m, 0 <= vars[Symbol.(var_names)] <= 9, Int) |> Text

# ╔═╡ 8f2623ed-e4b4-4501-8cc4-54b76fe39496
begin
	start_vars = unique([[s[1] for s = str_sums]..., str_result[1]])
	for s in Symbol.(start_vars)
		@constraint(m, vars[s] >= 1)
	end
end

# ╔═╡ ab374d36-79ed-44dd-bb84-2c27bbe5a893
@constraint(m, vars in CS.AllDifferent())

# ╔═╡ f15f1a14-8495-45fb-bb67-dbc68bbbfe1a
function word_to_sum(word, vars, m)
	word_vars = [Symbol(c) for c=word]
	times = reverse(collect(take(iterated(x -> 10x, 1), length(word))))

	mults = map(zip(word_vars, times)) do (w, t)
		vars[w]*t
	end
	sum(mults)
end

# ╔═╡ 4f5498f2-70d4-4283-a307-7cb551b54abb
left = [word_to_sum(str, vars, m) for str in str_sums]

# ╔═╡ 807fc2d8-f472-4920-8380-a6abb1a654e6
right = word_to_sum(str_result, vars, m)

# ╔═╡ d36a8f20-97ee-4828-b09b-ee320fc11bee
@constraint(m, sum(left) == sum(right))

# ╔═╡ cfa9df6f-81d5-40a5-b029-4a83311e74d4
begin
	optimize!(m)
	var_values = [round(Int, value(vars[Symbol(name)])) for name = var_names]
	var_names_s = string.(var_names)
	replacements = Pair.(var_names, var_values)

	let result = description
	for sub = replacements
		result = replace(result, sub)
	end

	HTML("""<pre>
			$description
			$result
			</pre>
			""")
	end
end

# ╔═╡ 881128de-513c-4588-b691-1efc836e6c96
status = JuMP.termination_status(m)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
ConstraintSolver = "e0e52ebd-5523-408d-9ca3-7641f1cd1405"
IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
CommonMark = "~0.8.3"
ConstraintSolver = "~0.7.0"
IterTools = "~1.3.0"
JuMP = "~0.21.10"
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[Arpack]]
deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
version = "0.4.0"

[[Arpack_jll]]
deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.0+3"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.0"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.7.2"

[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[ConstraintProgrammingExtensions]]
deps = ["LinearAlgebra", "MathOptInterface", "Test"]
git-tree-sha1 = "989acd05e57112a5fd1bc8849ebc398a614aa010"
uuid = "b65d079e-ed98-51d9-b0db-edee61a5c5f8"
version = "0.3.2"

[[ConstraintSolver]]
deps = ["ConstraintProgrammingExtensions", "DataStructures", "Formatting", "JSON", "JuMP", "LightGraphs", "MathOptInterface", "MatrixNetworks", "Random", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "47a4c07efa6890e7da286495e22742eb19f3eed9"
uuid = "e0e52ebd-5523-408d-9ca3-7641f1cd1405"
version = "0.7.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "7220bc21c33e990c14f4a9a319b1d242ebc5b269"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.1"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "c4203b60d37059462af370c4f3108fb5d155ff13"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.20"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSONSchema]]
deps = ["HTTP", "JSON", "URIs"]
git-tree-sha1 = "2f49f7f86762a0fbbeef84912265a1ae61c4ef80"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.4"

[[JuMP]]
deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Printf", "Random", "SparseArrays", "SpecialFunctions", "Statistics"]
git-tree-sha1 = "4358b7cbf2db36596bdbbe3becc6b9d87e4eb8f5"
uuid = "4076af6c-e467-56ae-b986-b466b2749572"
version = "0.21.10"

[[KahanSummation]]
deps = ["Test"]
git-tree-sha1 = "1f01068b28d3ad83d4d1212a0ce8d7ecacb33482"
uuid = "8e2b3108-d4c1-50be-a7a2-16352aec75c3"
version = "0.1.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "575644e3c05b258250bb599e57cf73bbf1062901"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.9.22"

[[MatrixNetworks]]
deps = ["Arpack", "DataStructures", "DelimitedFiles", "IterTools", "KahanSummation", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "f6794d987d426bab570aabb2d03e2f40ca4fc438"
uuid = "4f449596-a032-5618-b826-5a251cb6dc11"
version = "1.0.2"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "95072ef1a22b057b1e80f73c2a89ad238ae4cfff"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.12"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═9c2efd8a-26b6-11ec-2927-f55ab3486629
# ╠═0bb13e45-85da-45e4-a17a-bc56baccf28d
# ╠═fc760fa6-b9b3-4eac-bf7c-4ad3c356caaa
# ╠═3c149796-1639-4af1-b2b9-f49078db9906
# ╠═7d409d6c-b9cc-4ea4-9932-6d003e0abe0b
# ╠═25c5e105-ded0-40ca-89e0-85e70eb32a35
# ╠═4e211e3d-cec8-48ed-a5b1-bc198a8c0fc2
# ╠═8f2623ed-e4b4-4501-8cc4-54b76fe39496
# ╠═ab374d36-79ed-44dd-bb84-2c27bbe5a893
# ╠═f15f1a14-8495-45fb-bb67-dbc68bbbfe1a
# ╠═4f5498f2-70d4-4283-a307-7cb551b54abb
# ╠═807fc2d8-f472-4920-8380-a6abb1a654e6
# ╠═d36a8f20-97ee-4828-b09b-ee320fc11bee
# ╠═cfa9df6f-81d5-40a5-b029-4a83311e74d4
# ╠═881128de-513c-4588-b691-1efc836e6c96
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
struct GradientLocalOptimizer <: LocalOptimizer
    algo::Symbol
    options::Dict{Symbol,Any}
end

need_gradient(::GradientLocalOptimizer) = true

function LBFGSOptimizer(; options...)
    GradientLocalOptimizer(:LD_LBFGS, Dict{Symbol,Any}(options))
end

function minimize(o::GradientLocalOptimizer, f, lb, ub, x0)
    opt = create_nlopt(o.algo, lb, ub; min_objective = f, o.options...)
    fx, x, ret = NLopt.optimize(opt, x0)
    ret == :FORCED_STOP && throw(InterruptException())
    return fx, x, ret
end

function maximize(o::GradientLocalOptimizer, f, lb, ub, x0)
    opt = create_nlopt(o.algo, lb, ub; max_objective = f, o.options...)
    fx, x, ret = NLopt.optimize(opt, x0)
    ret == :FORCED_STOP && throw(InterruptException())
    return fx, x, ret
end


text:
abstract type AbstractStatus end
struct Optimal <: AbstractStatus end
struct Infeasible <: AbstractStatus end
struct Unknown <: AbstractStatus end
struct Unbounded <: AbstractStatus end


text:
# poly kernel: (s'*x+c)^d

type Poly <: Kernel; n; s; w; v; x; y; dy; c; d;
# TODO: remove PercLoss, integrate into kernel
# TODO: use daume trick for kernel perceptron
# TODO: define all in kperceptron (mirroring perceptron.jl) and have kernel fn as a field
    function Poly(nclass;o...)
        l = setparam!(new(nclass); o...)
        isdefined(l,:c) || (l.c = 0)
        isdefined(l,:d) || (l.d = 1)
        return l
    end
end

# copy(l::Poly; o...)=Poly(s=copy(l.s), w=copy(l.w), v=copy(l.v), c=l.c, d=l.d)
setparam!(l::Poly; o...)=(for (n,v) in o; l.(n)=v; end; l)
kernel(l::Poly)=((l.s' * l.x + l.c) .^ l.d)