text:
######################################################################
# Unit Tests for Lockables submodule
# -----
# Licensed under Apache License 2.0
module TestLockables
using Test
using Pia.Lockables

macro test_return(cond)
    esc(:(if !$cond return false end))
end

@testset "Lockables" begin
    @testset "lock" begin
        let lockable = Lockable(42)
            @test !islocked(lockable)
            lock(lockable) do
                @test islocked(lockable)
            end
            @test !islocked(lockable)
        end
    end

    @testset "sync" begin
        lockable = Lockable(42)
        t0 = t1 = t2 = time()
        success = false
        dt1 = 0.1
        dt2 = 0.2

        @sync begin
            @async begin
                lock(lockable) do
                    sleep(dt2)
                    lockable[] = 69
                end
            end

            @async begin
                sleep(dt1)
                t0 = time()
                lock(lockable) do
                    t1 = time()
                    success = lockable[] == 69
                end
                t2 = time()
            end
        end

        @test success
        @test t2-t0 < dt2
        @test isapprox(t1, t2)
    end
end

end # module TestLockables


text:
abstract type Execution end
"""
    Serial

Factory object for [`LShapedSolvers.SerialExecution`](@ref)/[`ProgressiveHedgingSolvers.SerialExecution`](@ref). Passed by default to `execution ` in the `LShapedSolver` and `ProgressiveHedgingSolver` factory functions.

"""
struct Serial <: Execution end

function str(::Serial)
    return ""
end

"""
    Synchronous

Factory object for [`LShapedSolvers.SynchronousExecution`](@ref)/[`ProgressiveHedgingSolvers.SynchronousExecution`](@ref). Pass to `execution` in the `LShapedSolver` and `ProgressiveHedgingSolver` factory functions.

"""
struct Synchronous <: Execution end

function str(::Synchronous)
    return "Synchronous "
end

"""
    Asynchronous

Factory object for [`LShapedSolvers.AsynchronousExecution`](@ref)/[`ProgressiveHedgingSolvers.AsynchronousExecution`](@ref). Pass to `execution` in the `LShapedSolver` and `ProgressiveHedgingSolver` factory functions

...
# Parameters
- `max_active::Int = 3`: Maximum number of active iterations that run asynchronously.
- `κ::T = 0.5`: Relative amount of finished subproblems required to start a new iterate. Governs the amount of asynchronicity.
...
"""
struct Asynchronous <: Execution
    max_active::Int
    κ::Float64

    function Asynchronous(; max_active::Int = 3, κ::AbstractFloat = 0.5)
        return new(max_active, κ)
    end
end

function str(::Asynchronous)
    return "Asynchronous "
end

ext: jl


text:
"""
ctmcst
"""

# struct CTMCModel{MatrixT,Tv}
#     Q::MatrixT{Tv}
#     initv::Vector{Tv}
#     reward::Dict{Symbol,Vector{Tv}}
#     states::Vector{Symbol}
#     index::Dict{Symbol,Ti}
# end

# function exrss(m::CTMCModel{MatrixT,Tv}; reward = :default) where {MatrixT,Tv}
# end

"""
macro
"""

mutable struct Markov
    Tv::DataType
    state::Set{Symbol}
    initial::Dict{Symbol,Any}
    trans::Dict{Tuple{Symbol,Symbol},Any}
    reward::Dict{Symbol,Dict{Symbol,Any}}
    ireward::Dict{Symbol,Dict{Tuple{Symbol,Symbol},Any}}
end

function Base.show(io::IO, x::Markov)
    Base.show(io, "Markov $(objectid(x))")
end

function Markov()
    Markov(Float64, Set{Symbol}(), Dict{Symbol,Any}(), Dict{Tuple{Symbol,Symbol},Any}(),
        Dict{Symbol,Dict{Symbol,Any}}(), Dict{Symbol,Dict{Tuple{Symbol,Symbol},Any}}())
end

function trans!(m::Markov, src::Symbol, dest::Symbol, tr)
    push!(m.state, src)
    push!(m.state, dest)
    m.Tv = promote_type(m.Tv, typeof(tr))
    m.trans[(src,dest)] = tr
end

function reward!(m::Markov, label::Symbol, s::Symbol, r)
    push!(m.state, s)
    d = get(m.reward, label) do
        m.reward[label] = Dict{Symbol,Any}()
    end
    m.Tv = promote_type(m.Tv, typeof(r))
    d[s] = r
end

function reward!(m::Markov, label, src, dest, r)
    push!(m.state, src)
    push!(m.state, dest)
    d = get(m.ireward, label) do
        m.ireward[label] = Dict{Tuple{Symbol,Symbol},Any}()
    end
    m.Tv = promote_type(m.Tv, typeof(r))
    d[(src,dest)] = r
end

function initial!(m::Markov, s::Symbol, p)
    push!(m.state, s)
    m.Tv = promote_type(m.Tv, typeof(p))
    m.initial[s] = p
end

function generate(m; modeltype::Symbol = :SparseCTMC)
    _generate(Val(modeltype), m)
end

function _generate(::Val{:SparseCTMC}, m::Markov)
    states = [x for x = m.state]
    index = Dict([states[i] => i for i = 1:length(states)]...)
    Q = spzeros(m.Tv, length(states), length(states))
    for ((src,dest),t) = m.trans
        Q[index[src],index[dest]] = t
        Q[index[src],index[src]] -= t
    end
    initv = zeros(m.Tv, length(states))
    for (s,p) = m.initial
        initv[index[s]] = p
    end
    rwd = Dict{Symbol,Vector{m.Tv}}()
    for (k,v) = m.reward
        rwdv = zeros(m.Tv, length(states))
        for (s,r) = v
            rwdv[index[s]] = r
        end
        rwd[k] = rwdv
    end
    Q, initv, rwd, states
end

function _generate(::Val{:DenseCTMC}, m::Markov)
    states = [x for x = m.state]
    index = Dict([states[i] => i for i = 1:length(states)]...)
    Q = zeros(m.Tv, length(states), length(states))
    for ((src,dest),t) = m.trans
        Q[index[src],index[dest]] = t
        Q[index[src],index[src]] -= t
    end
    initv = zeros(m.Tv, length(states))
    for (s,p) = m.initial
        initv[index[s]] = p
    end
    rwd = Dict{Symbol,Vector{m.Tv}}()
    for (k,v) = m.reward
        rwdv = zeros(m.Tv, length(states))
        for (s,r) = v
            rwdv[index[s]] = r
        end
        rwd[k] = rwdv
    end
    Q, initv, rwd, states
end

macro tr(m, block)
    if Meta.isexpr(block, :block)
        body = [_gentrans(x, m) for x = block.args]
        esc(Expr(:block, body...))
    else
        esc(_gentrans(block, m))
    end
end

function _gentrans(x::Any, m)
    x
end

function _gentrans(x::Expr, m)
    if Meta.isexpr(x, :tuple) && Meta.isexpr(x.args[1], :call) && x.args[1].args[1] == :(=>) && length(x.args) == 2
        src = x.args[1].args[2]
        dest = x.args[1].args[3]
        t = x.args[2]
        :(trans!($m, $src, $dest, $t))
    else
        throw(TypeError(x, "Invalid format for the transition"))
    end
end

macro init(m, block)
    if Meta.isexpr(block, :block)
        body = [_geninitial(x, m) for x = block.args]
        esc(Expr(:block, body...))
    else
        esc(_geninitial(block, m))
    end
end

function _geninitial(x::Any, m)
    x
end

function _geninitial(x::Expr, m)
    if Meta.isexpr(x, :tuple) && length(x.args) == 2
        s = x.args[1]
        p = x.args[2]
        :(initial!($m, $s, $p))
    else
        throw(TypeError(x, "Invalid format for the initial probability"))
    end
end

macro reward(m, label, block)
    if Meta.isexpr(block, :block)
        body = [_genreward(x, label, m) for x = block.args]
        esc(Expr(:block, body...))
    else
        esc(_genreward(block, label, m))
    end
end

function _genreward(x::Any, label, m)
    x
end

function _genreward(x::Expr, label, m)
    if Meta.isexpr(x, :tuple) && length(x.args) == 2
        s = x.args[1]
        r = x.args[2]
        :(reward!($m, $label, $s, $r))
    elseif Meta.isexpr(x, :tuple) && Meta.isexpr(x.args[1], :call) && x.args[1].args[1] == :(=>) && length(x.args) == 2
        src = x.args[1].args[2]
        dest = x.args[1].args[3]
        r = x.args[2]
        :(reward!($m, $label, $src, $dest, $r))
    else
        throw(TypeError(x, "Invalid format for the reward vector"))
    end
end

###

macro markov(f, block)
    body = []
    push!(body, Expr(:(=), :tmp, Expr(:call, :Markov)))
    if Meta.isexpr(block, :block)
        for x = block.args
            push!(body, _replace_macro(x))
        end
    end
    push!(body, :tmp)
    esc(Expr(:function, f, Expr(:block, body...)))
end

function _replace_macro(x::Any)
    x
end

function _replace_macro(x::Expr)
    if Meta.isexpr(x, :macrocall) && (x.args[1] == Symbol("@tr") || x.args[1] == Symbol("@init") || x.args[1] == Symbol("@reward"))
        Expr(:macrocall, x.args[1], x.args[2], :tmp, [_replace_macro(u) for u = x.args[3:end]]...)
    else
        Expr(x.head, [_replace_macro(u) for u = x.args]...)
    end
end

macro s(x...)
    esc(Expr(:call, :Symbol, x...))
end


text:
include("show.jl")
include("normalderivative.jl")


text:
module BiComplex

include("implementation.jl")

export Bicomplex, derivative, jm, second

end


text:
#function DemoLSTMsimple(PlotResults=true)
# Simple memory problem:
# If either of the first two inputs in x[t] is 1, then output to y[t] the previous memory m[t-1] and update the memory m[t] to store x[t]
# This can be achieved using a simple LSTM style approach.

useproc("CPU")

PlotResults=true

T=400
m=Array(ADnode,T)
alpha=Array(ADnode,T)
beta=Array(ADnode,T)
x=Array(ADnode,T)
y=Array(ADnode,T)
yout=Array(ADnode,T)
loss=Array(ADnode,T)

StartCode()
Win=ADvariable()
Wout=ADvariable()
Bin=ADvariable()
Bout=ADvariable()
x[1]=ADnode()
alpha[1]=sigmoid(Win*x[1])
m[1]=1*x[1]
y[1]=ADnode()
beta[1]=sigmoid(Wout*x[1])
for t=2:T
    x[t]=ADnode()
    y[t]=ADnode()
    alpha[t]=sigmoid(Win*x[t]+Bin)
    m[t]=alpha[t]*x[t]+(1-alpha[t])*m[t-1]
    beta[t]=sigmoid(Wout*x[t]+Bout)
    yout[t]=beta[t]*m[t-1]
    #loss[t]=meanSquareLoss(yout[t],y[t])
    loss[t]=BinaryKullbackLeiblerLoss(y[t],yout[t])
end
totalloss=mean(loss[2:T])
net=EndCode() # defines the graph

# instantiate parameter nodes and inputs:
X=4
# parameters:
net.value[Win]=0.01*randn(1,X)
net.value[Wout]=0.01*randn(1,X)
net.value[Bin]=0.01*randn(1,1)
net.value[Bout]=0.01*randn(1,1)

# training data:
# make close to zero or 1 (strictly 0 or 1 can give NAN in KL loss)
net.value[x[1]]=zeros(X,1)+0.001; net.value[x[1]][1]=0.999;
net.value[y[1]]=zeros(X,1)+0.001

xstore=copy(net.value[x[1]])
for t=2:T
    net.value[x[t]]=0.001+0.998*(rand(X,1).>0.5)
    net.value[x[t]][1]=0.001+0.998*(rand()>0.95)

    if sum(extract(net.value[x[t]][1]))>0.5
        net.value[y[t]]=xstore
        xstore=net.value[x[t]]
    else
        net.value[y[t]]=zeros(X,1)+0.001
    end
end

net=compile(net) # compile and preallocate memory
@gpu CUDArt.init([0])
@gpu net=convert(net,"GPU")

#gradcheck(net;showgrad=true) # use a small number of datapoints and small network to check the gradient, otherwise

# Training:
println("Training: using $(net.gpu==true? "GPU" : "CPU") ")
ParsToUpdate=Parameters(net)
nupdates=1000
error=Array(Float64,0)
    velo=NesterovInit(net) # Nesterov velocity
    for i=1:nupdates
    LearningRate=0.5
    ADforward!(net)
    ADbackward!(net)
    push!(error,extract(net.value[net.FunctionNode])[1])
    printover("iteration $i: training loss = $(error[i])")
    for par in ParsToUpdate
        #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
        NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,i)
    end
end
println()

if PlotResults==true
    yy=zeros(X,T); yyout=zeros(X,T); xx=zeros(X,T)
    for t=2:T
        yy[:,t]=net.value[y[t]]; yyout[:,t]=net.value[yout[t]]
        xx[:,t]=net.value[x[t]];
    end

    fig,ax=PyPlot.subplots(3,1,sharex=true)
    ax[1,1][:imshow](xx,interpolation="none");ax[1,1][:set_aspect](15)
    ax[2,1][:imshow](yy,interpolation="none"); ax[2,1][:set_aspect](15);
    ax[3,1][:imshow](yyout,interpolation="none"); ax[3,1][:set_aspect](15)
    ax[1,1][:set_title]("x: train input")
    ax[2,1][:set_title]("y: train output")
    ax[3,1][:set_title]("yout: net predicted output given x")
    fig[:canvas][:draw]()

    figure(); plot(error); title("training loss")
end
#end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
export Window,
       include_asset

immutable Window
    alive::Signal
    dimension::Signal
    route::Signal
    dir::AbstractString
    assets::Signal
end

Window(;
    alive=true,
    dimension=(0px, 0px),
    route="",
    dir="ltr") =
    Window(Signal(Bool, alive), Signal(Any, dimension), Signal(Any, route), "ltr", Signal(Any, "basics"))

resolve_asset(tup :: (@compat Tuple{AbstractString, AbstractString}), prefix ="pkg", joinfn=(x, y) -> x * "/" * y) = begin
    pkg = tup[1]
    slug = tup[2]
    path = Pkg.dir(pkg, "assets", slug * ".html")
    if isfile(path)
        return joinfn(prefix, joinfn( pkg, "$slug.html") )
    else
        error("Asset file $path doesn't exist")
    end
end

resolve_asset(slug::AbstractString) = resolve_asset(("Escher", slug))


text:
module ParallelUtilities

# package code goes here
include("macros.jl")

end # module


text:
"""
    couple_models(aero::LiftingLine, stru::GEBT, flap::LiftingLineFlaps)

Create a coupled model using a lifting line aerodynamic model, a
geometrically exact beam theory model, and a lifting line control surface model.
This model introduces additional parameters corresponding to the freestream
air density ``\\rho_\\infty``, followed by the external loads ``F_{x,i}, F_{y,i},
F_{z,i}, M_{x,i}, M_{y,i}, M_{z,i}`` or displacements ``u_{x,i}, u_{y,i}, u_{z,i},
\\theta_{x,i}, \\theta_{y,i}, \\theta_{z,i}`` for each node, followed by the constant
distributed loads ``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i}, m_{z,i}``
applied on each beam element (excluding aerodynamic loads), followed by the body
frame linear and angular velocities ``u, v, w, p, q, r`` and the control surface
deflections ``\\delta_1, \\delta_2, \\dots, \\delta_N``.

**NOTE: When using this model, the local frame for each beam element should be
oriented with the x-axis along the beam's axis, the y-axis forward, and the
z-axis normal to the surface**
"""
function couple_models(aero::LiftingLine, stru::GEBT, flap::LiftingLineFlaps)
    return (aero, stru, flap)
end

# --- traits --- #

function number_of_additional_parameters(aero::LiftingLine, stru::GEBT, flap::LiftingLineFlaps)

    return 1 + 6*length(stru.icol_point) + 6*length(stru.icol_elem) + 6 + length(flap.gains)
end

function coupling_inplaceness(::Type{<:LiftingLine}, ::Type{<:GEBT}, ::Type{<:LiftingLineFlaps})
    return InPlace()
end

function coupling_rate_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(islinear.(coupling_rate_jacobian_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

function coupling_state_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(islinear.(coupling_state_jacobian_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_parameter_jacobian_type(::Type{<:LiftingLine}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps}) = Nonlinear()

function coupling_time_gradient_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT},
    ::Type{LiftingLineFlaps{NF,NG,TF,TG}}) where {NA,NF,NG,TA,TF,TG}

    aero_model_types = TA.parameters
    flap_model_types = TF.parameters
    if all(isempty.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Empty()
    elseif all(iszero.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Zeros()
    elseif all(isconstant.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Constant()
    elseif all(islinear.(coupling_time_gradient_type.(aero_model_types, Ref(LiftingLineSection),
        flap_model_types, Ref(LiftingLineSectionControl))))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- methods --- #

function get_coupling_inputs!(y, aero::LiftingLine{NA,TA}, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, dx, x, p, t) where {NA,NF,NG,TA,TF,TG}

    # extract number of state variables, inputs, and parameters
    nua = number_of_states(aero) # number of aerodynamic states
    nus = number_of_states(stru) # number of structural states
    nuf = number_of_states(flap) # number of control surface states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nys = number_of_inputs(stru) # number of structural inputs
    nyf = number_of_inputs(flap) # number of control surface inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    nps = number_of_parameters(stru) # number of structural parameters
    npf = number_of_parameters(flap) # number of control surface parameters
    npadd = number_of_additional_parameters(aero, stru, flap) # number of additional parameters

    # get indices for state variables, inputs, and parameters
    iua = 1:nua # indices of aerodynamic states
    iya = 1:nya # indices of aerodynamic inputs
    ipa = 1:npa # indices of aerodynamic parameters
    ius = nua + 1 : nua + nus # indices of structural states
    iys = nya + 1 : nya + nys # indices of structural inputs
    ips = npa + 1 : npa + nps # indices of structural parameters
    iuf = nua + nus + 1 : nua + nus + nuf # indices of control surface states
    iyf = nya + nys + 1 : nya + nys + nyf # indices of control surface inputs
    ipf = npa + nps + 1 : npa + nps + npf # indices of control surface parameters
    ipadd = npa + nps + npf + 1 : npa + nps + npf + npadd # indices of additional parameters
    iuas = state_indices(aero.models) # indices of aerodynamic states for each section
    iyas = input_indices(aero.models) # indices of aerodynamic inputs for each section
    ipas = parameter_indices(aero.models) # indices of aerodynamic parameters for each section
    iufs = state_indices(flap.models) # indices of control surface states for each section
    iyfs = input_indices(flap.models) # indices of control surface inputs for each section
    ipfs = parameter_indices(flap.models) # indices of control surface parameters for each section

    # separate state variables, inputs, and parameters
    dua = view(dx, iua) # aerodynamic rates
    ua = view(x, iua) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dus = view(dx, ius) # structural rates
    us = view(x, ius) # structural states
    ys = view(y, iys) # structural inputs
    ps = view(p, ips) # structural parameters
    duf = view(dx, iuf) # control surface rates
    uf = view(x, iuf) # control surface states
    yf = view(y, iyf) # control surface inputs
    pf = view(p, ipf) # control surface parameters
    padd = view(p, ipadd) # additional parameters for coupled model
    duas = view.(Ref(dua), iuas) # aerodynamic rates for each section
    uas = view.(Ref(ua), iuas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section
    dufs = view.(Ref(duf), iufs) # control surface rates for each section
    ufs = view.(Ref(uf), iufs) # control surface states for each section
    yfs = view.(Ref(yf), iyfs) # control surface inputs for each section
    pfs = view.(Ref(pf), ipfs) # control surface parameters for each section

    # number of points and elements
    npoint = length(stru.icol_point)
    nelem = length(stru.icol_elem)

    # global parameters
    ρ = padd[1]
    ur = padd[1 + 6*npoint + 6*nelem + 1]
    vr = padd[1 + 6*npoint + 6*nelem + 2]
    wr = padd[1 + 6*npoint + 6*nelem + 3]
    pr = padd[1 + 6*npoint + 6*nelem + 4]
    qr = padd[1 + 6*npoint + 6*nelem + 5]
    rr = padd[1 + 6*npoint + 6*nelem + 6]
    δ = padd[SVector{NG}(1 + 6*npoint + 6*nelem + 6 + 1 : 1 + 6*npoint + 6*nelem + 6 + NG)]
    dδ = zero(δ)

    # body linear and angular velocity
    V = SVector(ur, vr, wr)
    Ω = SVector(pr, qr, rr)

    # construct assembly from structural parameters
    assembly = gxbeam_assembly(ps, npoint, nelem, stru.start, stru.stop)

    # save prescribed point loads/displacements
    for ip = 1:npoint
        yoff = 6*(ip-1)
        poff = 1 + 6*(ip-1)
        y[iys[yoff+1]] = padd[poff+1]
        y[iys[yoff+2]] = padd[poff+2]
        y[iys[yoff+3]] = padd[poff+3]
        y[iys[yoff+4]] = padd[poff+4]
        y[iys[yoff+5]] = padd[poff+5]
        y[iys[yoff+6]] = padd[poff+6]
    end

    # save aerodynamic inputs and distributed loads
    for i = 1:NA
        # models for this section
        section_aero = aero.models[i]
        section_stru = LiftingLineSection()
        section_flap = flap.models[i]
        section_ctrl = LiftingLineSectionControl()
        section_models = (section_aero, section_stru, section_flap, section_ctrl)

        # model dimensions for this section
        Nuai = number_of_states(section_aero)
        Nyai = number_of_inputs(section_aero)
        Npai = number_of_parameters(section_aero)
        Nusi = number_of_states(section_stru)
        Nysi = number_of_inputs(section_stru)
        Npsi = number_of_parameters(section_stru)
        Nufi = number_of_states(section_flap)
        Nyfi = number_of_inputs(section_flap)
        Npfi = number_of_parameters(section_flap)
        Nuci = number_of_states(section_ctrl)
        Nyci = number_of_inputs(section_ctrl)
        Npci = number_of_parameters(section_ctrl)

        # local section properties
        icol = stru.icol_elem[i]
        element = assembly.elements[i]
        u_elem = SVector(us[icol], us[icol+1], us[icol+2]) # linear displacement
        θ_elem = SVector(us[icol+3], us[icol+4], us[icol+5]) # angular displacement
        F_elem = SVector(us[icol+6], us[icol+7], us[icol+8]) .* stru.force_scaling # internal force
        M_elem = SVector(us[icol+9], us[icol+10], us[icol+11]) .* stru.force_scaling # internal moment
        P_elem = SVector(us[icol+12], us[icol+13], us[icol+14]) .* stru.mass_scaling # linear momentum
        H_elem = SVector(us[icol+15], us[icol+16], us[icol+17]) .* stru.mass_scaling # angular momentum

        dF_elem = SVector(dus[icol+6], dus[icol+7], dus[icol+8]) .* stru.force_scaling # internal force
        dM_elem = SVector(dus[icol+9], dus[icol+10], dus[icol+11]) .* stru.force_scaling # internal moment
        dP_elem = SVector(dus[icol+12], dus[icol+13], dus[icol+14]) .* stru.mass_scaling # linear momentum
        dH_elem = SVector(dus[icol+15], dus[icol+16], dus[icol+17]) .* stru.mass_scaling # angular momentum

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(θ_elem)
        θ_elem *= scaling # angular displacement (Wiener-Milenkovic parameters)

        # local to body transformation
        CtCab = GXBeam.get_C(θ_elem)'*element.Cab

        # local structural to aerodynamic transformation
        R = @SMatrix [0 -1 0; 1 0 0; 0 0 1]

        # local freestream linear and angular velocities
        vi = -R*(CtCab'*V + GXBeam.element_linear_velocity(element, P_elem, H_elem)) # local linear freestream velocity
        ωi = R*GXBeam.element_angular_velocity(element, P_elem, H_elem) # local angular freestream velocity

        # local freestream linear and angular accelerations
        dvi = -R*GXBeam.element_linear_velocity(element, dP_elem, dH_elem)
        dωi = R*GXBeam.element_angular_velocity(element, dP_elem, dH_elem)

        # NOTE: We assume that the local beam frame y-axis is oriented towards
        # the leading edge and the z-axis is oriented up
        # NOTE: We assume the local aerodynamic frame is oriented with the x-axis
        # in the chordwise direction and the z-axis in the (airfoil) normal
        # direction.

        # section state variable rates
        duai = SVector{Nuai}(duas[i]) # aerodynamic state variables
        dusi = vcat(dvi, dωi) # structural state variables
        dufi = SVector{Nufi}(dufs[i]) # control surface state variables
        duci = SMatrix{1,NG}(getindex.(flap.gains, i))* dδ # controller state variables
        dui = vcat(duai, dusi, dufi, duci)

        # section state variables
        uai = SVector{Nuai}(uas[i]) # aerodynamic state variables
        usi = vcat(vi, ωi) # structural state variables
        ufi = SVector{Nufi}(ufs[i]) # control surface state variables
        uci = SMatrix{1,NG}(getindex.(flap.gains, i)) * δ # controller state variables
        ui = vcat(uai, usi, ufi, uci)

        # section parameters
        pai = SVector{Npai}(pas[i]) # aerodynamic parameters
        psi = SVector(ρ) # structural parameters
        pfi = SVector{Npfi}(pfs[i]) # control surface parameters
        pci = SVector{0,Float64}() # controller parameters
        pi = vcat(pai, psi, pfi, pci)

        # section inputs
        yi = get_coupling_inputs(section_models, dui, ui, pi, t)

        # separate inputs
        yai = view(yi, 1:Nyai)
        ysi = view(yi, Nyai + 1 : Nyai + Nysi)
        yfi = view(yi, Nyai + Nysi + 1 : Nyai + Nysi + Nyfi)
        yci = view(yi, Nyai + Nysi + Nyfi + 1 : Nyai + Nysi + Nyfi + Nyci)

        # save local inputs
        y[iya[iyas[i]]] .= yai # aerodynamic inputs
        y[iyf[iyfs[i]]] .= yfi # control surface inputs

        # section aerodynamic loads (in body frame)
        fi = CtCab*R'*SVector(yi[Nyai+1], yi[Nyai+2], yi[Nyai+3])
        mi = CtCab*R'*SVector(yi[Nyai+4], yi[Nyai+5], yi[Nyai+6])

        # add constant distributed loads (in body frame)
        poff = 1 + 6*npoint + 6*(i-1)
        fi += SVector(padd[poff+1], padd[poff+2], padd[poff+3])
        mi += SVector(padd[poff+4], padd[poff+5], padd[poff+6])

        # save distributed loads for this element (in body frame)
        yoff = 6*npoint + 6*(i-1)
        y[iys[yoff+1 : yoff+3]] = fi
        y[iys[yoff+4 : yoff+6]] = mi
    end

    # save body frame linear/angular velocities
    yoff = 6*npoint + 6*nelem
    y[iys[yoff+1]] = ur
    y[iys[yoff+2]] = vr
    y[iys[yoff+3]] = wr
    y[iys[yoff+4]] = pr
    y[iys[yoff+5]] = qr
    y[iys[yoff+6]] = rr

    return y
end

# --- performance overloads --- #

# TODO

# --- convenience methods --- #

function set_additional_parameters!(padd, model::LiftingLine, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}; rho, point_conditions,
    element_loads, u, v, w, p, q, r, delta) where {NF,NG,TF,TG}

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    padd[1] = rho

    for ip = 1:np
        padd[1+6*(ip-1)+1] = point_conditions[6*(ip-1)+1]
        padd[1+6*(ip-1)+2] = point_conditions[6*(ip-1)+2]
        padd[1+6*(ip-1)+3] = point_conditions[6*(ip-1)+3]
        padd[1+6*(ip-1)+4] = point_conditions[6*(ip-1)+4]
        padd[1+6*(ip-1)+5] = point_conditions[6*(ip-1)+5]
        padd[1+6*(ip-1)+6] = point_conditions[6*(ip-1)+6]
    end

    for ie = 1:ne
        padd[1+6*np+6*(ie-1)+1] = element_loads[6*(ie-1)+1]
        padd[1+6*np+6*(ie-1)+2] = element_loads[6*(ie-1)+2]
        padd[1+6*np+6*(ie-1)+3] = element_loads[6*(ie-1)+3]
        padd[1+6*np+6*(ie-1)+4] = element_loads[6*(ie-1)+4]
        padd[1+6*np+6*(ie-1)+5] = element_loads[6*(ie-1)+5]
        padd[1+6*np+6*(ie-1)+6] = element_loads[6*(ie-1)+6]
    end

    padd[1 + 6*np + 6*ne + 1] = u
    padd[1 + 6*np + 6*ne + 2] = v
    padd[1 + 6*np + 6*ne + 3] = w
    padd[1 + 6*np + 6*ne + 4] = p
    padd[1 + 6*np + 6*ne + 5] = q
    padd[1 + 6*np + 6*ne + 6] = r
    padd[1 + 6*np + 6*ne + 6 + 1 : 1 + 6*np + 6*ne + 6 + NG] = delta

    return padd
end

function separate_additional_parameters(model::LiftingLine, stru::GEBT,
    flap::LiftingLineFlaps{NF,NG,TF,TG}, padd) where {NF,NG,TF,TG}

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    rho = padd[1]

    point_conditions = reshape(view(padd, 1 + 1 : 1 + 6*np), 6, np)

    element_loads = reshape(view(padd, 1 + 6*np + 1 : 1 + 6*np + 6*ne), 6, ne)

    u = padd[1 + 6*np + 6*ne + 1]
    v = padd[1 + 6*np + 6*ne + 2]
    w = padd[1 + 6*np + 6*ne + 3]
    p = padd[1 + 6*np + 6*ne + 4]
    q = padd[1 + 6*np + 6*ne + 5]
    r = padd[1 + 6*np + 6*ne + 6]
    delta = view(padd, 1 + 6*np + 6*ne + 6 + 1 : 1 + 6*np + 6*ne + 6 + NG)

    return (rho = rho, point_conditions = point_conditions, element_loads = element_loads,
        u = u, v = v, w = w, p = p, q = q, r = r, delta = delta)
end


text:
using Base.Collections

"Generate k-nearest neighborhood graph with distances"
function find_nn{T<:Real}(X::AbstractMatrix{T}, k::Int=12; excluding=true)
    m, n = size(X)
    r = Array(T, (n, n))
    d = Array(T, excluding ? k : k+1, n)
    e = Array(Int, excluding ? k : k+1, n)

    At_mul_B!(r, X, X)
    sa2 = sum(X.^2, 1)
    idx_range = excluding ? (2:k+1) : (1:k+1)

    @inbounds for j = 1 : n
        @inbounds for i = 1 : j-1
             r[i,j] = r[j,i]
        end
        r[j,j] = 0
        @inbounds for i = j+1 : n
            v = sa2[i] + sa2[j] - 2 * r[i,j]
            r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
        e[:, j] = sortperm(r[:,j])[idx_range]
        d[:, j] = r[e[:, j],j]
    end

    return (d, e)
end

"Find connected components for a undirected graph given its adjacency matrix"
function components(E::AbstractMatrix{Int})
    m, n = size(E)
    cmap = zeros(Int, n)
    cc = Array(Vector{Int}, 0)
    queue = Int[]

    for v in 1 : n
        if cmap[v] == 0
            # Start BFS with coloring
            c = Int[]
            push!(queue, v)
            while !isempty(queue)
                w = shift!(queue)
                for k in E[:, w]
                    if cmap[k] == 0
                        cmap[k] = 1
                        push!(queue, k)
                    end
                end
                cmap[w] = 2
                push!(c, w)
            end
            # Save component
            push!(cc, c)
        end
    end

    return cc
end

"Dijkstra's algorithm for single-source shortest path"
function dijkstra{T<:Real}(D::AbstractMatrix{T}, E::AbstractMatrix{Int}, src::Int, dst::Int =-1)
    m, n = size(D)
    path = zeros(Int, n)
    dist = fill(Inf, n)
    dist[src] = 0

    q = PriorityQueue(1:n, dist)
    while !isempty(q)
        u = dequeue!(q)
        if u == dst
            break
        end
        for i in 1:m
            v = E[i,u]
            alt = dist[u] + D[i,u]
            if haskey(q, v) && alt < q[v]
                dist[v] = q[v] = alt
                path[v] = u
            end
        end
    end

    return (path, dist)
end

"Find strongly connected components for undirected graph given its adjacency matrix"
function scomponents(E::AbstractMatrix{Int})

    function tarjan(v::Int)
        global E, I, index, lowlink, stack, components
        index[v] = I
        lowlink[v] = I
        I += 1
        push!(stack, v)

        for w in E[:,v]
            if index[w] == -1
                tarjan(w)
                lowlink[v] = min(lowlink[v], lowlink[w])
            elseif w in stack
                lowlink[v] = min(lowlink[v], lowlink[w])
            end
        end

        if lowlink[v] == index[v]
            component = Int[]
            while true
                w = pop!(stack)
                push!(component, w)
                if w == v
                    break
                end
            end
            push!(components, component)
        end
    end

    m, n = size(E)
    I = 0
    index = fill(-1, n)
    lowlink = fill(-1, n)
    stack = Int[]
    components = Array{Int,1}[]

    for v = 1 : n
        if index[v] == -1
            tarjan(v)
        end
    end

    return components
end

"Generate a swiss roll dataset"
function swiss_roll(n::Int = 1000, noise::Float64=0.05)
    t = (3 * pi / 2) * (1 .+ 2 * rand(n, 1))
    height = 30 * rand(n, 1)
    X = [t .* cos(t) height t .* sin(t)] + noise * randn(n, 3)
    labels = vec(rem(sum([round(Int, t / 2) round(Int, height / 12)], 2), 2))
    return X', labels
end

"Perform spectral decomposition for Ax=λI"
function decompose{T<:Real}(M::AbstractMatrix{T}, d::Int)
    W = isa(M, AbstractSparseArray) ? Symmetric(full(M)) : Symmetric(M)
    F = eigfact!(W)
    idx = sortperm(F[:values])[2:d+1]
    return F[:values][idx], F[:vectors][:,idx]
end

"Perform spectral decomposition for Ax=λB"
function decompose{T<:Real}(A::AbstractMatrix{T}, B::AbstractMatrix{T}, d::Int)
    AA = isa(A, AbstractSparseArray) ? Symmetric(full(A)) : Symmetric(A)
    BB = isa(B, AbstractSparseArray) ? Symmetric(full(B)) : Symmetric(B)
    F = eigfact(AA, BB)
    idx = sortperm(F[:values])[2:d+1]
    return F[:values][idx], F[:vectors][:,idx]
end


text:
using WeakRefStrings
using Base.Test

if !isdefined(Core, :String)
    typealias String UTF8String
end

data = "hey there sailor".data

str = WeakRefString(pointer(data), 3)
str2 = WeakRefString(pointer(data), 3, 1)

@test length(str) == 3
for (i,c) in enumerate(str)
    @test data[i] == c % UInt8
end
@test string(str) == "hey"
@test String(str) == "hey"

@test str == str2

io = IOBuffer()
show(io, str)
@test takebuf_string(io) == "\"hey\""

show(io, typeof(str))
@test takebuf_string(io) == "WeakRefString{UInt8}"

# simulate UTF16 data
data = reinterpret(UInt16, [0x68, 0x00, 0x65, 0x00, 0x79, 0x00])

str = WeakRefString(pointer(data), 3)
@test typeof(str) == WeakRefString{UInt16}
@test String(str) == "hey"
@test length(str) == 3
for (i,c) in enumerate(str)
    @test data[i] == c % UInt8
end

# simulate UTF32 data with trailing null byte
data = reinterpret(UInt32, [0x68, 0x00, 0x00, 0x00, 0x65, 0x00, 0x00, 0x00, 0x79, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00])

str = WeakRefString(pointer(data), 4)
@test typeof(str) == WeakRefString{UInt32}
@test String(str) == "hey"
@test length(str) == 4
for (i,c) in enumerate(str)
    @test data[i] == c % UInt8
end


text:
#===============================================================================

The palindromic number 595 is interesting because it can be written as the sum
of consecutive squares: 6^2 + 7^2 + 8^2 + 9^2 + 10^2 + 11^2 + 12^2.

There are exactly eleven palindromes below one-thousand that can be written as
consecutive square sums, and the sum of these palindromes is 4164. Note that 1 =
0^2 + 1^2 has not been included as this problem is concerned with the squares of
positive integers.

Find the sum of all the numbers less than 108 that are both palindromic and can
be written as the sum of consecutive squares.

--------------------------------------------------------------------------------

The solution is easy: we generate all i^2 + ... + j^2 below the limit and we
check that the number obtained is smaller than the limit.

Trap: a same number can be the sum of different sums of squares...

===============================================================================#

include("Euler.jl")

function euler(limit)
    i = 1
    ns = Set{Int}()

    while true
        n = i^2
        if limit <= n
            break
        end

        j = i+1
        n += j^2

        while n < limit
            if is_palindromic(n)
                push!(ns, n)
            end
            j += 1
            n += j^2
        end

        i += 1
    end

    sum(ns)
end

println("10^3: $(euler(10^3))")
println("10^8: $(euler(10^8))")

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
"""
This creates an enum which details how extrapolation from the interpolation domain should be done.
"""
@enum Schumaker_ExtrapolationSchemes begin
    Curve = 0
    Linear = 1
    Constant = 2
end

"""
    Schumaker(x::Array{T,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
    Schumaker(x::Array{Int,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
    Schumaker(x::Array{Date,1},y::Array{T,1} ; gradients::Union{Missing,Array{T,1}} = missing, extrapolation::Schumaker_ExtrapolationSchemes = Curve,
                  left_gradient::Union{Missing,T} = missing, right_gradient::Union{Missing,T} = missing)
Creates a Schumaker spline.
### Takes
* x - A vector of x coordinates.
* y - A vector of y coordinates.
* extrapolation (optional) - This should be Curve, Linear or Constant specifying how to interpolate outside of the sample domain. By default it is Curve which extends out the first and last quadratic curves. The other options are Linear which extends the line (from first and last curve) out from the first and last point and Constant which extends out the y value at the first and last point.
* gradients (optional)- A vector of gradients at each point. If not supplied these are imputed from x and y.
* left_gradient - The gradient at the lowest value of x in the domain. This will override the gradient imputed or submitted in the gradients optional argument (if it is submitted there)
* right_gradient - The gradient at the highest value of x in the domain. This will override the gradient imputed or submitted in the gradients optional argument (if it is submitted there)

### Returns
* A Schumaker object which details the spline. This object can then be evaluated with evaluate or evaluate_integral.
 """
struct Schumaker{T<:AbstractFloat}
    IntStarts_::Array{T,1}
    coefficient_matrix_::Array{T,2}
    function Schumaker(x::Array{<:Real,1},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, left_gradient::Union{Missing,<:Real} = missing, right_gradient::Union{Missing,<:Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve))
        T = typeof(AbstractFloat(0.0)) # TODO This is pretty dodgy. Try to not have to do this.
        if length(x) == 0
            error("Zero length x vector is insufficient to create Schumaker Spline.")
        elseif length(x) == 1
            IntStarts = Array{T,1}(x)
            SpCoefs = [0 0 y[1]]
            # Note that this hardcodes in constant extrapolation. This is only
            # feasible one as we do not have derivative or curve information.
            return new{T}(IntStarts, SpCoefs)
        elseif length(x) == 2
            IntStarts = Array{T,1}([x[1]])
            IntEnds   = Array{T,1}([x[2]])
            linear_coefficient = (y[2]- y[1]) / (x[2]-x[1])
            SpCoefs = [0 linear_coefficient y[1]]
            # In this case it defaults to curve extrapolation (which is same as linear here)
            # So we just alter in case constant is specified.
            if extrapolation[1] == Constant || extrapolation[2] == Constant
                matrix_without_extrapolation = hcat(IntStarts , SpCoefs)
                matrix_with_extrapolation    = extrapolate(matrix_without_extrapolation, extrapolation, x[2], y)
                return new{T}(matrix_with_extrapolation[:,1], matrix_with_extrapolation[:,2:4])
            else
                return new{T}(IntStarts, SpCoefs)
            end
        end
        if ismissing(gradients)
           gradients = imputeGradients(x,y)
        end
        if !ismissing(left_gradient)
            gradients[1] = left_gradient
        end
        if !ismissing(right_gradient)
            gradients[length(gradients)] = right_gradient
        end
        IntStarts, SpCoefs = getCoefficientMatrix(x,y,gradients, extrapolation)
        return new{T}(IntStarts, SpCoefs)
     end
    function Schumaker(IntStarts_::Array{T,1}, coefficient_matrix_::Array{T,2}) where {T<:Real}
        return new{T}(IntStarts_, coefficient_matrix_)
    end
    function Schumaker{Q}(IntStarts_::Array{<:Real,1}, coefficient_matrix_::Array{<:Real,2}) where Q<:Real
        return new{Q}(Q.(IntStarts_), Q.(coefficient_matrix_))
    end
    function Schumaker(x::Array{Date,1},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, left_gradient::Union{Missing,<:Real} = missing, right_gradient::Union{Missing,<:Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve))
        days_as_ints = Dates.days.(x)
        return Schumaker(days_as_ints , y; gradients = gradients , extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
    end
    function Schumaker(x::Union{Array{T,1},Array{Union{Missing,T},1}},y::Union{Array{R,1},Array{Union{Missing,R},1}} ; gradients::Union{Missing,Array{<:Real,1}} = missing, left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where T<:Real where R<:Real
        got_both = (!).(ismissing.(x) .| ismissing.(y))
        promo_type = promote_type(R,T)
        new_x = convert.(promo_type, x[got_both])
        new_y = convert.(promo_type, y[got_both])
        new_gradients = ismissing(gradients)        ? missing : convert.(promo_type, gradients[got_both])
        converted_left  = ismissing(left_gradient)  ? missing : convert(promo_type, left_gradient)
        converted_right = ismissing(right_gradient) ? missing : convert(promo_type, right_gradient)
        new_left  = got_both[1]                     ? converted_left  : missing
        new_right = got_both[length(got_both)]      ? converted_right : missing
        if length(new_x) == 0 error("After removing missing elements there are no points left to estimate schumaker spline") end
        return Schumaker(new_x , new_y; gradients = new_gradients , extrapolation = extrapolation, left_gradient = new_left, right_gradient = new_right)
    end
    function Schumaker{Q}(x::AbstractArray, y::AbstractArray ; gradients::Union{Missing,AbstractArray} = missing, left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing,
                          extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where Q<:Real
        got_both = (!).(ismissing.(x) .| ismissing.(y))
        new_x = convert.(Q, x[got_both])
        new_y = convert.(Q, y[got_both])
        new_gradients = ismissing(gradients)        ? missing : convert.(Q, gradients[got_both])
        converted_left  = ismissing(left_gradient)  ? missing : convert(Q, left_gradient)
        converted_right = ismissing(right_gradient) ? missing : convert(Q, right_gradient)
        new_left  = got_both[1]                     ? converted_left  : missing
        new_right = got_both[length(got_both)]      ? converted_right : missing
        if length(new_x) == 0 error("After removing missing elements there are no points left to estimate schumaker spline") end
        return Schumaker(new_x , new_y; gradients = new_gradients , extrapolation = extrapolation, left_gradient = new_left, right_gradient = new_right)
    end
    function Schumaker(x::Union{AbstractArray{T,1},AbstractArray{Union{Missing,T},1}},y::Union{AbstractArray{R,1},AbstractArray{Union{Missing,R},1}} ; gradients::Union{Missing,AbstractArray{<:Real,1}} = missing, left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing,
                       extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve)) where T<:Real where R<:Real
       promo_type = promote_type(T,R)
       return Schumaker{promo_type}(x, y; gradients = gradients, extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
   end
end
Base.broadcastable(e::Schumaker) = Ref(e)

"""
Evaluates the spline at a point. The point can be specified as a Real number (Int, Float, etc) or a Date.
Derivatives can also be taken.
### Takes
 * spline - A Schumaker type spline
 * PointToExamine - The point at which to evaluate the integral
 * derivative - The derivative being sought. This should be 0 to just evaluate the spline, 1 for the first derivative or 2 for a second derivative.
 Higher derivatives are all zero (because it is a quadratic spline). Negative values do not give integrals. Use evaluate_integral instead.
### Returns
 * A value of the spline or appropriate derivative in the same format as specified in the spline.
"""
function evaluate(spline::Schumaker, PointToExamine::T,  derivative::Integer = 0) where T<:Real
    # Derivative of 0 means normal spline evaluation.
    # Derivative of 1, 2 are first and second derivatives respectively.
    IntervalNum = searchsortedlast(spline.IntStarts_, PointToExamine)
    IntervalNum = max(IntervalNum, 1)
    xmt = PointToExamine - spline.IntStarts_[IntervalNum]
    Coefs = spline.coefficient_matrix_[IntervalNum,:]
    if derivative == 0
        return reshape(Coefs' * [xmt^2 xmt 1]', 1)[1]
    elseif derivative == 1
        return reshape(Coefs' * [2*xmt 1 0]', 1)[1]
    elseif derivative == 2
        return reshape(Coefs' * [2 0 0]', 1)[1]
    elseif derivative < 0
        error("This function cannot do integrals. Use evaluate_integral instead")
    else
        return 0.0
    end
end
function evaluate(spline::Schumaker, PointToExamine::Date,  derivative::Int = 0)
    days_as_int = Dates.days.(PointToExamine)
    return evaluate(spline,days_as_int,  derivative)
end
function (s::Schumaker)(x::Union{Real,Date}, deriv::Integer = 0)
    return evaluate(s, x, deriv)
end

"""
Estimates the integral of the spline between lhs and rhs. These end points can be input
as Reals or Dates.
### Takes
 * spline - A Schumaker type spline
 * lhs - The left hand limit of the integral
 * rhs - The right hand limit of the integral

### Returns
 * A Float64 value of the integral.
"""
function evaluate_integral(spline::Schumaker, lhs::Real, rhs::Real)
    first_interval = searchsortedlast(spline.IntStarts_, lhs)
    last_interval = searchsortedlast(spline.IntStarts_, rhs)
    number_of_intervals = last_interval - first_interval
    if number_of_intervals == 0
        return section_integral(spline , lhs, rhs)
    elseif number_of_intervals == 1
        first = section_integral(spline , lhs , spline.IntStarts_[first_interval + 1])
        last  = section_integral(spline , spline.IntStarts_[last_interval]  , rhs)
        return first + last
    else
        interior_areas = 0.0
        first = section_integral(spline , lhs , spline.IntStarts_[first_interval + 1])
        for i in 1:(number_of_intervals-1)
            interior_areas = interior_areas + section_integral(spline ,  spline.IntStarts_[first_interval + i] , spline.IntStarts_[first_interval + i+1] )
        end
        last  = section_integral(spline , spline.IntStarts_[last_interval]  , rhs)
        return first + interior_areas + last
    end
end
function evaluate_integral(spline::Schumaker, lhs::Date, rhs::Date)
    return evaluate_integral(spline, Dates.days.(lhs) , Dates.days.(rhs))
end
function section_integral(spline::Schumaker, lhs::Real,  rhs::Real)
    # Note that the lhs is used to infer the interval.
    IntervalNum = searchsortedlast(spline.IntStarts_, lhs)
    IntervalNum = max(IntervalNum, 1)
    Coefs = spline.coefficient_matrix_[ IntervalNum , :]
    r_xmt = rhs - spline.IntStarts_[IntervalNum]
    l_xmt = lhs - spline.IntStarts_[IntervalNum]
    return reshape(Coefs' * [(1/3)*r_xmt^3 0.5*r_xmt^2 r_xmt]', 1)[1] - reshape(Coefs' * [(1/3)*l_xmt^3 0.5*l_xmt^2 l_xmt]', 1)[1]
end

"""
find_derivative_spline(spline::Schumaker)
    Returns a SchumakerSpline that is the derivative of the input spline
"""
function find_derivative_spline(spline::Schumaker)
    coefficient_matrix = Array{Float64,2}(undef, size(spline.coefficient_matrix_)...   )
    coefficient_matrix[:,3] .= spline.coefficient_matrix_[:,2]
    coefficient_matrix[:,2] .= 2 .* spline.coefficient_matrix_[:,1]
    coefficient_matrix[:,1] .= 0.0
    return Schumaker(spline.IntStarts_, coefficient_matrix)
end

"""
    imputeGradients(x::Array{T,1}, y::Array{T,1})
Imputes gradients based on a vector of x and y coordinates.
"""
function imputeGradients(x::Array{<:Real,1}, y::Array{<:Real,1})
     n = length(x)
     # Judd (1998), page 233, second last equation
     L = sqrt.( (x[2:n]-x[1:(n-1)]).^2 + (y[2:n]-y[1:(n-1)]).^2)
     # Judd (1998), page 233, last equation
     d = (y[2:n]-y[1:(n-1)])./(x[2:n]-x[1:(n-1)])
     # Judd (1998), page 234, Eqn 6.11.6
     Conditionsi = d[1:(n-2)].*d[2:(n-1)] .> 0
     MiddleSiwithoutApplyingCondition = (L[1:(n-2)].*d[1:(n-2)]+L[2:(n-1)].* d[2:(n-1)]) ./ (L[1:(n-2)]+L[2:(n-1)])
     sb = Conditionsi .* MiddleSiwithoutApplyingCondition
     # Judd (1998), page 234, Second Equation line plus 6.11.6 gives this array of slopes.
     ff = [((-sb[1]+3*d[1])/2);  sb ;  ((3*d[n-1]-sb[n-2])/2)]
     return ff
 end

"""
Splits an interval into 2 subintervals and creates the quadratic coefficients
### Takes
 * s - A 2 entry vector with gradients at either end of the interval
 * z - A 2 entry vector with y values at either end of the interval
 * Smallt - A 2 entry vector with x values at either end of the interval

### Returns
 * A 2 x 5 matrix. The first column is the x values of start of the two subintervals. The second column is the ends. The last 3 columns are quadratic coefficients in two subintervals.
"""
function schumakerIndInterval(s::Array{<:Real,1}, z::Array{<:Real,1}, Smallt::Array{<:Real,1}) where T<:Real
   # The SchumakerIndInterval function takes in each interval individually
   # and returns the location of the knot as well as the quadratic coefficients in each subinterval.

   # Judd (1998), page 232, Lemma 6.11.1 provides this if condition:
   if (sum(s)*(Smallt[2]-Smallt[1]) == 2*(z[2]-z[1]))
     tsi = Smallt[2]
   else
     # Judd (1998), page 233, Algorithm 6.3 along with equations 6.11.4 and 6.11.5 provide this whole section
     delta = (z[2] -z[1])/(Smallt[2]-Smallt[1])
     Condition = ((s[1]-delta)*(s[2]-delta) >= 0)
     Condition2 = abs(s[2]-delta) < abs(s[1]-delta)
     if (Condition)
       tsi = sum(Smallt)/2
     elseif (Condition2)
       tsi = (Smallt[1] + (Smallt[2]-Smallt[1])*(s[2]-delta)/(s[2]-s[1]))
     else
       tsi = (Smallt[2] + (Smallt[2]-Smallt[1])*(s[1]-delta)/(s[2]-s[1]))
     end
   end

   # Judd (1998), page 232, 3rd last equation of page.
   alpha = tsi-Smallt[1]
   beta = Smallt[2]-tsi
   # Judd (1998), page 232, 4th last equation of page.
   sbar = (2*(z[2]-z[1])-(alpha*s[1]+beta*s[2]))/(Smallt[2]-Smallt[1])
   # Judd (1998), page 232, 3rd equation of page. (C1, B1, A1)
   Coeffs1 = [ (sbar-s[1])/(2*alpha)  s[1]  z[1] ]
   if (beta == 0)
     Coeffs2 = Coeffs1
   else
     # Judd (1998), page 232, 4th equation of page. (C2, B2, A2)
     Coeffs2 = [ (s[2]-sbar)/(2*beta)  sbar  Coeffs1 * [alpha^2, alpha, 1] ]
   end
   Machine4Epsilon = 4*eps()
     if (tsi  <  Smallt[1] + Machine4Epsilon )
         return [Smallt[1] Coeffs2]
     elseif (tsi + Machine4Epsilon > Smallt[2] )
         return [Smallt[1]  Coeffs1]
     else
         return [Smallt[1] Coeffs1 ; tsi Coeffs2]
     end
 end

 """
 Calls SchumakerIndInterval many times to get full set of spline intervals and coefficients. Then calls extrapolation for out of sample behaviour
### Takes
 * gradients - A vector of gradients at each point
 * x - A vector of x coordinates
 * y - A vector of y coordinates
 * extrapolation - A string in ("Curve", "Linear", "Constant") that gives behaviour outside of interpolation range.

### Returns
 * A vector of interval starts
 * A vector of interval ends
 * A matrix of all coefficients
  """
 function getCoefficientMatrix(x::Array{<:Real,1}, y::Array{<:Real,1}, gradients::Array{<:Real,1},
                               extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes})
   n = length(x)
   fullMatrix = schumakerIndInterval([gradients[1], gradients[2]], [y[1], y[2]], [x[1], x[2]] )
    for intrval = 2:(n-1)
      Smallt = [ x[intrval] , x[intrval + 1] ]
      s = [ y[intrval], y[intrval + 1] ]
      z = [ gradients[intrval], gradients[intrval + 1] ]
      intMatrix = schumakerIndInterval(z,s,Smallt)
      fullMatrix = vcat(fullMatrix,intMatrix)
    end
    fullMatrix = extrapolate(fullMatrix, extrapolation, x[n], y)
   return fullMatrix[:,1], fullMatrix[:,2:4]
 end

"""
 Adds a row on top and bottom of coefficient matrix to give out of sample prediction.
### Takes
 * fullMatrix - output from GetCoefficientMatrix first few lines
 * extrapolation - A tuple with two enums in (Curve, Linear, Constant) that gives behaviour outside of interpolation range.
 * x - A vector of x coordinates
 * y - A vector of y coordinates

### Returns
  * A new version of fullMatrix with out of sample prediction built into it.
"""
function extrapolate(fullMatrix::Array{<:Real,2}, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes}, Topx::Real, y::Array{<:Real,1})
    # In the fullMatrix the first column are the x values and then the next 3 columns are a, b, c in the expression a(x-start)^2 + b(x-start) + c
  if (extrapolation[1] == Curve) && (extrapolation[2] == Curve)
    return fullMatrix
  end
  # Preliminaries used throughout
  dim = size(fullMatrix)[1]
  Botx   = fullMatrix[1,1]
  Boty   = y[1]
  Topy = y[length(y)]
  # Initialising variable so their domain is not restricted to the if statement spaces.
  BotB = 0.0
  BotC = 0.0
  TopA = 0.0
  TopB = 0.0
  TopC = 0.0
  # Now creating the extrapolation to the left.
  if extrapolation[1] == Linear
    BotB = fullMatrix[1 , 3]
    BotC   = Boty - BotB
  elseif extrapolation[1] == Constant
    BotB = 0.0
    BotC = Boty
  end # Note for the curve case we will simply not append the new block.
  BotRow = [Botx-1e-10, 0.0, BotB, BotC]

  # Now doing the extrapolation to the right.
  if extrapolation[2] == Linear # This is a bit more complicated than before because the
                                # coefficients by themselves give the gradients at the left
                                # of the interval. Here we want the gradient at the right.
    last_interval_width = Topx - fullMatrix[dim , 1]
    grad_at_right = 2 * fullMatrix[dim , 2] * last_interval_width + fullMatrix[dim , 3]
    TopB = grad_at_right
    TopC = Topy
  elseif extrapolation[2] == Constant
    TopB = 0.0
    TopC = Topy
  elseif extrapolation[2] == Curve # We are just going to add this one on regardless as otherwise end of data interval information is lost.
    Gap  =  Topx - fullMatrix[dim ,1]
    TopA = fullMatrix[dim ,2]
    TopB = fullMatrix[dim ,3] + 2 * TopA * Gap
    TopC = fullMatrix[dim ,4] + TopB * Gap - TopA * Gap*Gap
  end
  TopRow = [ Topx, TopA ,TopB ,TopC]
  # Appending blocks and returning.
  if extrapolation[1] != Curve fullMatrix = vcat(BotRow' , fullMatrix) end
  fullMatrix = vcat(fullMatrix,  TopRow')
  return fullMatrix
end


text:
abstract type AtmosLinearModel <: BalanceLaw
end

vars_state(lm::AtmosLinearModel, FT) = vars_state(lm.atmos,FT)
vars_gradient(lm::AtmosLinearModel, FT) = @vars()
vars_diffusive(lm::AtmosLinearModel, FT) = @vars()
vars_aux(lm::AtmosLinearModel, FT) = vars_aux(lm.atmos,FT)
vars_integrals(lm::AtmosLinearModel, FT) = @vars()

update_aux!(dg::DGModel, lm::AtmosLinearModel, Q::MPIStateArray, auxstate::MPIStateArray, t::Real) = nothing
integrate_aux!(lm::AtmosLinearModel, integ::Vars, state::Vars, aux::Vars) = nothing
flux_diffusive!(lm::AtmosLinearModel, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real) = nothing
function wavespeed(lm::AtmosLinearModel, nM, state::Vars, aux::Vars, t::Real)
  ref = aux.ref_state
  return soundspeed_air(ref.T)
end

function boundary_state!(nf::Rusanov, lm::AtmosLinearModel, x...)
  atmos_boundary_state!(nf, NoFluxBC(), lm.atmos, x...)
end
function boundary_state!(nf::CentralNumericalFluxDiffusive, lm::AtmosLinearModel, x...)
  nothing
end
init_aux!(lm::AtmosLinearModel, aux::Vars, geom::LocalGeometry) = nothing
init_state!(lm::AtmosLinearModel, state::Vars, aux::Vars, coords, t) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
  atmos::M
end
function flux_nondiffusive!(lm::AtmosAcousticLinearModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  ref = aux.ref_state
  e_pot = gravitational_potential(lm.atmos.orientation, aux)

  flux.ρ = state.ρu
  # TODO: use MoistThermodynamics.linearized_air_pressure
  # need to avoid dividing then multiplying by ρ
  pL = state.ρ * FT(R_d) * FT(T_0) + FT(R_d) / FT(cv_d) * (state.ρe - state.ρ * e_pot)
  flux.ρu += pL*I
  flux.ρe = ((ref.ρe + ref.p)/ref.ρ - e_pot)*state.ρu
  nothing
end
function source!(lm::AtmosAcousticLinearModel, source::Vars, state::Vars, aux::Vars, t::Real)
  nothing
end

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
  atmos::M
end
function flux_nondiffusive!(lm::AtmosAcousticGravityLinearModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  ref = aux.ref_state
  e_pot = gravitational_potential(lm.atmos.orientation, aux)

  flux.ρ = state.ρu
  pL = state.ρ * FT(R_d) * FT(T_0) + FT(R_d) / FT(cv_d) * (state.ρe - state.ρ * e_pot)
  flux.ρu += pL*I
  flux.ρe = ((ref.ρe + ref.p)/ref.ρ)*state.ρu
  nothing
end
function source!(lm::AtmosAcousticGravityLinearModel, source::Vars, state::Vars, aux::Vars, t::Real)
  ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
  source.ρu = state.ρ * ∇Φ
  nothing
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#2 - echelon system with production

using Distributions
using InventoryManagement

#define network connectivity
adj_matrix = [0 0 1;
              0 0 1;
              0 0 0]
net = MetaDiGraph(adj_matrix) # 1 (plant) & 2 (distribution center) supply 3 (market)
materials = [:A, :B]
bom = [0 0; # B -> A
      -1 0]
set_prop!(net, :materials, materials)
set_prop!(net, :bill_of_materials, bom)

#specify parameters, holding costs and capacity, market demands and penalty for unfilfilled demand
set_props!(net, 1, Dict(:initial_inventory => Dict(:A => 75, :B => 100),
                        :holding_cost => Dict(:A => 0, :B => 0),
                        :production_cost => Dict(:A => 0.01)))

set_props!(net, 2, Dict(:initial_inventory => Dict(:A => 125),
                        :holding_cost => Dict(:A => 0)))

set_props!(net, 3, Dict(:initial_inventory => Dict(:A => 100),
                        :holding_cost => Dict(:A => 0.01),
                        :demand_distribution => Dict(:A => Normal(5,0.5)),
                        :demand_frequency => Dict(:A => 2),
                        :sales_price => Dict(:A => 3),
                        :demand_penalty => Dict(:A => 0.01),
                        :supplier_priority => Dict(:A => [1,2])))

#specify sales prices, transportation costs, lead time
set_props!(net, 1, 3, Dict(:sales_price => Dict(:A => 2),
                          :transportation_cost => Dict(:A => 0.1),
                          :lead_time => Dict(:A => Poisson(3))))

set_props!(net, 2, 3, Dict(:sales_price => Dict(:A => 1),
                          :transportation_cost => Dict(:A => 0.1),
                          :lead_time => Dict(:A => Poisson(7))))

#define reorder policy parameters
policy_type = :sS #(s, S) policy
review_period = 1 #continuous review
s = Dict((3,:A) => 50) #lower bound on inventory
S = Dict((3,:A) => 100) #base stock level

#create environment and run simulation with reorder policy
num_periods = 100
env = SupplyChainEnv(net, num_periods, backlog = true, reallocate = true)
simulate_policy!(env, s, S; policy_type, review_period)

#make plots
using DataFrames, StatsPlots
#profit
node_profit = groupby(env.profit, :node)
profit = transform(node_profit, :value => cumsum)
fig1 = @df profit plot(:period, :value_cumsum, group={Node = :node}, legend = :topleft,
                    xlabel="period", ylabel="cumulative profit")

#inventory position
inv_position = filter(i -> i.node > 1 ? i.material == :A : i.material in [:A,:B], env.inv_position)
fig2 = @df inv_position plot(:period, :level, group={Node = :node, Mat = :material}, linetype=:steppost, legend = :bottomleft,
                    xlabel="period", ylabel="inventory position", yticks = 0:25:125)


text:
function PolynomialFunctionNewtonForm(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return PolynomialFunctionNewtonForm((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function degree(obj::PolynomialFunctionNewtonForm)
    return jcall(obj, "degree", jint, ())
end

function evaluate(::Type{PolynomialFunctionNewtonForm}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(PolynomialFunctionNewtonForm, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function get_centers(obj::PolynomialFunctionNewtonForm)
    return jcall(obj, "getCenters", Vector{jdouble}, ())
end

function get_coefficients(obj::PolynomialFunctionNewtonForm)
    return jcall(obj, "getCoefficients", Vector{jdouble}, ())
end

function get_newton_coefficients(obj::PolynomialFunctionNewtonForm)
    return jcall(obj, "getNewtonCoefficients", Vector{jdouble}, ())
end

function value(obj::PolynomialFunctionNewtonForm, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::PolynomialFunctionNewtonForm, arg0::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement,), arg0)
end

function value(obj::PolynomialFunctionNewtonForm, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
"""
Case 28:
This case studies a OMIB with a simple Periodic Variable Source with a single frequency component.
There is no fault, however, the system is not in steady state due to the changing source.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test28.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

tspan = (0.0, 1.0);
step = 1e-1
tsteps = tspan[1]:step:tspan[2]

@testset "Test 28 Periodic Variable Source ResidualModel" begin
    path = (joinpath(pwd(), "test-28"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, # system
            path,
            tspan,
        )

        #Solve problem
        @test execute!(sim, IDA(), saveat = tsteps) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)
        pvs = collect(get_components(PeriodicVariableSource, sys))[1]

        #get time domain results from PVS data
        Vt = zeros(length(tsteps)) .+ get_internal_voltage_bias(pvs)
        ω_V = get_internal_voltage_frequencies(pvs)
        for (ix, A) in enumerate(get_internal_voltage_coefficients(pvs))
            Vt += A[1] * sin.(ω_V[ix] .* tsteps) + A[2] * cos.(ω_V[ix] .* tsteps)
        end
        θt = zeros(length(tsteps)) .+ get_internal_angle_bias(pvs)
        ω_θ = get_internal_angle_frequencies(pvs)
        for (ix, A) in enumerate(get_internal_angle_coefficients(pvs))
            θt += A[1] * sin.(ω_θ[ix] .* tsteps) + A[2] * cos.(ω_θ[ix] .* tsteps)
        end

        #Obtain simulation data
        Vt_sim = get_state_series(results, ("InfBus", :Vt))
        θt_sim = get_state_series(results, ("InfBus", :θt))

        @test LinearAlgebra.norm(Vt .- Vt_sim[2]) <= 5e-3
        @test LinearAlgebra.norm(θt .- θt_sim[2]) <= 5e-3
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 28 Periodic Variable Source MassMatrixModel" begin
    path = (joinpath(pwd(), "test-28"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            sys, # system
            path,
            tspan,
        )

        #Solve problem
        @test execute!(sim, Rodas4(), saveat = tsteps) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)
        pvs = collect(get_components(PeriodicVariableSource, sys))[1]

        #get time domain results from PVS data
        Vt = zeros(length(tsteps)) .+ get_internal_voltage_bias(pvs)
        ω_V = get_internal_voltage_frequencies(pvs)
        for (ix, A) in enumerate(get_internal_voltage_coefficients(pvs))
            Vt += A[1] * sin.(ω_V[ix] .* tsteps) + A[2] * cos.(ω_V[ix] .* tsteps)
        end
        θt = zeros(length(tsteps)) .+ get_internal_angle_bias(pvs)
        ω_θ = get_internal_angle_frequencies(pvs)
        for (ix, A) in enumerate(get_internal_angle_coefficients(pvs))
            θt += A[1] * sin.(ω_θ[ix] .* tsteps) + A[2] * cos.(ω_θ[ix] .* tsteps)
        end

        #Obtain simulation data
        Vt_sim = get_state_series(results, ("InfBus", :Vt))
        θt_sim = get_state_series(results, ("InfBus", :θt))

        @test LinearAlgebra.norm(Vt .- Vt_sim[2]) <= 5e-3
        @test LinearAlgebra.norm(θt .- θt_sim[2]) <= 5e-3
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
BoundaryW(v1CG,v2CG,CG,Param) =
  BoundaryW(v1CG,v2CG,CG,Param.cache.dXdxIC)

function BoundaryW(v1CG,v2CG,CG,dXdxIC::AbstractArray)
wCG = -(dXdxIC[:,:,:,1,3,1].*v1CG[:,:,:,1]+
  dXdxIC[:,:,:,1,3,2].*v2CG[:,:,:,1])./
  dXdxIC[:,:,:,1,3,3];
  return wCG
end


text:
include("./src/version.jl")

vers = string("v", SYMATA_VERSION)
com = `git tag $vers`
run(com)


text:
"""

This example learns to classify images of fashion products(trousers, shirts, bags...)
from the [Fashion-MNIST](https://github.com/zalandoresearch/fashion-mnist) dataset.
There are 60000 training and 10000 test examples. Each input x
consists of 784 pixels representing a 28x28 image. The pixel values are
normalized to [0,1]. Each output y is converted to a ten-dimensional
one-hot vector (a vector that has a single non-zero component) indicating
the correct class (0-9) for a given image. 10 is used instead of 0.
Labels and descriptions are shown below.

    Label   Description
    1       T-shirt/top
    2       Trouser
    3       Pullover
    4       Dress
    5       Coat
    6       Sandal
    7       Shirt
    8       Sneaker
    9       Bag
    10      Ankle boot

You can run the demo using `julia fashion-mnist.jl` on the command line or
by first including `julia> include("fashion-mnist.jl")` and typing `julia> FashionMNIST.main()`
at the Julia prompt.  Options can be used like `julia fashion-mnist.jl --epochs 3`
or `julia> FashionMNIST.main("--epochs 3")`. Use `julia fashion-mnist.jl --help`
for a list of options.  The dataset will be automatically downloaded.
By default a softmax model will be trained for 10 epochs. You can also
train a multi-layer perceptron by specifying one or more --hidden sizes.
The accuracy for the training and test sets will be printed at every epoch
and optimized parameters will be returned.

"""
module FashionMNISTMLP
using Knet,CUDA,MLDatasets,ArgParse

function predict(w,x; pdrop=0)
    x = mat(x)
    for i=1:2:length(w)
        x = w[i]*dropout(x, pdrop) .+ w[i+1]
        if i<length(w)-1
            x = relu.(x) # max(0,x)
        end
    end
    return x
end

loss(w,x,ygold;pdrop=0) = nll(predict(w,x;pdrop=pdrop), ygold)

lossgradient = grad(loss)

function train(w, dtrn; lr=.5, epochs=10, pdrop=0)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = lossgradient(w, x, y; pdrop=pdrop)
            update!(w, g, lr=lr)
        end
    end
    return w
end

function weights(h...; atype=Array{Float32}, winit=0.1)
    w = Any[]
    x = 28*28
    for y in [h..., 10]
        push!(w, convert(atype, winit*randn(y,x)))
        push!(w, convert(atype, zeros(y, 1)))
        x = y
    end
    return w
end

function main(args="")
    s = ArgParseSettings()
    s.description="fashion-mnist.jl (c) 2017 Adapted by Emre Unal based on Deniz Yuret’s MNIST example https://github.com/denizyuret/Knet.jl/tree/master/examples/mnist-mlp/mlp.jl.\nMulti-layer perceptron model on the Fashion-MNIST dataset from https://github.com/zalandoresearch/fashion-mnist.\n"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=100; help="minibatch size")
        ("--epochs"; arg_type=Int; default=10; help="number of epochs for training")
        ("--hidden"; nargs='*'; arg_type=Int; help="sizes of hidden layers, e.g. --hidden 128 64 for a net with two hidden layers")
        ("--lr"; arg_type=Float64; default=0.15; help="learning rate")
        ("--winit"; arg_type=Float64; default=0.1; help="w initialized with winit*randn()")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--atype"; default="$(Knet.array_type[])"; help="array type: Array for cpu, KnetArray for gpu")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
        ("--dropout"; arg_type=Float64; default=0.5; help="Dropout probability.")
        # These are to experiment with sparse arrays
        # ("--xtype"; help="input array type: defaults to atype")
        # ("--ytype"; help="output array type: defaults to atype")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    o = parse_args(args, s; as_symbols=true)
    if !o[:fast]
        println(s.description)
        println("opts=",[(k,v) for (k,v) in o]...)
    end
    o[:seed] > 0 && Knet.seed!(o[:seed])
    atype = eval(Meta.parse(o[:atype]))
    w = weights(o[:hidden]...; atype=atype, winit=o[:winit])
    xtrn,ytrn = FashionMNIST.traindata()
    xtst,ytst = FashionMNIST.testdata()
    global dtrn = minibatch(xtrn, ytrn, o[:batchsize]; xtype=atype, xsize=(size(xtrn,1),size(xtrn,2),1,o[:batchsize]))
    global dtst = minibatch(xtst, ytst, o[:batchsize]; xtype=atype, xsize=(size(xtrn,1),size(xtrn,2),1,o[:batchsize]))
    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))
    if o[:fast]
        @time (train(w, dtrn; lr=o[:lr], epochs=o[:epochs], pdrop=o[:dropout]); CUDA.functional() && CUDA.synchronize())
    else
        report(0)
        @time for epoch=1:o[:epochs]
            train(w, dtrn; lr=o[:lr], epochs=1)
            report(epoch)
            if o[:gcheck] > 0
                gradcheck(loss, w, first(dtrn)...; gcheck=o[:gcheck], verbose=true)
            end
        end
    end
    return w
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia mnist.jl --epochs 10
# julia> FashionMNIST.main("--epochs 10")
PROGRAM_FILE == "fashion-mnist.jl" && main(ARGS)

end # module


text:
include("paper/N_plots.jl")
Random.seed!(7)

model, = TrajectoryOptimization.Dynamics.dubinscar!
n, m = model.n,model.m

# Objective
x0 = [2.5;2.5;0.]
xf = [7.5;2.5;0.]
tf =  3.0
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

r = 0.5
s1 = 30; s2 = 50; s3 = 15
circles = []
for i in range(0,stop=5,length=s1)
    push!(circles,(0.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(5.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(10.,i,r))
end
for i in range(0,stop=10,length=s2)
    push!(circles,(i,0.,r))
end
for i in range(0,stop=3,length=s3)
    push!(circles,(i,5.,r))
end
for i in range(5,stop=8,length=s3)
    push!(circles,(i,5.,r))
end

n_circles = 3*s1 + s2 + 2*s3

# plt = plot(title="Escape",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
# plot_obstacles(circles,:red)
# display(plt)

function cI(c,x,u)
    for i = 1:n_circles
        c[i] = circle_constraint(x,circles[i][1],circles[i][2],circles[i][3])
    end
    c
end

u_bnd = 6
obj_con_obstacles = TrajectoryOptimization.ConstrainedObjective(obj,cI=cI,u_min=-u_bnd,u_max=u_bnd)

opts = SolverOptions()
opts.resolve_feasible = false
opts.R_infeasible = 1.
opts.live_plotting = false
opts.verbose = false
N = 101

solver = Solver(model, obj_con_obstacles, integration=:rk3, N=N, opts=opts)
solver_foh = Solver(model, obj_con_obstacles, integration=:rk3_foh, N=N, opts=opts)
solver_dircol = Solver(model, obj_con_obstacles, integration=:rk3_foh, N=N)

# -Initial state and control trajectories
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = interp_rows(N, tf, Array(X_guess'))
U0 = rand(m,N)

@time res, stats = TrajectoryOptimization.solve(solver,X0,U0)
@time res_foh, stats_foh = TrajectoryOptimization.solve(solver_foh,X0,U0)

solver_dircol.opts.verbose = false
res_d,stats_d = solve_dircol(solver_dircol, X0, U0)

# Escape
plt = plot(title="Escape",aspect_ratio=:equal)
plot_obstacles(circles)
plot_trajectory!(to_array(res.X),width=2,color=:blue,label="rk3",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot_trajectory!(to_array(res_foh.X),width=2,color=:purple,label="rk3_foh",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot!(X0[1,:],X0[2,:],label="Infeasible Initialization",width=1,color=:purple,linestyle=:dash)
plot_trajectory!(res_d.X,width=2,color=:black,label="DIRCOL",linestyle=:dash)

eval_f, = gen_usrfun_ipopt(solver_dircol,:hermite_simpson)
var_zoh = DircolVars(res)
var_foh = DircolVars(res_foh)
eval_f(var_zoh.Z)
eval_f(var_foh.Z)
eval_f(res_d.Z)

@show stats["iterations"]
@show stats_foh["iterations"]

group = "dubinscar/escape"
solver_truth = Solver(model, obj_con_obstacles, N=501, integration=:rk3_foh)
dt_truth = solver_truth.dt
run_dircol_truth(solver_truth, Array(res_d.X), Array(res_d.U), group::String)

run_step_size_comparison(model, obj_con_obstacles, U0, group::String, [101,201,301]; integrations=[:rk3,:rk3_foh],dt_truth=dt_truth,opts=opts,infeasible=true,X0=X0)
plot_stat("runtime",group,legend=:topleft,title="Constrained Parallel Park")
plot_stat("iterations",group,legend=:left,title="Constrained Parallel Park")
plot_stat("error",group,yscale=:log10,legend=:right,title="Constrained Parallel Park")
plot_stat("c_max",group,yscale=:log10,legend=:right,title="Constrained Parallel Park")


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

using Markdown
using InteractiveUtils

# ╔═╡ 0f225480-3b85-11ec-2de2-f9c51a8642aa
md"First of all, thank you for having a look into Nbdev. If you are reading this then it's safe to say that you have taken Nbdev for a spin and you want to contribute to the project. May be you found somethign that can add value to Nbdev or may be something is missing which you want ot add. Whatever the reason maybe, here are some of the ways in which you can contribute 👇

* You can write a tutorial.
* You can contribute to documentations.
* You can provide feedback.
* You can contribute code.

Before moving into the details about how you can do the above mentioned contributions, let's see how to setup Nbdev for the first time."

# ╔═╡ 6fe29a04-31b0-4e32-ac20-80eac2693b2c
md"""## Setting things up

Working on Nbdev or any other open source project starts with a "fork". If you are new to open source or this is your first time contribution to an open source project then read the following [tutorial](https://docs.github.com/en/get-started/quickstart/fork-a-repo) to get a quick understanding of what a "fork" is and how you can get started with forking a project.

After you fork the project, clone the repo somewhere in your local machine.

Nbdev is built on itself i.e. the starting point of all the code and documentation in Nbdev is a Pluto notebook.

To get started, create a project directory and move into the directory. Here, you need to start a new environment like this 👇

```
mkdir ~/nbdev_dev
cd ~/nbdev_dev
julia
```

This would create a directory namedb "nbdev_dev" and take you inside that directory. After which the Julia terminal would fire up. Once inside Julia terminal type `]` to enter the package manager then type the following 👇

```
(v1.0) pkg> activate .
(nbdev_dev) pkg> dev <path/to/cloned/repo
```

After this move back to the Julia terminal and type 👇

```
julia> import Pluto
julia> Pluto.run()
```
This gives you a fresh Pluto environment to work in. From this point onwards, you can either start a new notebook or work on an existing notebook if it's a bug fix or documentation change.
"""

# ╔═╡ 518a5202-b1fc-4baf-b0ba-6b6ecd672ea8
md"## How to get started 🏃♀
Here are some of the way to get started with contributing 👇

* Pick a piece of code in nbdev and checkout the test for it. If there is a missing test then add it. If an existing test doesn't provide sufficient coverage then add it.

* You noticed somethign missing in the documentation? go ahead and add it.

* Feel like the tutorial needs imrpovement or a few more tutorials could be useful? go ahead and write a tutorial.
"

# ╔═╡ 4ab7ac54-6f5a-400e-91ce-1469e73e936c
md"## Oh! you found a bug 🖊🐛
* First ensure that such bug was not reported earlier in [issues](https://github.com/sapal6/nbdev.jl/issues)

* If you're unable to find an open issue addressing the problem, open a [new one](https://github.com/sapal6/nbdev.jl/issues/new). Be sure to include a title and clear description, as much relevant information as possible, and a code sample or an executable test case demonstrating the expected behavior that is not occurring.

* Be sure to add the complete error messages.
"

# ╔═╡ 95facb52-4681-4126-b128-8f769b3f2d61
md"## Intending to add a new feature or enhancing an old one? 👩🔧

* You can start a discussion on the [discussion tab](https://github.com/sapal6/nbdev.jl/discussions). this will help you to find if others are instrested in same thing as you and maybe you can find collaborators.

* While working on a new feature first create a detailed notebook about the feature that you would like to add. The notebook should contain a detailed description about the feature, why you think this feature is needed and a step-by-step description about what your code is doing. Once you think that you are ready with your example notebook, post a gist in the discussion forum.

* Once your feature request is dicussed and finalized then you can open a PR. Make sure to include a complete description of the feature and an example of it's usage. Make sure that your code has a good documentation and follows the style guide.

* Ensure that your PR includes tests that exercise not only your feature, but also any other code that might be impacted.
"

# ╔═╡ 7e6f9df9-1cc3-41c9-8cb7-24a8f86f75b9
md"Submitting a PR 🎈

* All Nbdev code starts with a Pluto notebook. The exported code from the notebook then appears as the final code under the src directory and the related documentation under the docs directory. If you are submitting your PR -->
  * Always start your developement with a Pluto notebook because if you modify code, add code or create a new code file under the src directory directly then the source notebooks won't have your modifications.

  * Once you add a new code in an existing notebook, update an existing code in a cell that is supposed to be exported or added a new code in a new notebook with cells marked with `#export` then make sure to run `notebook2script` in the same notebook so as to export the code to the respective files under the src directory.

  * After everything is finalized, make sure to run `make.jl docs` at the terminal in order to export the documentation.

* Make sure that your PR is accompanied by sufficient test coverage. Either in the notebook itself or under the test dicrectory. Currently Nbdev has poor test coverage so you may need to include test for the impacted code as well.

* Do not combine several unrelated fixes together. Create separate branches for separate fixes, separate features etc.

* Preserve the original code styling of an existing file.
"

# ╔═╡ 25981537-b077-44cb-8aed-078db3fa7fad
md"## Want to contribute to the docs? 📓
* The docs are generated from the notebooks in the`nbs` directory.

* If you are planning to edit an already existing notebook then make sure not to add the `#export` flag to the existing documentation cells or any new documentation cell that you might add.

* If you are planning to add a new notebook for documentation then make sure not to add `#export` to any of the cell in this new notebook.
"

# ╔═╡ Cell order:
# ╠═0f225480-3b85-11ec-2de2-f9c51a8642aa
# ╠═6fe29a04-31b0-4e32-ac20-80eac2693b2c
# ╠═518a5202-b1fc-4baf-b0ba-6b6ecd672ea8
# ╠═4ab7ac54-6f5a-400e-91ce-1469e73e936c
# ╠═95facb52-4681-4126-b128-8f769b3f2d61
# ╠═7e6f9df9-1cc3-41c9-8cb7-24a8f86f75b9
# ╠═25981537-b077-44cb-8aed-078db3fa7fad


text:
@testset "with_logabsdet_jacobian" begin
    derivative(f, x) = ChainRulesTestUtils.frule((ChainRulesTestUtils.NoTangent(), 1), f, x)[2]

    x = randexp()

    ChangesOfVariables.test_with_logabsdet_jacobian(log1pexp, x, derivative)
    ChangesOfVariables.test_with_logabsdet_jacobian(log1pexp, -x, derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(logexpm1, x, derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(log1mexp, -x, derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(log2mexp, log(2) - x, derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(logistic, -x, derivative)
    ChangesOfVariables.test_with_logabsdet_jacobian(logistic, x, derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(logit, rand(), derivative)

    ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, x, derivative)
    ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, -x, derivative)
end


text:
module SmartPlay

import PlayingCards: Card, createdeck

import Cribbage: Player
import Cribbage.RandomPlay: discardtocrib, pickcardtoplay
import Cribbage.Rules: validhands, validplays
import Cribbage.Scoring: scorehand, scoreplay

struct SmartPlayer <: Player
    name :: String
end

function SmartPlayer()
    return SmartPlayer(UUIDs.uuid4())
end

function averagescore(cards::Tuple{Vararg{Card}}, others::Array{Card,1})
    cards = collect(cards)
    scores = Int8[]
    for starter in others
        push!(scores, scorehand(cards, starter))
    end
    return sum(scores) / length(scores)
end

DECK = createdeck()

function discardtocrib(player::SmartPlayer, cards::Tuple{Vararg{Card}})
    subsets = validhands(cards)
    others = setdiff(DECK, cards)
    sort!(subsets, by=x -> averagescore(x, others), rev=true)
    hand = subsets[1]
    crib = setdiff(cards, hand)
    @debug "keeping $(string(hand)), discarding $(string(crib))"
    return hand, crib
end

function pickcardtoplay(player::SmartPlayer, prior::Tuple{Vararg{Card}}, cards::Tuple{Vararg{Card}})::Union{Missing,Card}
    possible_plays = validplays(prior, cards)
    if length(possible_plays) == 0
        return missing
    end
    if length(prior) == 0
        priorarray = Card[]
    else
        priorarray = collect(prior)
    end
    sort!(possible_plays, by=x -> scoreplay(priorarray, x), rev=true)
    return possible_plays[1]
end

end  # module SmartPlay


text:
#using Lint

function readSMAT(filename::AbstractString)
    (rows,header) = readdlm(filename;header=true)
    A = sparse(
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),1])+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),2])+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),3]),
               parse(Int,header[1]),
               parse(Int,header[2])
               )
    return A
end

type MatrixNetworkMetadata
    A::SparseMatrixCSC{Int64,Int64}
    labels::Vector{AbstractString}
    xy::Array{Float64,2}
    source::AbstractString
end

function load_matrix_network_all(name::AbstractString)
    A = load_matrix_network(name)
    pathname = joinpath(Pkg.dir("MatrixNetworks"),"data")

    meta_source = joinpath(pathname,"$(name).source")
    if isfile(meta_source)
        source = open(readall, meta_source)
    else
        source = "(None given)"
    end

    meta_xy = joinpath(pathname,"$(name).xy")
    if isfile(meta_xy)
        xy = readdlm(meta_xy)
    else
        xy = zeros(0,2)
    end

    meta_labels = joinpath(pathname,"$(name).labels")
    if isfile(meta_labels)
        labels = open(readlines, meta_labels)
    else
        labels = map(x -> @sprintf("%i",x), 1:size(A,1))
    end

    return MatrixNetworkMetadata(A,labels,xy,source)
end

function load_matrix_network(name::AbstractString)
    pathname = joinpath(Pkg.dir("MatrixNetworks"),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    if isfile(smatfile)
        return readSMAT(smatfile)
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function load_matrix_network_metadata(name::AbstractString)
    pathname = joinpath(Pkg.dir("MatrixNetworks"),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    meta_xy = joinpath(pathname,"$(name).xy")
    meta_labels = joinpath(pathname,"$(name).labels")
    if isfile(smatfile)
        if isfile(meta_xy) && isfile(meta_labels)
            xy = readdlm(meta_xy)
            labels = readdlm(meta_labels)
            return (readSMAT(smatfile),xy,labels)
        end
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function matrix_network_datasets()
    datasets_location = joinpath(Pkg.dir("MatrixNetworks"),"data")
    content = readdir(datasets_location)
    smat_files = filter(x->contains(x,".smat"),content)
    for i = 1:length(smat_files)
        # no need for this lintpragma anymore
        # @lintpragma( "Ignore use of undeclared variable end")
        smat_files[i] = smat_files[i][1:end-5]
    end
    return smat_files
end


text:
function load_matctr_instance!(db::SQLite.DB, name, scenario, source_type, out)
    println("Loading matctr: ($name, $scenario)")
    instance_name = source_type == "RAWGO" ? "$(name)_$(scenario)" : "$(name)"
    mat_path = abspath(joinpath(out, "$(instance_name)_sdp_mat.txt"))
    ctr_path = abspath(joinpath(out, "$(instance_name)_sdp_ctr.txt"))
    isfile_mat = isfile(mat_path)
    isfile_ctr = isfile(ctr_path)
    query = ""
    if isfile_mat && isfile_ctr
        query = """
        UPDATE instances SET mat_path='$mat_path', ctr_path='$ctr_path'
        WHERE name = '$name' AND scenario = $scenario
        """
    elseif  !isfile_mat || !isfile_ctr
        if isfile_mat
            println("No ctr file found.")
            query = """
            UPDATE instances SET mat_path='$mat_path'
            WHERE name = '$name' AND scenario = $scenario
            """
        elseif isfile_ctr
            println("No mat file found.")
            query = """
            UPDATE instances SET ctr_path = '$ctr_path'
            WHERE name = '$name' AND scenario = $scenario
            """
        else
            println("No ctr or mat file found.")
            return
        end
    end
    println(query)
    execute_query(db, query)
end

function load_matctr_instance_dfrow!(db::SQLite.DB, row, out)
    load_matctr_instance!(db, row[:name], row[:scenario], row[:source_type], out)
end

function load_matctr_instances!(db::SQLite.DB; out, recompute=false, subset=nothing, min_nv=typemin(Int), max_nv=typemax(Int))
    query = "SELECT name, scenario, source_type FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !recompute
        query *= " AND mat_path IS NULL OR ctr_path IS NULL"
    end
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Loading matctr paths: out=$out, recompute=$recompute")
    println("subset\n$subset\nend subset")
    load_matctr_instance_function!(row) = load_matctr_instance_dfrow!(db, row, out)
    load_matctr_instance_function!.(eachrow(results[!, [:name, :scenario, :source_type]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
export BrickletAnalogInV2VoltageCallbackThreshold
struct BrickletAnalogInV2VoltageCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAnalogInV2AnalogValueCallbackThreshold
struct BrickletAnalogInV2AnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAnalogInV2Identity
struct BrickletAnalogInV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAnalogInV2
"""
Measures DC voltage between 0V and 42V
"""
mutable struct BrickletAnalogInV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAnalogInV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_analog_in_v2")
        deviceInternal = package.BrickletAnalogInV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the measured voltage. The resolution is approximately 10mV.

If you want to get the voltage periodically, it is recommended to use the
:cb:`Voltage` callback and set the period with
:func:`Set Voltage Callback Period`.
"""
function get_voltage(device::BrickletAnalogInV2)
    return device.deviceInternal.get_voltage()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletAnalogInV2)
    return device.deviceInternal.get_analog_value()
end

export set_voltage_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Voltage` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Voltage` callback is only triggered if the voltage has changed since
the last triggering.
"""
function set_voltage_callback_period(device::BrickletAnalogInV2, period)
    device.deviceInternal.set_voltage_callback_period(period)
end

export get_voltage_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Voltage Callback Period`.
"""
function get_voltage_callback_period(device::BrickletAnalogInV2)
    return device.deviceInternal.get_voltage_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletAnalogInV2, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletAnalogInV2)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_voltage_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Voltage Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the voltage is *outside* the min and max values"
 "'i'",    "Callback is triggered when the voltage is *inside* the min and max values"
 "'<'",    "Callback is triggered when the voltage is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the voltage is greater than the min value (max is ignored)"
"""
function set_voltage_callback_threshold(device::BrickletAnalogInV2, option, min, max)
    device.deviceInternal.set_voltage_callback_threshold(option, min, max)
end

export get_voltage_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Voltage Callback Threshold`.
"""
function get_voltage_callback_threshold(device::BrickletAnalogInV2)
    return device.deviceInternal.get_voltage_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletAnalogInV2, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletAnalogInV2)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Voltage Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Voltage Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletAnalogInV2, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletAnalogInV2)
    return device.deviceInternal.get_debounce_period()
end

export set_moving_average
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the voltage.

Setting the length to 1 will turn the averaging off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletAnalogInV2, average)
    device.deviceInternal.set_moving_average(average)
end

export get_moving_average
"""
    $(SIGNATURES)

Returns the length of the moving average as set by :func:`Set Moving Average`.
"""
function get_moving_average(device::BrickletAnalogInV2)
    return device.deviceInternal.get_moving_average()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAnalogInV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAnalogInV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using ClassicalOrthogonalPolynomials, Test
import ClassicalOrthogonalPolynomials: orthogonalityweight

@testset "Laguerre" begin
    @testset "Laguerre weight" begin
        @test sum(LaguerreWeight()) ≡ 1.0
        @test LaguerreWeight()[0.1] ≈ orthogonalityweight(Laguerre())[0.1] ≈ exp(-0.1)
        @test LaguerreWeight(1/2)[0.1] ≈ sqrt(0.1)exp(-0.1)
    end

    @testset "Laguerre L" begin
        x = 0.1
        @test laguerrel.(0:3, 0, x) ≈ laguerrel.(0:3, x) ≈ [1, 1-x, (2 - 4x + x^2)/2, (6 - 18x + 9x^2 - x^3)/6]
        @test laguerrel.(0:3, 1/2, x) ≈ [1, 3/2 - x, (15 - 20x + 4x^2)/8, (105 - 210x + 84x^2 - 8x^3)/48]
    end
    @testset "Derivatives" begin
        L = Laguerre()
        D = Derivative(axes(L,1))
        x = 0.1
        @test (D*L)[x,1:4] ≈ [0,-1,x-2,-x^2/2 + 3x - 3]

        L = Laguerre(1/2)
        @test (D*L)[x,1:4] ≈ [0,-1,x-20/8,-x^2/2 + 7/2*x - 35/8]

        x = axes(L,1)
        X = L \ (x .* L)
        @test 0.1*L[0.1,1:10]' ≈ L[0.1,1:11]'*X[1:11,1:10]
    end
end


text:
using Distributions
using JSON

n = 15
p = 0.4
d = Binomial( n,p )

x = [ -3, -1, 0, 1, 3 ]

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("n", n),
	("p", p),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/partial.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

max_stars_repo_name: distributions-io/binomial-mgf

max_issues_repo_name: distributions-io/binomial-mgf

max_forks_repo_name: distributions-io/binomial-mgf