text:
module Extract
using Interpolations, DSP, AxisAlgorithms, SpecialMatrices
export PositiveThreshold, NegativeThreshold, DoubleThreshold, extract_wf, extract_wf_teo, noise_autocovariance

immutable PositiveThreshold{T<:Real}
    thr::T
end
apply_threshold(t::PositiveThreshold, x) = x >= t.thr
find_alignment_point(t::PositiveThreshold, v::AbstractVector) = indmax(v)

immutable NegativeThreshold{T<:Real}
    thr::T
    NegativeThreshold(thr) = new(-abs(thr))
end
NegativeThreshold(thr::Real) = NegativeThreshold{typeof(thr)}(thr)
apply_threshold(t::NegativeThreshold, x) = x <= t.thr
find_alignment_point(t::NegativeThreshold, v::AbstractVector) = indmin(v)

immutable DoubleThreshold{T<:Real}
    posthr::T
    negthr::T
    DoubleThreshold(thr) = new(abs(thr), -abs(thr))
    DoubleThreshold(posthr, negthr) = new(abs(posthr), -abs(negthr))
end
DoubleThreshold(posthr::Real, negthr::Real) = DoubleThreshold{Base.promote_typeof(posthr, negthr)}(posthr, negthr)
DoubleThreshold(thr::Real) = DoubleThreshold(thr, thr)
apply_threshold(t::DoubleThreshold, x) = (x >= t.posthr) | (x <= t.negthr)

function find_closest_extremum(x, istart, dt, firstallowed, lastallowed)
    irextremum = istart+zero(dt)
    sgn = sign(x[istart-dt] - x[istart])
    while irextremum < lastallowed && sign(x[irextremum] - x[irextremum+dt]) == sgn
        irextremum += dt
    end
    ilextremum = istart+zero(dt)
    sgn = sign(x[istart+dt] - x[istart])
    while ilextremum > firstallowed && sign(x[ilextremum] - x[ilextremum-dt]) == sgn
        ilextremum -= dt
    end
    (istart - ilextremum) < (irextremum - istart) ? ilextremum : irextremum
end

immutable FastInterpolationUpdater{T,TCoefs,S}
    interp::Interpolations.BSplineInterpolation{T,1,TCoefs,BSpline{Cubic{Flat}},OnGrid,1}
    prefiltering_system::S
    tmp::TCoefs
end

function FastInterpolationUpdater(x)
    interp = interpolate(x, BSpline(Cubic(Flat())), OnGrid())
    prefiltering_system = Interpolations.prefiltering_system(Interpolations.tweight(x), Interpolations.tcoef(x), length(interp.coefs), Cubic{Flat}, OnGrid)
    FastInterpolationUpdater(interp, prefiltering_system, similar(interp.coefs))
end

function update_interpolation!(updater::FastInterpolationUpdater, x, rg)
    interp = updater.interp
    copy!(interp.coefs, 2:length(rg)+1, x, rg)
    interp.coefs[1] = interp.coefs[end] = 0
    W, b = updater.prefiltering_system
    R1 = R2 = CartesianRange(())
    dest = src = interp.coefs
    Interpolations._A_ldiv_B_md!(dest, W.A, src, R1, R2, b)
    tmp = updater.tmp
    fill!(tmp, 0)
    tmp[1] = W.Cp[1,1]*src[3]+W.Cp[1,2]*src[end-2]
    tmp[end] = W.Cp[2,1]*src[3]+W.Cp[2, 2]*src[end-2]
    Interpolations._A_ldiv_B_md!(tmp, W.A, tmp, R1, R2)
    AxisAlgorithms.sub!(dest, tmp)
    interp
end

function extract_wf(x, thr_type::Union{PositiveThreshold,NegativeThreshold}, pre, post, align; upsample::Int=1, refractory::Int=post)
    nx = length(x)
    i = pre+1
    nwf = 0
    @inbounds while i <= nx
        if apply_threshold(thr_type, x[i])
            nwf += 1
            i += refractory - 1
        end
        i += 1
    end
    tmp = zeros(eltype(x), upsample*(pre+post)+1)
    updater = FastInterpolationUpdater(zeros(eltype(x), 2*pre+2*post+1))
    interpolated = updater.interp
    wf = zeros(eltype(x), pre+post+1, nwf)
    times = zeros(Int, nwf)
    oldnwf = nwf # Might detect fewer waveforms due to alignment
    nwf = 0
    i = pre+1
    while i <= nx-post
        if apply_threshold(thr_type, x[i])
            nwf += 1

            oldi = -1
            if align
                oldi = i
                for iter = 1:10
                    if upsample != 1
                        interp_start = i-pre-post
                        interp_end = i+pre+post
                        update_interpolation!(updater, x, interp_start:interp_end)
                        for isample = 1:upsample*(pre+post+1)
                            tmp[isample] = interpolated[post+1+(isample-1)/upsample]
                        end
                        pt = interp_start+post+1+(find_alignment_point(thr_type, tmp)-1)/upsample
                        for iwf = 1:pre+post+1
                            wf[iwf, nwf] = interpolated[pt+iwf-1-pre]
                        end
                        i = round(Int, interp_start+pt-1)
                    else
                        copy!(tmp, 1, x, i-pre, pre+post+1)
                        i += find_alignment_point(thr_type, tmp) - 1
                    end

                    if i == oldi || i < max(initiali - pre, lastspike + refractory)
                        i = oldi
                        break
                    end
                end
            end

            if !align || upsample == 1
                for iwf = 1:pre+post+1
                    wf[iwf, nwf] = x[i+iwf-1-pre]
                end
            end

            times[nwf] = i
            i += refractory - 1
            while i <= nx-post && apply_threshold(thr_type, x[i])
                i += 1
            end
        end
        i += 1
    end
    if nwf < oldnwf
        wf = wf[:, 1:nwf]
        resize!(times, nwf)
    end
    wf, times
end

function extract_wf(x, thr_type::DoubleThreshold, pre, post, align; upsample::Int=1, refractory::Int=post)
    nx = length(x)
    i = pre+1
    nwf = 0
    @inbounds while i <= nx
        if apply_threshold(thr_type, x[i])
            nwf += 1
            i += refractory - 1
        end
        i += 1
    end
    updater = FastInterpolationUpdater(zeros(eltype(x), 2*pre+2*post+1))
    wf = zeros(eltype(x), pre+post+1, nwf)
    times = zeros(Int, nwf)
    grad = zeros(eltype(x), 1)
    oldnwf = nwf # Might detect fewer waveforms due to alignment
    nwf = 0
    lastspike = 0
    i = pre+1
    interpolated = updater.interp
    extrema = Tuple{Float64,eltype(x)}[]
    while i <= nx-post
        if apply_threshold(thr_type, x[i])
            nwf += 1

            if align
                initiali = i
                pt = 1/one(upsample)
                # Find the local extremum that corresponds to the highest TEO
                best = -Inf
                secondbest = -Inf
                for iter = 1:10
                    oldi = i
                    interp_start = i-pre-post
                    interp_end = i+pre+post
                    update_interpolation!(updater, x, interp_start:interp_end)
                    empty!(extrema)
                    npos = 0
                    nneg = 0
                    for iext = post+1:1/upsample:2*post+pre+1
                        v = interpolated[iext]
                        !apply_threshold(thr_type, v) && continue
                        gradient!(grad, interpolated, iext-1/(2*upsample))
                        g1 = grad[1]
                        gradient!(grad, interpolated, iext+1/(2*upsample))
                        g2 = grad[1]
                        if sign(g1) != sign(g2)
                            push!(extrema, (iext, v))
                            if v > 0
                                npos += 1
                            else
                                nneg += 1
                            end
                        end
                    end
                    if npos == 1 && nneg > 1
                        # If multiple extrema on one side but only one on another, use the one
                        for (iext, v) in extrema
                            if v > 0
                                pt = iext
                                break
                            end
                        end
                    elseif nneg >= 1
                        # Use largest negative threshold crossing
                        mext = 0.0
                        for (iext, v) in extrema
                            if v < mext
                                pt = iext
                                mext = v
                            end
                        end
                    else
                        # Use the first positive crossing
                        pt = extrema[1][1]
                    end

                    i = round(Int, interp_start+pt-1)
                    if i == oldi || i < max(initiali - pre, lastspike + refractory)
                        i = oldi
                        break
                    end
                end

                for iwf = 1:pre+post+1
                    wf[iwf, nwf] = interpolated[pt+iwf-1-pre]
                end
            else
                for iwf = 1:pre+post+1
                    wf[iwf, nwf] = x[i+iwf-1-pre]
                end
            end

            lastspike = i
            times[nwf] = i
            i += refractory
            while i <= nx-post && apply_threshold(thr_type, x[i])
                i += 1
            end
        end
        i += 1
    end
    if nwf < oldnwf
        wf = wf[:, 1:nwf]
        resize!(times, nwf)
    end
    wf, #=wfteo, extrema,=# times
end

# Find segments of noise suitably far from a spike
function noise_autocovariance(x, thr_type, n, gap)
    acov = zeros(eltype(x), n)
    acovn = zeros(Int, n)
    i = 1
    nx = length(x)
    @inbounds while i <= nx
        j = i
        while j < nx && !apply_threshold(thr_type, x[j+1])
            j += 1
        end
        j -= gap
        for lag = 0:min(n-1, j-i)
            acov[lag+1] += dot(x, i:j-lag, x, i+lag:j)
            acovn[lag+1] += (j-i+1-lag)
        end
        i = j + 2*gap + 1
    end
    acov./acovn
end

whiten(fwf, acov) = sqrtm(full(Toeplitz([reverse(acov[2:end]); acov])))\fwf

end # module


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = Float64.([1 1 1 1 1])

star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))

st_inst = NewStencilInstance(star_stencil, 4, coefs)

radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

ext: jl


text:
module PUCP_Functions

using Distributions

function cdf_Wᵣ(Val::Real,Qₜ::Real,σ :: Real,t::Real)
    if Qₜ <= 0
        error("Qₜ es menor a 0, ello no es posible.")
    α=1/(log(σ + 1))
    Cₜ=(-log(1-t))^(1/α)
    β = Qₜ/Cₜ
    W = Weibull(α,β)
    Wᵣ = cdf(Val,W)
    return Wᵣ
end

function rand_Wᵣ(n::Int,Qₜ::Real,σ::Real,t::Real)
    if t >= 1
        error("t debe estar entre 0 y 1")
    end
    α = 1/(log(σ+1))
    Cₜ=(-log(1-t))^(1/α)
    β = Qₜ/Cₜ
    W = Weibull(α,β)
    random = rand(W,n)
    return random
end

function reg_Wᵣ(DF::DataFrame,ϕᵢ::Int,ϕᵢ₊₁::Int,t::Real)

    function likelihood(Parameters::AbstractVector,ϕᵢ::Int,ϕᵢ₊₁::Int,t::Real)
        β = Parameters[1]
        σ = Parameters[2]
        Cₜ = -log(1-t)
        X = DF ## Agregar 1 para col de diseño y quitar limites
        Qₜ = exp(X*β)

        ## Inicialización de veroslimilitud
        Sₗₖ = 0
        ## Identificar la forma matricial de la verosimilitud
    end

end

end

max_stars_repo_name: jamanrique/MSc-Tesis-PUCP

max_issues_repo_name: jamanrique/MSc-Tesis-PUCP

max_forks_repo_name: jamanrique/MSc-Tesis-PUCP


text:
@testset "n-Queens" begin

    N = 8
    P = 100
    generatePositions(N::Int) = collect(1:N)[randperm(N)]

    # Vector of N cols filled with numbers from 1:N specifying row position
    function nqueens(queens::Vector{Int})
        n = length(queens)
        fitness = 0
        for i=1:(n-1)
            for j=(i+1):n
                k = abs(queens[i] - queens[j])
                if (j-i) == k || k == 0
                    fitness += 1
                end
                # println("$(i),$(queens[i]) <=> $(j),$(queens[j]) : $(fitness)")
            end
        end
        return fitness
    end
    @test nqueens([2,4,1,3]) == 0
    @test nqueens([3,1,2]) == 1

    # Testing: GA solution with various mutations
    for muts in [inversion, insertion, swap2, scramble, shifting]
        result, fitness, cnt = optimize(nqueens, GA(N=N;
            initPopulation = generatePositions,
            populationSize = P,
            selection = sus,
            crossover = pmx,
            mutation = muts))
        println("GA:PMX:$(string(muts))(N=$(N), P=$(P)) => F: $(fitness), C: $(cnt), OBJ: $(result)")
        @test nqueens(result) == 0
        @test cnt < 800 # Test early stopping
    end

    # Testing: ES
    for muts in [inversion, insertion, swap2, scramble, shifting]
        result, fitness, cnt = optimize(nqueens, ES(N=N;
            initPopulation = generatePositions,
            mutation = mutationwrapper(muts),
            μ = 15, ρ = 1, λ = P))
        println("(15+$(P))-ES:$(string(muts)) => F: $(fitness), C: $(cnt), OBJ: $(result)")
        @test nqueens(result) == 0
    end

end


text:
using DifferentialEquations
using DynamicalNetworks
import SimpleWeightedGraphs: add_edge!

# Callback
function condition(u, t, integrator)
    tswitch = integrator.sol.prob.p[2].tswitch
    integrator.t ∈ tswitch
end

function action!(integrator)
    net, params = integrator.sol.prob.p
    if popfirst!(params.bits)
        add_edge!(net.graph, params.srcnode, params.dstnode, params.high)
    else
        add_edge!(net.graph, params.srcnode, params.dstnode, params.low)
    end
end

# Update problem
function prob_func(prob, i, repeat)
    net, params = prob.p
    @info "Solving problem $i / $(length(params.noiserange))"
    net.noisemat[params.ctrlidx] .= params.noiserange[i]
    prob
end

# Output function
function output_func(sol, i)
    getindex.(sol.u, 1), false
end


text:
"""
    AbstractTrajectory

An abstract representation of a state and control trajectory.

| **Required methods**  | **Brief description**  |
|:--------------------- |:---------------------- |
| `getstate(Z, t)` | Get the state at time `t`. |
| `getcontrol(Z, t)` | Get the control at time `t`. |
| `getinitialtime(Z)` | Get the initial time of the trajectory. |
| `getfinaltime(Z)` | Get the final time of the trajectory. |
"""
abstract type AbstractTrajectory end

getstate(Z::AbstractTrajectory, t) = throw(NotImplementedError("getstate not implemented for $(typeof(Z))."))
getcontrol(Z::AbstractTrajectory, t) = throw(NotImplementedError("getcontrol not implemented for $(typeof(Z))."))
getinitialtime(Z::AbstractTrajectory) = throw(NotImplementedError("getinitialtime not implemented for $(typeof(Z))."))
getfinaltime(Z::AbstractTrajectory) = throw(NotImplementedError("getfinaltime not implemented for $(typeof(Z))."))

"""
    SampledTrajectory{Nx,Nu,T,KP} <: AbstractTrajectory

A trajectory represented by a sample of knot points.
A vector of `AbstractKnotPoint`s of type `KP` with state dimension `Nx`,
control dimension `Nu`, and value type `T`.

Supports iteration and indexing.

# Constructors
    SampledTrajectory{Nx,Nu}(n, m; [equal, dt, tf, N])
    SampledTrajectory{Nx,Nu}(X, U; [dt, tf, N=length(X)])

    SampledTrajectory(n, m; [equal, dt, tf, N])
    SampledTrajectory(X, U; [dt, tf, N=length(X)])

where at least 2 of `dt`, `tf`, or `N` must be specified. The length `N` is automatically
inferred when passing in state and control trajectories `X` and `U`, which can either be
vectors of vectors or a 2D matrix whose 2nd dimension is time.
"""
struct SampledTrajectory{n,m,T,KP} <: AbstractTrajectory
    data::Vector{KP}
    times::Vector{T}
    function SampledTrajectory(Z::Vector{<:AbstractKnotPoint{n,m,<:Any,T}}) where {n,m,T}
        times = zeros(T, length(Z))
        for k = 1:length(Z)
            times[k] = time(Z[k])
        end
        new{n,m,T,eltype(Z)}(Z, times)
    end
end

#############################################
# AbstractTrajectory Interface
#############################################
getstate(Z::SampledTrajectory, t) = state(Z.data[getk(Z, t)])
getcontrol(Z::SampledTrajectory, t) = control(Z.data[getk(Z, t)])
getinitialtime(Z::SampledTrajectory) = time(Z.data[1])
getfinaltime(Z::SampledTrajectory) = time(Z.data[end])

#############################################
# Iteration and Indexing
#############################################
@inline Base.iterate(Z::SampledTrajectory) = iterate(Z.data)
@inline Base.iterate(Z::SampledTrajectory, k::Int) = iterate(Z.data, k)
Base.keys(Z::SampledTrajectory) = Base.keys(Z.data)
Base.IteratorSize(Z::SampledTrajectory) = Base.HasLength()
Base.IteratorEltype(Z::SampledTrajectory) = Base.IteratorEltype(Z.data)
@inline Base.eltype(Z::SampledTrajectory) = eltype(Z.data)
@inline Base.length(Z::SampledTrajectory) = length(Z.data)
@inline Base.size(Z::SampledTrajectory) = size(Z.data)
@inline Base.getindex(Z::SampledTrajectory, i) = Z.data[i]
@inline Base.setindex!(Z::SampledTrajectory, v, i) = Z.data[i] = v
@inline Base.firstindex(Z::SampledTrajectory) = 1
@inline Base.lastindex(Z::SampledTrajectory) = lastindex(Z.data)
Base.IndexStyle(::SampledTrajectory) = IndexLinear()

#############################################
# Constructors
#############################################
function SampledTrajectory{Nx,Nu}(n::Int, m::Int; equal=false, kwargs...) where {Nx,Nu}
    times = _gettimeinfo(; kwargs...)
    dt = push!(diff(times), 0.0)
    Z = [KnotPoint{Nx,Nu}(n, m, fill(NaN, n + m), times[k], dt[k]) for k = 1:length(dt)]
    if equal
        Z[end].dt = Inf
    end
    return SampledTrajectory(Z)
end

function SampledTrajectory{Nx,Nu}(X::Matrix, U::Matrix; kwargs...) where {Nx,Nu}
    Xvec = [Vector(x) for x in eachcol(X)]
    Uvec = [Vector(u) for u in eachcol(U)]
    SampledTrajectory{Nx,Nu}(Xvec, Uvec; kwargs...)
end

function SampledTrajectory{Nx,Nu}(X::Vector, U::Vector; kwargs...) where {Nx, Nu}
    N = length(X)
    if !haskey(kwargs, :tf) && !haskey(kwargs, :dt)
        error("Must specify either the time step or the total time.")
    end
    if haskey(kwargs, :N) && kwargs[:N] != N
        error("Specified an N inconsistent with the number of state vectors.")
    end
    times = _gettimeinfo(; N=N, kwargs...)
    dt = diff(times)
    Z = [KnotPoint{Nx,Nu}(length(X[k]), length(U[k]), [X[k]; U[k]], times[k], dt[k]) for k = 1:N-1]
    if length(U) == length(X)
        push!(Z, KnotPoint{Nx,Nu}(length(X[end]), length(U[end]), [X[end]; U[end]], times[N], Inf))
    else
        push!(Z, KnotPoint{Nx,Nu}(length(X[end]), length(U[end]), [X[end]; U[end]*0], times[N], 0.0))
    end
    return SampledTrajectory(Z)
end

function SampledTrajectory(X::Vector{<:StaticVector{Nx}}, U::Vector{<:StaticVector{Nu}};
              kwargs...) where {Nx,Nu}
    SampledTrajectory{Nx,Nu}(X, U; kwargs...)
end

SampledTrajectory(Z::SampledTrajectory) = Z

SampledTrajectory(args...; kwargs...) = SampledTrajectory{Any,Any}(args...; kwargs...)

#############################################
# Getters
#############################################
vectype(Z::SampledTrajectory{<:Any,<:Any,<:Any,KP}) where KP = vectype(KP)
has_terminal_control(Z::SampledTrajectory) = !RobotDynamics.is_terminal(Z[end])
state_dim(Z::SampledTrajectory{n}) where {n} = n
control_dim(Z::SampledTrajectory{<:Any,m}) where {m} = m
state_dim(Z::SampledTrajectory, k::Integer) = state_dim(Z[k])
control_dim(Z::SampledTrajectory, k::Integer) = control_dim(Z[k])
dims(Z::SampledTrajectory) = (state_dim.(Z), control_dim.(Z), length(Z))
getk(Z::SampledTrajectory, t::Real) = searchsortedfirst(Z.times, t)

"""
    num_vars(Z)

Total number of states and controls in a trajectory `Z`.
"""
function num_vars(Z::SampledTrajectory)
    mapreduce(+, Z) do z
        state_dim(z) + (is_terminal(z) ? 0 : control_dim(z))
    end
end

function num_vars(n::Int, m::Int, N::Int, equal::Bool=false)
    Nu = equal ? N : N-1
    return N*n + Nu*m
end

"""
    eachcontrol(Z)

Get the range of indices for valid controls.
"""
eachcontrol(Z::SampledTrajectory) = has_terminal_control(Z) ? Base.OneTo(length(Z)) : Base.OneTo(length(Z)-1)

"""
    states(Z)
    states(Z, i::Integer)
    states(Z, inds)

Get a list of all the state vectors for the trajectory `Z`. Passing
an integer extracts a vector of the `i`th state. Passing a vector of
integers provides a list of `N`-dimensional vectors, containing the
time history for each state index in the vector.
"""
@inline states(Z::SampledTrajectory) = state.(Z)
states(Z::SampledTrajectory, inds::AbstractVector{<:Integer}) = [states(Z, i) for i in inds]
states(Z::SampledTrajectory, ind::Integer) = [state(z)[ind] for z in Z]

"""
    controls(Z)
    controls(Z, i::Integer)
    controls(Z, inds)

Get a list of all the control vectors for the trajectory `Z`. Passing
an integer extracts a vector of the `i`th control. Passing a vector of
integers provides a list of `N`-dimensional vectors, containing the
time history for each control index in the vector.
"""
function controls(Z::SampledTrajectory)
    return [control(Z[k]) for k in eachcontrol(Z) ]
end
controls(Z::SampledTrajectory, inds::AbstractVector{<:Integer}) = [controls(Z, i) for i in inds]
controls(Z::SampledTrajectory, ind::Integer) = [control(Z[k])[ind] for k in eachcontrol(Z)]

"""
    gettimes(Z)

Get a vector of times for the entire trajectory.
"""
function gettimes(Z::SampledTrajectory)
    [time(z) for z in Z]
end

"""
    getdata(Z)

Get a list of the concatenated state and control vectors.
"""
getdata(Z::SampledTrajectory) = get_data.(Z)

#############################################
# Setters
#############################################
"""
    setstates!(Z, X)

Set the states of a trajectory `Z`, where `X` can be a vector of vectors or a matrix
of size `(n,N)`.
"""
function setstates!(Z::SampledTrajectory, X)
    for k in eachindex(Z)
        setstate!(Z[k], X[k])
    end
end

function setstates!(Z::SampledTrajectory, X::AbstractMatrix)
    for k in eachindex(Z)
        setstate!(Z[k], view(X, :, k))
    end
end

"""
    setcontrols!(Z, U)

Set the controls of a trajectory `Z`, where `U` can be a vector of vectors or a matrix
of size `(m,N)` or a single vector of size `(m,)`, which will be copied to all the
time steps.
"""
function setcontrols!(Z::SampledTrajectory, U)
    for k in 1:length(Z)-1
        setcontrol!(Z[k], U[k])
    end
end

function setcontrols!(Z::SampledTrajectory, U::AbstractMatrix)
    for k in 1:length(Z)-1
        setcontrol!(Z[k], view(U, :, k))
    end
end

function setcontrols!(Z::SampledTrajectory, u::AbstractVector{<:Real})
    for k in 1:length(Z)-1
        setcontrol!(Z[k], u)
    end
end

"""
    setdata!(Z, V)

Set the concatenated state and control vector for each knot point in the trajectory `Z`.
`V` may be either a vector of vectors or a 2D matrix.
"""
function setdata!(Z::SampledTrajectory, V)
    for k in 1:length(Z)-1
        setdata!(Z[k], V[k])
    end
end

function setdata!(Z::SampledTrajectory, V::AbstractMatrix)
    for k in 1:length(Z)-1
        setdata!(Z[k], view(V, :, k))
    end
end

"""
    settimes!(Z, ts)

Set the times for the entire trajectory. The time steps are automatically updated.
"""
function settimes!(Z::SampledTrajectory, ts)
    for k in eachindex(ts)
        Z[k].t = ts[k]
        k < length(ts) && (Z[k].dt = ts[k+1] - ts[k])
    end
end

"""
    set_dt!(Z, dt)

Set a constant time step for the entire trajectory.
"""
function set_dt!(Z::SampledTrajectory, dt::Real)
    t = Z[1].t
    for z in Z
        z.t = t
        if !is_terminal(z)
            z.dt = dt
            t += dt
        end
    end
    return t
end

"""
    setinitialtime!(Z, t0)

Set the initial time of the trajectory, shifting all of the times by the required amount.
"""
function setinitialtime!(Z::SampledTrajectory, t0)
    t0_prev = time(Z[1])
    Δt = t0 - t0_prev
    for z in Z
        t = time(z)
        settime!(z, t + Δt)
    end
    return Z
end

#############################################
# Copying and Comparison
#############################################
function Base.copy(Z::SampledTrajectory{Nx,Nu}) where {Nx,Nu}
    SampledTrajectory([KnotPoint{Nx,Nu}(z.n, z.m, copy(z.z), z.t, z.dt) for z in Z])
end

function Base.isapprox(Z1::SampledTrajectory, Z2::SampledTrajectory)
    all(zs->zs[1] ≈ zs[2], zip(Z1,Z2))
end

function Base.copyto!(Z::SampledTrajectory, Z0::SampledTrajectory)
    @assert length(Z) == length(Z0)
    V = vectype(Z)
    N = length(Z)
    for k = 1:N-1
        if V <: StaticVector
            Z[k].z = Z0[k].z
        else
            copyto!(Z[k].z, Z0[k].z)
        end
        Z[k].t = Z0[k].t
        Z[k].dt = Z0[k].dt
    end
    if is_terminal(Z[end])
        setstate!(Z[end], state(Z0[end]))
    else
        if V <: StaticVector
            Z[N].z = Z0[N].z
        else
            copyto!(Z[N].z, Z0[N].z)
        end
    end
    Z[N].t = Z0[N].t
    Z[N].dt = Z0[N].dt
    Z
end

#############################################
# Misc
#############################################
function _gettimeinfo(;t0=0.0, tf=NaN, dt=NaN, N=0)
    if dt isa Vector
        return _gettimeinfo(dt, t0=t0, tf=tf, N=N)
    end
    Δt = zero(t0)
    if isnan(tf) + isnan(dt) + (N == 0) > 1
        error("Must specify at least two of the following: dt, tf, N")
    end

    # Given tf, dt
    if (N==0) && !isnan(tf) && !isnan(dt)
        N = round(Int, tf / dt + 1)
        Δt = tf - t0
        dt = Δt / (N - 1)
    end

    # Given N, tf
    if !isnan(tf) && isnan(dt)
        dt = tf / (N-1)
    end

    # Given N, dt
    if !isnan(dt) && isnan(tf)
        tf = dt * (N - 1)
    end
    Δt = tf - t0

    # Check consistency
    @assert tf ≈ dt * (N-1) "Inconsistent time step and final time."
    @assert dt ≈ Δt / (N - 1) "Inconsistent time step and trajectory length."
    return range(t0, tf, length=N)
end

function _gettimeinfo(dt::AbstractVector; t0=0, tf=NaN, N=0)
    @assert all(x->x>zero(x), dt) "All time steps must be positive"
    if (N == 0)
        N = length(dt) + 1
    end
    if isnan(tf)
        tf = sum(dt) + t0
    end
    Δt = tf - t0
    @assert length(dt) == N-1 "Time step vector must have length of N-1, got $(length(dt)) with N=$N."
    @assert Δt ≈ sum(dt) "Elapsed time ($Δt) must be equal to the sum of the time steps (got $(sum(dt)))."
    return [t0; cumsum(dt)]
end

function shift_fill!(Z::SampledTrajectory, n=1)
    N = length(Z)
    isterm = is_terminal(Z[end])
    for k = 1+n:N
        Z[k-n] = copy(Z[k])
    end
    xf = state(Z[N-n])
    uf = control(Z[N-n])
    dt = Z[N-n-1].dt
    for k = N-n:N
        setstate!(Z[k], xf)
        Z[k].t = Z[k-1].t + dt
        if k == N && is_terminal(Z[k])
            Z[k].dt = 0
        else
            setcontrol!(Z[k], uf)
            Z[k].dt = dt
        end
    end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~ FUNCTIONS ON TRAJECTORIES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

# state_diff_jacobian!(model::AbstractModel, G, Z::AbstractTrajectory) =
#     state_diff_jacobian!(statevectortype(model), model, G, Z)
# state_diff_jacobian!(::EuclideanState, model::AbstractModel, G, Z::AbstractTrajectory) = nothing

# function state_diff_jacobian!(::RotationState, model::AbstractModel, G, Z::AbstractTrajectory)
#     for k in eachindex(Z)
#         G[k] .= 0
#         state_diff_jacobian!(RotationState(), model, G[k], Z[k])
#     end
# end

function rollout!(sig::FunctionSignature, model::DiscreteDynamics, Z::AbstractTrajectory, x0=state(Z[1]))
    setstate!(Z[1], x0)
    for k = 2:length(Z)
        RobotDynamics.propagate_dynamics!(sig, model, Z[k], Z[k-1])
    end
end


text:
MAIN_KEYS = ["query_history", "subspace_grids", "subspaces", "data", "labels", "params"]
PARAM_KEYS = ["C", "gamma", "classifier", "query_strategy"]

function check_keys(r)
    missing_keys = []
    for k ∈ MAIN_KEYS
        if k ∉ keys(r)
            push!(missing_keys, k)
        end
    end
    if "params" ∈ keys(r)
        for k ∈ PARAM_KEYS
            if k ∉ keys(r["params"])
                push!(missing_keys, k)
            end
        end
    end
    if !isempty(missing_keys)
        throw(ArgumentError("Missing keys in request: $(missing_keys)."))
    end
    return nothing
end

function convert_value!(r, name, target_type, func=identity)
    try
        r[name] = func(convert(target_type, r[name]))
    catch e
        throw(ArgumentError("Failed parsing '$(name)': $(r[name])."))
    end
end

function convert_value!(r, name, target_type::Type{Vector{Symbol}}, func=identity)
    try
        r[name] = func([Symbol.(x) for x in r[name]])
    catch e
        throw(ArgumentError("Failed parsing '$(name)': $(r[name])."))
    end
end

function check_and_parse_values!(r)
    convert_value!(r, "data", Vector{Vector{Float64}}, x -> hcat(x...))
    convert_value!(r, "labels", Vector{Symbol})
    convert_value!(r, "params", Dict{String, Any})
    convert_value!(r, "query_history", Vector{Vector{Int}})
    convert_value!(r, "subspaces", Vector{Vector{Int}})
    convert_value!(r, "subspace_grids", Vector{Vector{Vector{Float64}}}, x -> [hcat(s...) for s in x])
    return nothing
end

function check_value_sizes(r)
    size(r["data"], 2) == length(r["labels"]) || throw(ArgumentError("Data size does not match number of labels."))
    length(r["subspaces"]) == length(r["subspace_grids"]) || throw(ArgumentError("Number of subspaces does not match requested number of grids."))
    return nothing
end

function parse_parameters!(r)
    check_keys(r)
    check_and_parse_values!(r)
    check_value_sizes(r)
    return nothing
end


text:
# A module that returns the zeroth and first persistent diagram
# of a given protein. It uses the ripser algorithm which utilizes the
# vietoris-rips complex filtration

#__precompile__()

module ProteinPersistent

export returndiagram, valuesvariable, coordpdb

using PyCall
using DataFrames

const ripser = PyNULL()
const pdb = PyNULL()
const re = PyNULL()
const warnings = PyNULL()

function __init__()
    copy!(warnings, pyimport("warnings"))
    warnings[:filterwarnings]("ignore")
    copy!(ripser, pyimport("ripser"))
    copy!(pdb, pyimport("Bio.PDB"))
    copy!(re, pyimport("re"))
end

# define the functions to load the pdb file as a dataframe with
# coordinates, atoms and van der waals radii

# return the shortname of the atom (e.g. CA returns C)
function atom_shortname(s)
    pos = re[:search]("[A-Z]", s)
    pos = pos[:start]()
    return s[pos+1]
end

# return van der Waals radius associated with 'atom_name'
function vanderWaalsRadius(atom_name)
    d = Dict('H' => 1.2, 'C' => 1.7, 'N' => 1.55, 'O' => 1.52,
             'S' => 1.8, 'P' => 1.8, 'M' => 1.73)
    return d[atom_name]
end

# return the coordinates, atom sequence and van der wall radii as a Dataframe
function coordpdb(filename; kwargs...)

    atoms = ["CA", "C", "N", "H"]

    for (p,v) in kwargs
        if p == :atoms
            atoms = v
        end
    end

    parser = pdb[:PDBParser]()
    structure = parser[:get_structure]("", filename)

    df = DataFrame(Any, 0, 5)

    for atom in structure[:get_atoms]()
        coordinates = atom[:get_coord]()
        name_atom = atom[:get_name]()
        if name_atom in atoms
            atom_short = atom_shortname(atom[:get_name]())
            radius = vanderWaalsRadius(atom_short)
            list_to_df = [coordinates[1], coordinates[2],
                          coordinates[3], name_atom, radius]
            push!(df, list_to_df)
        end
    end
    name = zip([:x1, :x2, :x3, :x4, :x5], [:coord1, :coord2, :coord3,
                                           :atom_name, :radius])
    rename!(df, f => t for (f,t) = name)
    return df
end

# returns the zeroth and first persistent diagram of the given pdb file
function returndiagram(pdbfile::String; kwargs...)

    maxdim = 1
    for (p,v) in kwargs
        if p == :maxdim
            maxdim = v
        end
    end

    # get the coordinates of the protein
    coordinates = coordpdb(pdbfile)
    # convert to an array
    coordinates = Matrix{Float64}(coordinates[[:coord1, :coord2, :coord3]])
    # calculate the 0th and 1st persistent diagram using Vietoris Rips
    diagrams = ripser[:ripser](coordinates, maxdim=maxdim)
    return diagrams["dgms"]
end

end


text:
__precompile__()

module XDiff

export
    _xdiff,
    xdiff,
    @scalardiff,
    @tensordiff,
    @specialdiff,
    # reexport from Espresso
    @get_or_create,
    sum_1,
    sum_2,
    squeeze_sum,
    squeeze_sum_1,
    squeeze_sum_2,
    # reexport from Einsum
    @einsum

include("core.jl")

end # module


text:
########################################################################
#                            GRAD RESULT                               #
########################################################################

function field_paths(tape::Tape)
    paths = Dict()
    for op in reverse(tape.ops)
        _op = op
        path = []
        while _op isa Call && _op.fn in (Base.getproperty, __getfield__)
            field_name = tape[_op.args[2]].val
            push!(path, field_name)
            _op_id = _op.args[1]
            _op = tape[_op_id]
        end
        if !isempty(path)
            struct_id = _op.id
            if !haskey(paths, struct_id)
                paths[struct_id] = Dict()
            end
            paths[struct_id][(reverse(path)...,)] = op.id
        end
    end
    return paths
end

struct GradResult
    tape::Tape
    gvars::Dict{Int, Any}  # gradient vars: argid -> gradient var
end

function GradResult(tape::Tape)
    tape.fieldpaths = field_paths(tape)
    gvars = Dict{Int,Any}()
    # struct fields
    for (argid, dct) in tape.fieldpaths
        tape[argid] isa Input || continue  # skip non-input struct fields
        gvars[argid] = Dict(field_path => tape.derivs[var_id]
                            for (field_path, var_id) in dct
                            if haskey(tape.derivs, var_id))  # not all fields may have derivatives
    end
    # other arguments
    struct_arg_ids = Set(keys(tape.fieldpaths))
    for op in tape
        if (op isa Input                    # if operation is an Input
            && !in(op.id, struct_arg_ids)   # and not a struct
            && haskey(tape.derivs, op.id))  # and there's derivative for this var
            gvars[op.id] = tape.derivs[op.id]
        end
    end
    return GradResult(tape, gvars)
end

Base.show(io::IO, g::GradResult) = print(io, "GradResult($(length(g.gvars)))")

function Base.getindex(g::GradResult, argid::Int)
    tape = g.tape
    gvar = g.gvars[argid]
    if isa(gvar, Dict)
        return Dict(f => tape[id].val for (f, id) in gvar)
    else
        return tape[gvar].val
    end
end

Base.length(g::GradResult) = length(g.gvars)
Base.iterate(g::GradResult) = length(g) > 0 ? (g[1], 2) : nothing
Base.iterate(g::GradResult, state) = length(g) >= state ? (g[state], state + 1) : nothing

########################################################################
#                              GRAD                                    #
########################################################################

const DEBUG_STATE = Any[]

getderiv(tape::Tape, id::Int) = haskey(tape.derivs, id) ? tape[tape.derivs[id]] : nothing
getderiv(tape::Tape, op::AbstractOp) = getderiv(tape, op.id)
setderiv!(tape::Tape, op_id::Int, grad_op_id::Int) = (tape.derivs[op_id] = grad_op_id)
setderiv!(tape::Tape, op::AbstractOp, grad_op::AbstractOp) = (tape.derivs[op.id] = grad_op.id)

Espresso.to_expr(tape::Tape, op::Call) = begin
    Expr(:call, op.fn, [Symbol("%$i") for i in op.args]...)
end

to_unbroadcast_expr(tape::Tape, op::Call) =
    Expr(:call, tape[op.args[1]].val, [Symbol("%$i") for i in op.args[2:end]]...)

function deriv!(tape::Tape, op::AbstractOp, i::Int, dy::AbstractOp)
    ex = to_expr(tape, op)
    dep_types = [tape[arg].typ for arg in op.args]
    dex = deriv_expr(ex, dep_types, i)
    st = Dict(Symbol("%$i") => i for i in op.args)
    st[:dy] = dy.id
    st[:y] = op.id
    ret_id = record_expr!(tape, dex; st=st)
    return tape[ret_id]
end

function deriv_broadcast!(tape::Tape, op::AbstractOp, i::Int, dy::AbstractOp)
    ex = to_unbroadcast_expr(tape, op)
    dep_eltypes = [eltype(tape[arg].typ) for arg in op.args[2:end]]
    dex = deriv_expr(ex, dep_eltypes, i-1)
    st = Dict(Symbol("%$id") => id for id in op.args)
    st[:dy] = dy.id
    st[:y] = op.id
    ret_id = record_expr!(tape, dex; st=st, bcast=true)
    return tape[ret_id]
end

function set_or_add_deriv!(tape::Tape, x::AbstractOp, dx::AbstractOp)
    if !haskey(tape.derivs, x.id)
        setderiv!(tape, x, dx)
    else
        old_dx = getderiv(tape, x)
        val = dx.val .+ old_dx.val
        dot_add_id = record!(tape, Constant, +)
        new_dx_id = record!(tape, Call, val, broadcast, [dot_add_id, dx.id, old_dx.id])
        new_dx = tape[new_dx_id]
        setderiv!(tape, x, new_dx)
    end
end

function step_back!(tape::Tape, op::Union{Call}, i::Int)
    y = op
    x = tape[op.args[i]]
    dy = getderiv(tape, y)
    dy == nothing && return   # op is not part of computation graph, e.g. range in loop
    # we handle broadcasting for a few built-in functions like normal derivatives
    # all other cases are handled by a generic bcast mechanism
    use_bcast_rules = (op.fn == broadcast) && !in(tape[op.args[1]].val, Set([+, *]))
    dx = try
        use_bcast_rules ? deriv_broadcast!(tape, op, i, dy) : deriv!(tape, op, i, dy)
    catch
        @error("Failed to find a derivative for $op at position $i, " *
               "current state of backpropagation saved to Yota.DEBUG_STATE")
        push!(DEBUG_STATE, (tape, op, i))
        rethrow()
    end
    set_or_add_deriv!(tape, x, dx)
end

"""
Backpropagate through the tape, record derivatives as new operations
"""
function back!(tape::Tape)
    # z - final variable (usually a loss)
    # y - resulting variable of current op
    # x - dependencies of y
    # dy - derivative of z w.r.t. y
    z = tape[tape.resultid]
    # using one() of type of the result for seed to keep type stability
    dy_id = record!(tape, Constant, one(tape[tape.resultid].val))
    dy = tape[dy_id]
    # set initial derivative value
    tape.derivs[z.id] = dy.id
    for op in reverse(tape.ops[1:end-1])
        if op isa Call
            if op.fn == Base.getproperty
                # since first argument of getproperty is a struct,
                # we treat it in a special way: instead of inventing some kind
                # of "struct derivative", we look for the call to __new__()
                # that created this struct; if found, we add current derivative to
                # the variable used as the field when instantiating the struct
                if haskey(tape.derivs, op.id)
                    dy_id = tape.derivs[op.id]
                    dy_id != nothing || continue
                    x = find_field_source_var(tape, op)
                    if x != nothing
                        set_or_add_deriv!(tape, x, tape[dy_id])
                        # tape.derivs[x_id] = dy_id
                    end
                end
            elseif op.fn == __getfield__
                # unstructuring of tuples is lowered into pretty weird code sequence
                # ending with __getfield__; similar to getproperty(), we find source var
                # for the corresponding tuple argument and backprop to it
                if haskey(tape.derivs, op.id)
                    dy_id = tape.derivs[op.id]
                    dy_id != nothing || continue
                    x = find_tuple_field_source_var(tape, op)
                    if x != nothing
                        set_or_add_deriv!(tape, x, tape[dy_id])
                        # tape.derivs[x_id] = dy_id
                    end
                end
            elseif op.fn in (__new__, Base.indexed_iterate)
                # we take care of these operations somewhere else
                # so just skip them here
            else
                # ordinary function call
                for i=1:length(op.args)
                    # backpropagate only non-constant vars
                    # note that it also prevents backprop on 1st param of broadcast
                    arg_op = tape[op.args[i]]
                    if !isa(arg_op, Constant) && !dont_diff(tape, op, i)
                        # println(op, " ", i)
                        step_back!(tape, op, i)
                    end
                end
            end
        end
    end
end

"""
For each input that has a derivative on this tape check if the derivative
has the same size as the input.
"""
function check_deriv_sizes(tape::Tape)
    for (var_id, grad_var_id) in tape.derivs
        if !isstruct(tape[var_id].val)
            var_size = size(tape[var_id].val)
            grad_var_size = size(tape[grad_var_id].val)
            if  var_size != grad_var_size
                @warn "Gradient %$grad_var_id has size $grad_var_size, " *
                    "but original variable %$var_id has size $var_size"
            end
        end
    end
end

"""
Calculate and record to the tape gradients of `tape[tape.resultid]` w.r.t. `Input` nodes
"""
function _grad(tape::Tape)
    # apply preprocessing transformations
    tape = preprocess(tape)
    # backpropagate gradients
    back!(tape)
    # consistency check
    check_deriv_sizes(tape)
    # apply postprocessing transformations
    tape = postprocess(tape)
    return tape
end

function _grad(f::Function, args...)
    val, tape = trace(f, args...)
    # calculate gradients
    tape = _grad(tape)
    # construct GradResult object that wraps tape and provides accessors for computed derivatives
    return val, GradResult(tape)
end

const DYNAMIC_GRAD_CACHE = Dict{Any, Tape}()

function _dynamic_grad(f::Function, args...)
    val, tape = trace(f, args...)
    if haskey(DYNAMIC_GRAD_CACHE, tape)
        # take already processed tape from the cache and just play it
        key = tape
        tape = DYNAMIC_GRAD_CACHE[key]
        # play to propagate both - value (should be unchanged) and gradients
        val = play!(tape)
        return val, GradResult(tape)
    else
        # copy tape just after tracing to use as key in cache later
        key = deepcopy(tape)
        # calculate gradients
        tape = _grad(tape)
        # save to cache
        DYNAMIC_GRAD_CACHE[key] = tape
        return val, GradResult(tape)
    end
end

const GRAD_CACHE = Dict{Any, Tape}()

"""
Find gradient of `f` w.r.t. its arguments.
Example:

    val, g = grad(sum, rand(3))

where:
  - val is the value of `f` at this point
  - g::GradResult is a collection of gradients

GradResult is indexed by argument index and contains gradients
in a format most suitable for that argument, namely:

  - for arrays: arrays of the same type and size
  - for reals: reals
  - for mutable structs: dictionary of {(:field, :path) => value} pairs.

All gradients can be applied to original variables using `update!()` function.
"""
function grad(f::Function, args...; dynamic=false)
    # key consists of function type and type of argument (for structs) or its size
    cache_key = (f, ([isstruct(arg) ? typeof(arg) : size(arg) for arg in args]...,))
    if dynamic
        return _dynamic_grad(f, args...)
    elseif haskey(GRAD_CACHE, cache_key)
        tape = GRAD_CACHE[cache_key]
        val = play!(tape, args...)
        return val, GradResult(tape)
    else
        val, g = _grad(f, args...)
        compile!(g.tape)
        GRAD_CACHE[cache_key] = g.tape
        return val, g
    end
end

function simplegrad(f, args...)
    val, g = _grad(f, args...)
    return compile(g.tape, bind=false, ret_grad=true)
end


text:
using Distributions, Interpolations

function simulate(
		model::CartpoleFriction,
		policy,
		K,
		z_nom, u_nom,
		Q, R,
		T_sim, Δt,
		z0, w;
        _norm = 2,
        ul = -Inf * ones(length(u_nom[1])),
        uu = Inf * ones(length(u_nom[1])),
        friction = false,
        μ = 0.1)

    T = length(K) + 1
    times = [(t-1) * Δt for t = 1:T-1]
    tf = Δt * T
    t_sim = range(0, stop = tf, length = T_sim)
    t_ctrl = range(0, stop = tf, length = T)
    dt_sim = tf / (T_sim - 1)

    p = 1:policy.output

    A_state = hcat(z_nom...)
    A_ctrl = hcat(u_nom...)

    z_rollout = [z0]
    u_rollout = []
    J = 0.0
    Jx = 0.0
    Ju = 0.0
    for tt = 1:T_sim-1
        t = t_sim[tt]
        k = searchsortedlast(times, t)

        z_cubic = zeros(model.n)
        for i = 1:model.n
            interp_cubic = CubicSplineInterpolation(t_ctrl, A_state[i,:])
            z_cubic[i] = interp_cubic(t)
        end

        z = z_rollout[end] + dt_sim * w[:,tt]
        u = eval_policy(policy, K[k], z, z_cubic, u_nom[k])

        # clip controls
        u = max.(u, ul[p])
        u = min.(u, uu[p])

        if friction
            _u = [u[1] - μ * sign(z_cubic[3]) * model.g * (model.mp + model.mc);
				  0.0;
				  0.0]
        else
            _u = u[1]
        end

        push!(z_rollout, rk3(model, z, _u, zeros(model.d), dt_sim))
        push!(u_rollout, u)

        if _norm == 2
            J += (z_rollout[end] - z_cubic)' * Q[k + 1] * (z_rollout[end] - z_cubic)

			J += (u_rollout[end] - u_nom[k][p])' * R[k][p,
                p] * (u_rollout[end] - u_nom[k][p])

			Jx += (z_rollout[end] - z_cubic)' * Q[k + 1] * (z_rollout[end] - z_cubic)

			Ju += (u_rollout[end] - u_nom[k][p])' * R[k][p,
                p] * (u_rollout[end] - u_nom[k][p])
        else
            J += norm(sqrt(Q[k + 1]) * (z_rollout[end] - z_cubic), _norm)
            J += norm(sqrt(R[k][p, p]) * (u - u_nom[k][p]), _norm)
            Jx += norm(sqrt(Q[k + 1]) * (z_rollout[end] - z_cubic), _norm)
            Ju += norm(sqrt(R[k][p, p]) * (u - u_nom[k][p]), _norm)
        end
    end
    return z_rollout, u_rollout, J / (T_sim - 1), Jx / (T_sim - 1), Ju / (T_sim - 1)
	nothing
end

# Nominal trajectories
X̄_nominal, Ū_nominal = unpack(Z̄_nominal, prob_nominal)
X̄_friction, Ū_friction = unpack(Z̄_friction, prob_friction)
X̄_dpo, Ū_dpo = unpack(Z[prob_dpo.prob.idx.nom], prob_dpo.prob.prob.nom)

# Policies
K_nominal = tvlqr(model, X̄_nominal, [Ū_nominal[t][1:1] for t = 1:T-1],
 	Q, [R[t][1:1, 1:1] for t = 1:T-1], h)
K_friction = tvlqr(model, X̄_friction, [Ū_friction[t][1:1] for t = 1:T-1],
 	Q, [R[t][1:1, 1:1] for t = 1:T-1], h)
θ = [reshape(Z[prob_dpo.prob.idx.policy[prob_dpo.prob.idx.θ[t]]],
		1, model.n) for t = 1:T-1]

# Simulation
T_sim = 10T
Δt = h
dt_sim = tf / (T_sim - 1)

W = Distributions.MvNormal(zeros(model.n), Diagonal(1.0e-5 * ones(model.n)))
w = rand(W, T_sim)

W0 = Distributions.MvNormal(zeros(model.n), Diagonal(1.0e-5 * ones(model.n)))
w0 = rand(W0, 1)

model_sim = model_friction
μ_sim = 0.1

t_sim = range(0, stop = tf, length = T_sim)

z_lqr, u_lqr, J_lqr, Jx_lqr, Ju_lqr = simulate(model_sim,
	policy,
	K_nominal,
    X̄_nominal, Ū_nominal,
	Q, R,
	T_sim, h,
	x1 + vec(w0),
	w,
	ul = ul_friction, uu = uu_friction,
    friction = true,
    μ = μ0)

z_lqr_fr, u_lqr_fr, J_lqr_fr, Jx_lqr_fr, Ju_lqr_fr = simulate(model_sim,
	policy,
	K_friction,
    X̄_friction, Ū_friction,
	Q, R,
	T_sim, h,
	x1 + vec(w0),
	w,
	ul = ul_friction, uu = uu_friction,
    friction = true,
    μ = μ0)

z_dpo, u_dpo, J_dpo, Jx_dpo, Ju_dpo = simulate(model_sim,
	policy,
	θ,
    X̄_dpo, Ū_dpo,
	Q, R,
	T_sim, h,
	x1 + vec(w0),
	w,
	ul = ul_friction, uu = uu_friction,
    friction = true,
    μ = μ0)

# objective value
J_lqr
J_lqr_fr
J_dpo

# state tracking
Jx_lqr
Jx_lqr_fr
Jx_dpo

# control tracking
Ju_lqr
Ju_lqr_fr
Ju_dpo

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Pkg
Pkg.activate(".")
using HPCManager

init_cluster()

@everywhere begin
	using Pkg
	Pkg.activate(".")
end


text:
using BenchmarkTools
using Dates: Period, Day, Hour, DateFormat, @dateformat_str
using TimeZones
using TimeZones.TZData: parse_components, build
using Test: GenericString

const SUITE = BenchmarkGroup()

function gen_range(num_units::Period, tz::TimeZone)
    zdt = ZonedDateTime(2000, tz)
    return StepRange(zdt, oneunit(num_units), zdt + num_units)
end

# Ensure that when comparing benchmarks
# That the compiled TZData is compatible with this version
build("2016j") # version consistent with tests

include("tzdata.jl")
include("interpret.jl")
include("zoneddatetime.jl")
include("arithmetic.jl")
include("parse.jl")


text:
@testset "Normal form f4 finite fields" begin

    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x", "y", "z"])

    Gs = [
        [x], [x, y], [x, y, z]
    ]

    @test Groebner.normalform(Gs[1], x) == 0
    @test Groebner.normalform(Gs[1], y) == y
    @test Groebner.normalform(Gs[3], x + y^2) == 0
    @test Groebner.normalform(Gs[2], 5x^2 - y^3 - 1) == -1
    @test Groebner.normalform(Gs[2], x + y + 3z) == 3z

    G = [
        x^2 + y*x + 1,
        y^2 - y
    ]

    @test Groebner.normalform(G, x + y + z) == x + y + z
    @test Groebner.normalform(G, x^2 + y^2) == y - y*x - 1
    @test Groebner.normalform(G, y^3) == Groebner.normalform(G, y^4) == y

    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:degrevlex)
    G = [
        x^2 + y,
        y^2 + x
    ]
    @test Groebner.normalform(G, x^2 + y^2) == -x - y

end


text:
using Revise, ApproxOperator

# elements, nodes = importmsh("./msh/bar.msh")
elements, nodes = importmsh("./msh/patchtest.msh")
nₚ = length(nodes[:x])

# type = (Node,:Quadratic2D,:□,:CubicSpline)
# type = (SNode,:Quadratic2D,:□,:QuinticSpline)
type = (SNode,:Cubic2D,:□,:QuinticSpline)
# s = 2.5/20*ones(nₚ)
s = 3.5/20*ones(nₚ)
push!(nodes,:s₁=>s,:s₂=>s,:s₃=>s)

sp = RegularGrid(nodes[:x],nodes[:y],nodes[:z],n = 3,γ = 5)
elements["Ω"] = ReproducingKernel{type...,:Tri3}(elements["Ω"],sp)
set𝓖!(elements["Ω"],:TriRK6)
set_memory_𝗠!(elements["Ω"],:∂1,:∂x,:∂y,:∂z,:∇̃)
set_memory_𝝭!(elements["Ω"],:∂1,:∂x,:∂y,:∂z,:∂x²,:∂x∂y,:∂y²,:∂x∂z,:∂y∂z,:∂z²)
# set𝓖!(elements["Ω"],:TriGI13,:∂1,:∂x,:∂y,:∂z,:∂x²,:∂x∂y,:∂y²,:∂z²,:∂x∂z,:∂y∂z,:∂x³,:∂x²∂y,:∂x∂y²,:∂y³)

# set𝝭!(elements["Ω"])
set∇𝝭!(elements["Ω"])
# set∇̃𝝭!(elements["Ω"])
set∇̃²𝝭!(elements["Ω"])
# set∇²𝝭!(elements["Ω"])
# set∇³𝝭!(elements["Ω"])
# set∇̃𝝭!(elements["Ω"])
# set∇̃²𝝭!(elements["Ω"])
# f = checkConsistency(elements["Ω"])
# f = checkConsistency(elements["Ω"],ApproxOperator.get∇𝝭,ApproxOperator.get∇𝒑)
f = checkConsistency(elements["Ω"],ApproxOperator.get∇²𝝭,ApproxOperator.get∇²𝒑)
# f = checkConsistency(elements["Ω"],ApproxOperator.get∇³𝝭,ApproxOperator.get∇³𝒑)

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
using LeastSquaresOptim, Statistics, Test

# StRD tests (in progress)
# http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml

function misra1a()
    name = "misra1a"
    data = [
    10.07E0      77.6E0;
    14.73E0     114.9E0;
    17.94E0     141.1E0;
    23.93E0     190.8E0;
    29.61E0     239.9E0;
    35.18E0     289.0E0;
    40.02E0     332.8E0;
    44.82E0     378.4E0;
    50.76E0     434.8E0;
    55.05E0     477.3E0;
    61.01E0     536.8E0;
    66.40E0     593.1E0;
    75.47E0     689.1E0;
    81.78E0     760.0E0;
    ]
    parameters =[
    5.0000000000E+02        2.5000000000E+02        2.3894212918E+02;
    1.0000000000E-04        5.0000000000E-04        5.5015643181E-04;
    ]
    f(x, beta) = beta[1]*(1-exp(-beta[2]*x))
    solution = [2.3894212918E+02, 5.5015643181E-04]

    name, data, parameters, f, solution
end

function Chwirut2()
    name = "Chwirut2"
    data = [
    92.9000E0     0.500E0;
    57.1000E0     1.000E0;
    31.0500E0     1.750E0;
    11.5875E0     3.750E0;
    8.0250E0     5.750E0;
    63.6000E0     0.875E0;
    21.4000E0     2.250E0;
    14.2500E0     3.250E0;
    8.4750E0     5.250E0;
    63.8000E0     0.750E0;
    26.8000E0     1.750E0;
    16.4625E0     2.750E0;
    7.1250E0     4.750E0;
    67.3000E0     0.625E0;
    41.0000E0     1.250E0;
    21.1500E0     2.250E0;
    8.1750E0     4.250E0;
    81.5000E0      .500E0;
    13.1200E0     3.000E0;
    59.9000E0      .750E0;
    14.6200E0     3.000E0;
    32.9000E0     1.500E0;
    5.4400E0     6.000E0;
    12.5600E0     3.000E0;
    5.4400E0     6.000E0;
    32.0000E0     1.500E0;
    13.9500E0     3.000E0;
    75.8000E0      .500E0;
    20.0000E0     2.000E0;
    10.4200E0     4.000E0;
    59.5000E0      .750E0;
    21.6700E0     2.000E0;
    8.5500E0     5.000E0;
    62.0000E0      .750E0;
    20.2000E0     2.250E0;
    7.7600E0     3.750E0;
    3.7500E0     5.750E0;
    11.8100E0     3.000E0;
    54.7000E0      .750E0;
    23.7000E0     2.500E0;
    11.5500E0     4.000E0;
    61.3000E0      .750E0;
    17.7000E0     2.500E0;
    8.7400E0     4.000E0;
    59.2000E0      .750E0;
    16.3000E0     2.500E0;
    8.6200E0     4.000E0;
    81.0000E0      .500E0;
    4.8700E0     6.000E0;
    14.6200E0     3.000E0;
    81.7000E0      .500E0;
    17.1700E0     2.750E0;
    81.3000E0      .500E0;
    28.9000E0     1.750E0;
    ]

    parameters = [
    0.1 0.15;
    0.01        0.008;
    0.02        0.010;
    ]
    f(x, beta) = exp(-beta[1]*x)/(beta[2]+beta[3]*x)
    solution =[ 1.6657666537E-01, 5.1653291286E-03, 1.2150007096E-02]

    name, data, parameters, f, solution
end

function Chwitrut1()
    name = "Chwirut1"
    data =[
    92.9000E0     0.5000E0;
    78.7000E0     0.6250E0;
    64.2000E0     0.7500E0;
    64.9000E0     0.8750E0;
    57.1000E0     1.0000E0;
    43.3000E0     1.2500E0;
    31.1000E0     1.7500E0;
    23.6000E0     2.2500E0;
    31.0500E0     1.7500E0;
    23.7750E0     2.2500E0;
    17.7375E0     2.7500E0;
    13.8000E0     3.2500E0;
    11.5875E0     3.7500E0;
    9.4125E0     4.2500E0;
    7.7250E0     4.7500E0;
    7.3500E0     5.2500E0;
    8.0250E0     5.7500E0;
    90.6000E0     0.5000E0;
    76.9000E0     0.6250E0;
    71.6000E0     0.7500E0;
    63.6000E0     0.8750E0;
    54.0000E0     1.0000E0;
    39.2000E0     1.2500E0;
    29.3000E0     1.7500E0;
    21.4000E0     2.2500E0;
    29.1750E0     1.7500E0;
    22.1250E0     2.2500E0;
    17.5125E0     2.7500E0;
    14.2500E0     3.2500E0;
    9.4500E0     3.7500E0;
    9.1500E0     4.2500E0;
    7.9125E0     4.7500E0;
    8.4750E0     5.2500E0;
    6.1125E0     5.7500E0;
    80.0000E0     0.5000E0;
    79.0000E0     0.6250E0;
    63.8000E0     0.7500E0;
    57.2000E0     0.8750E0;
    53.2000E0     1.0000E0;
    42.5000E0     1.2500E0;
    26.8000E0     1.7500E0;
    20.4000E0     2.2500E0;
    26.8500E0     1.7500E0;
    21.0000E0     2.2500E0;
    16.4625E0     2.7500E0;
    12.5250E0     3.2500E0;
    10.5375E0     3.7500E0;
    8.5875E0     4.2500E0;
    7.1250E0     4.7500E0;
    6.1125E0     5.2500E0;
    5.9625E0     5.7500E0;
    74.1000E0     0.5000E0;
    67.3000E0     0.6250E0;
    60.8000E0     0.7500E0;
    55.5000E0     0.8750E0;
    50.3000E0     1.0000E0;
    41.0000E0     1.2500E0;
    29.4000E0     1.7500E0;
    20.4000E0     2.2500E0;
    29.3625E0     1.7500E0;
    21.1500E0     2.2500E0;
    16.7625E0     2.7500E0;
    13.2000E0     3.2500E0;
    10.8750E0     3.7500E0;
    8.1750E0     4.2500E0;
    7.3500E0     4.7500E0;
    5.9625E0     5.2500E0;
    5.6250E0     5.7500E0;
    81.5000E0      .5000E0;
    62.4000E0      .7500E0;
    32.5000E0     1.5000E0;
    12.4100E0     3.0000E0;
    13.1200E0     3.0000E0;
    15.5600E0     3.0000E0;
    5.6300E0     6.0000E0;
    78.0000E0      .5000E0;
    59.9000E0      .7500E0;
    33.2000E0     1.5000E0;
    13.8400E0     3.0000E0;
    12.7500E0     3.0000E0;
    14.6200E0     3.0000E0;
    3.9400E0     6.0000E0;
    76.8000E0      .5000E0;
    61.0000E0      .7500E0;
    32.9000E0     1.5000E0;
    13.8700E0     3.0000E0;
    11.8100E0     3.0000E0;
    13.3100E0     3.0000E0;
    5.4400E0     6.0000E0;
    78.0000E0      .5000E0;
    63.5000E0      .7500E0;
    33.8000E0     1.5000E0;
    12.5600E0     3.0000E0;
    5.6300E0     6.0000E0;
    12.7500E0     3.0000E0;
    13.1200E0     3.0000E0;
    5.4400E0     6.0000E0;
    76.8000E0      .5000E0;
    60.0000E0      .7500E0;
    47.8000E0     1.0000E0;
    32.0000E0     1.5000E0;
    22.2000E0     2.0000E0;
    22.5700E0     2.0000E0;
    18.8200E0     2.5000E0;
    13.9500E0     3.0000E0;
    11.2500E0     4.0000E0;
    9.0000E0     5.0000E0;
    6.6700E0     6.0000E0;
    75.8000E0      .5000E0;
    62.0000E0      .7500E0;
    48.8000E0     1.0000E0;
    35.2000E0     1.5000E0;
    20.0000E0     2.0000E0;
    20.3200E0     2.0000E0;
    19.3100E0     2.5000E0;
    12.7500E0     3.0000E0;
    10.4200E0     4.0000E0;
    7.3100E0     5.0000E0;
    7.4200E0     6.0000E0;
    70.5000E0      .5000E0;
    59.5000E0      .7500E0;
    48.5000E0     1.0000E0;
    35.8000E0     1.5000E0;
    21.0000E0     2.0000E0;
    21.6700E0     2.0000E0;
    21.0000E0     2.5000E0;
    15.6400E0     3.0000E0;
    8.1700E0     4.0000E0;
    8.5500E0     5.0000E0;
    10.1200E0     6.0000E0;
    78.0000E0      .5000E0;
    66.0000E0      .6250E0;
    62.0000E0      .7500E0;
    58.0000E0      .8750E0;
    47.7000E0     1.0000E0;
    37.8000E0     1.2500E0;
    20.2000E0     2.2500E0;
    21.0700E0     2.2500E0;
    13.8700E0     2.7500E0;
    9.6700E0     3.2500E0;
    7.7600E0     3.7500E0;
    5.4400E0     4.2500E0;
    4.8700E0     4.7500E0;
    4.0100E0     5.2500E0;
    3.7500E0     5.7500E0;
    24.1900E0     3.0000E0;
    25.7600E0     3.0000E0;
    18.0700E0     3.0000E0;
    11.8100E0     3.0000E0;
    12.0700E0     3.0000E0;
    16.1200E0     3.0000E0;
    70.8000E0      .5000E0;
    54.7000E0      .7500E0;
    48.0000E0     1.0000E0;
    39.8000E0     1.5000E0;
    29.8000E0     2.0000E0;
    23.7000E0     2.5000E0;
    29.6200E0     2.0000E0;
    23.8100E0     2.5000E0;
    17.7000E0     3.0000E0;
    11.5500E0     4.0000E0;
    12.0700E0     5.0000E0;
    8.7400E0     6.0000E0;
    80.7000E0      .5000E0;
    61.3000E0      .7500E0;
    47.5000E0     1.0000E0;
    29.0000E0     1.5000E0;
    24.0000E0     2.0000E0;
    17.7000E0     2.5000E0;
    24.5600E0     2.0000E0;
    18.6700E0     2.5000E0;
    16.2400E0     3.0000E0;
    8.7400E0     4.0000E0;
    7.8700E0     5.0000E0;
    8.5100E0     6.0000E0;
    66.7000E0      .5000E0;
    59.2000E0      .7500E0;
    40.8000E0     1.0000E0;
    30.7000E0     1.5000E0;
    25.7000E0     2.0000E0;
    16.3000E0     2.5000E0;
    25.9900E0     2.0000E0;
    16.9500E0     2.5000E0;
    13.3500E0     3.0000E0;
    8.6200E0     4.0000E0;
    7.2000E0     5.0000E0;
    6.6400E0     6.0000E0;
    13.6900E0     3.0000E0;
    81.0000E0      .5000E0;
    64.5000E0      .7500E0;
    35.5000E0     1.5000E0;
    13.3100E0     3.0000E0;
    4.8700E0     6.0000E0;
    12.9400E0     3.0000E0;
    5.0600E0     6.0000E0;
    15.1900E0     3.0000E0;
    14.6200E0     3.0000E0;
    15.6400E0     3.0000E0;
    25.5000E0     1.7500E0;
    25.9500E0     1.7500E0;
    81.7000E0      .5000E0;
    61.6000E0      .7500E0;
    29.8000E0     1.7500E0;
    29.8100E0     1.7500E0;
    17.1700E0     2.7500E0;
    10.3900E0     3.7500E0;
    28.4000E0     1.7500E0;
    28.6900E0     1.7500E0;
    81.3000E0      .5000E0;
    60.9000E0      .7500E0;
    16.6500E0     2.7500E0;
    10.0500E0     3.7500E0;
    28.9000E0     1.7500E0;
    28.9500E0     1.7500E0]

    parameters = [
    0.1         0.15;
    0.01        0.008;
    0.02        0.010;
    ]
    f(x, beta) = exp(-beta[1]*x)/(beta[2]+beta[3]*x)
    solution = [1.9027818370E-01, 6.1314004477E-03, 1.0530908399E-02]

    name, data, parameters, f, solution
end

function Lanczos3()
    name = "Lanczos3"
    data = [
    2.5134E+00  0.00000E+00;
    2.0443E+00  5.00000E-02;
    1.6684E+00  1.00000E-01;
    1.3664E+00  1.50000E-01;
    1.1232E+00  2.00000E-01;
    0.9269E+00  2.50000E-01;
    0.7679E+00  3.00000E-01;
    0.6389E+00  3.50000E-01;
    0.5338E+00  4.00000E-01;
    0.4479E+00  4.50000E-01;
    0.3776E+00  5.00000E-01;
    0.3197E+00  5.50000E-01;
    0.2720E+00  6.00000E-01;
    0.2325E+00  6.50000E-01;
    0.1997E+00  7.00000E-01;
    0.1723E+00  7.50000E-01;
    0.1493E+00  8.00000E-01;
    0.1301E+00  8.50000E-01;
    0.1138E+00  9.00000E-01;
    0.1000E+00  9.50000E-01;
    0.0883E+00  1.00000E+00;
    0.0783E+00  1.05000E+00;
    0.0698E+00  1.10000E+00;
    0.0624E+00  1.15000E+00]

    parameters = [
    1.2         0.5;
    0.3         0.7;
    5.6         3.6;
    5.5         4.2;
    6.5         4  ;
    7.6         6.3;
    ]
    f(x, beta) = beta[1] * exp(-beta[2] * x) + beta[3] * exp(-beta[4]*x) + beta[5] * exp(-beta[6]*x)
    solution = [
 8.6816414977E-02,
 9.5498101505E-01,
 8.4400777463E-01,
 2.9515951832E+00,
 1.5825685901E+00,
 4.9863565084E+00
 ]

    name, data, parameters, f, solution
end

function Gauss1()
    name = "Gauss1"
    data = [
    97.62227    1.000000;
    97.80724    2.000000;
    96.62247    3.000000;
    92.59022    4.000000;
    91.23869    5.000000;
    95.32704    6.000000;
    90.35040    7.000000;
    89.46235    8.000000;
    91.72520    9.000000;
    89.86916   10.000000;
    86.88076    11.00000;
    85.94360    12.00000;
    87.60686    13.00000;
    86.25839    14.00000;
    80.74976    15.00000;
    83.03551    16.00000;
    88.25837    17.00000;
    82.01316    18.00000;
    82.74098    19.00000;
    83.30034    20.00000;
    81.27850    21.00000;
    81.85506    22.00000;
    80.75195    23.00000;
    80.09573    24.00000;
    81.07633    25.00000;
    78.81542    26.00000;
    78.38596    27.00000;
    79.93386    28.00000;
    79.48474    29.00000;
    79.95942    30.00000;
    76.10691    31.00000;
    78.39830    32.00000;
    81.43060    33.00000;
    82.48867    34.00000;
    81.65462    35.00000;
    80.84323    36.00000;
    88.68663    37.00000;
    84.74438    38.00000;
    86.83934    39.00000;
    85.97739    40.00000;
    91.28509    41.00000;
    97.22411    42.00000;
    93.51733    43.00000;
    94.10159    44.00000;
    101.91760    45.00000;
    98.43134    46.00000;
    110.4214    47.00000;
    107.6628    48.00000;
    111.7288    49.00000;
    116.5115    50.00000;
    120.7609    51.00000;
    123.9553    52.00000;
    124.2437    53.00000;
    130.7996    54.00000;
    133.2960    55.00000;
    130.7788    56.00000;
    132.0565    57.00000;
    138.6584    58.00000;
    142.9252    59.00000;
    142.7215    60.00000;
    144.1249    61.00000;
    147.4377    62.00000;
    148.2647    63.00000;
    152.0519    64.00000;
    147.3863    65.00000;
    149.2074    66.00000;
    148.9537    67.00000;
    144.5876    68.00000;
    148.1226    69.00000;
    148.0144    70.00000;
    143.8893    71.00000;
    140.9088    72.00000;
    143.4434    73.00000;
    139.3938    74.00000;
    135.9878    75.00000;
    136.3927    76.00000;
    126.7262    77.00000;
    124.4487    78.00000;
    122.8647    79.00000;
    113.8557    80.00000;
    113.7037    81.00000;
    106.8407    82.00000;
    107.0034    83.00000;
    102.46290    84.00000;
    96.09296    85.00000;
    94.57555    86.00000;
    86.98824    87.00000;
    84.90154    88.00000;
    81.18023    89.00000;
    76.40117    90.00000;
    67.09200    91.00000;
    72.67155    92.00000;
    68.10848    93.00000;
    67.99088    94.00000;
    63.34094    95.00000;
    60.55253    96.00000;
    56.18687    97.00000;
    53.64482    98.00000;
    53.70307    99.00000;
    48.07893   100.00000;
    42.21258   101.00000;
    45.65181   102.00000;
    41.69728   103.00000;
    41.24946   104.00000;
    39.21349   105.00000;
    37.71696    106.0000;
    36.68395    107.0000;
    37.30393    108.0000;
    37.43277    109.0000;
    37.45012    110.0000;
    32.64648    111.0000;
    31.84347    112.0000;
    31.39951    113.0000;
    26.68912    114.0000;
    32.25323    115.0000;
    27.61008    116.0000;
    33.58649    117.0000;
    28.10714    118.0000;
    30.26428    119.0000;
    28.01648    120.0000;
    29.11021    121.0000;
    23.02099    122.0000;
    25.65091    123.0000;
    28.50295    124.0000;
    25.23701    125.0000;
    26.13828    126.0000;
    33.53260    127.0000;
    29.25195    128.0000;
    27.09847    129.0000;
    26.52999    130.0000;
    25.52401    131.0000;
    26.69218    132.0000;
    24.55269    133.0000;
    27.71763    134.0000;
    25.20297    135.0000;
    25.61483    136.0000;
    25.06893    137.0000;
    27.63930    138.0000;
    24.94851    139.0000;
    25.86806    140.0000;
    22.48183    141.0000;
    26.90045    142.0000;
    25.39919    143.0000;
    17.90614    144.0000;
    23.76039    145.0000;
    25.89689    146.0000;
    27.64231    147.0000;
    22.86101    148.0000;
    26.47003    149.0000;
    23.72888    150.0000;
    27.54334    151.0000;
    30.52683    152.0000;
    28.07261    153.0000;
    34.92815    154.0000;
    28.29194    155.0000;
    34.19161    156.0000;
    35.41207    157.0000;
    37.09336    158.0000;
    40.98330    159.0000;
    39.53923    160.0000;
    47.80123    161.0000;
    47.46305    162.0000;
    51.04166    163.0000;
    54.58065    164.0000;
    57.53001    165.0000;
    61.42089    166.0000;
    62.79032    167.0000;
    68.51455    168.0000;
    70.23053    169.0000;
    74.42776    170.0000;
    76.59911    171.0000;
    81.62053    172.0000;
    83.42208    173.0000;
    79.17451    174.0000;
    88.56985    175.0000;
    85.66525    176.0000;
    86.55502    177.0000;
    90.65907    178.0000;
    84.27290    179.0000;
    85.72220    180.0000;
    83.10702    181.0000;
    82.16884    182.0000;
    80.42568    183.0000;
    78.15692    184.0000;
    79.79691    185.0000;
    77.84378    186.0000;
    74.50327    187.0000;
    71.57289    188.0000;
    65.88031    189.0000;
    65.01385    190.0000;
    60.19582    191.0000;
    59.66726    192.0000;
    52.95478    193.0000;
    53.87792    194.0000;
    44.91274    195.0000;
    41.09909    196.0000;
    41.68018    197.0000;
    34.53379    198.0000;
    34.86419    199.0000;
    33.14787    200.0000;
    29.58864    201.0000;
    27.29462    202.0000;
    21.91439    203.0000;
    19.08159    204.0000;
    24.90290    205.0000;
    19.82341    206.0000;
    16.75551    207.0000;
    18.24558    208.0000;
    17.23549    209.0000;
    16.34934    210.0000;
    13.71285    211.0000;
    14.75676    212.0000;
    13.97169    213.0000;
    12.42867    214.0000;
    14.35519    215.0000;
    7.703309    216.0000;
    10.234410    217.0000;
    11.78315    218.0000;
    13.87768    219.0000;
    4.535700    220.0000;
    10.059280    221.0000;
    8.424824    222.0000;
    10.533120    223.0000;
    9.602255    224.0000;
    7.877514    225.0000;
    6.258121    226.0000;
    8.899865    227.0000;
    7.877754    228.0000;
    12.51191    229.0000;
    10.66205    230.0000;
    6.035400    231.0000;
    6.790655    232.0000;
    8.783535    233.0000;
    4.600288    234.0000;
    8.400915    235.0000;
    7.216561    236.0000;
    10.017410    237.0000;
    7.331278    238.0000;
    6.527863    239.0000;
    2.842001    240.0000;
    10.325070    241.0000;
    4.790995    242.0000;
    8.377101    243.0000;
    6.264445    244.0000;
    2.706213    245.0000;
    8.362329    246.0000;
    8.983658    247.0000;
    3.362571    248.0000;
    1.182746    249.0000;
    4.875359    250.0000;
    ]

    parameters = [
    97.0        94.0;
    0.009       0.0105;
    100.0        99.0;
    65.0        63.0;
    20.0        25.0;
    70.0        71.0 ;
    178.0       180.0;
    16.5        20.0;
    ]

    solution = [
    9.8778210871E+01,
    1.0497276517E-02,
    1.0048990633E+02,
    6.7481111276E+01,
    2.3129773360E+01,
    7.1994503004E+01,
    1.7899805021E+02,
    1.8389389025E+01
    ]
    f(x, beta) = beta[1] * exp(-beta[2]x) + beta[3]*exp(-(x - beta[4])^2 / beta[5]^2) + beta[6]*exp(-(x - beta[7])^2 / beta[8]^2)

    return name, data, parameters, f, solution
end

function Gauss2()
    name = "Gauss2"
    data = [
    97.58776    1.000000;
    97.76344    2.000000;
    96.56705    3.000000;
    92.52037    4.000000;
    91.15097    5.000000;
    95.21728    6.000000;
    90.21355    7.000000;
    89.29235    8.000000;
    91.51479    9.000000;
    89.60966   10.000000;
    86.56187    11.00000;
    85.55316    12.00000;
    87.13054    13.00000;
    85.67940    14.00000;
    80.04851    15.00000;
    82.18925    16.00000;
    87.24081    17.00000;
    80.79407    18.00000;
    81.28570    19.00000;
    81.56940    20.00000;
    79.22715    21.00000;
    79.43275    22.00000;
    77.90195    23.00000;
    76.75468    24.00000;
    77.17377    25.00000;
    74.27348    26.00000;
    73.11900    27.00000;
    73.84826    28.00000;
    72.47870    29.00000;
    71.92292    30.00000;
    66.92176    31.00000;
    67.93835    32.00000;
    69.56207    33.00000;
    69.07066    34.00000;
    66.53983    35.00000;
    63.87883    36.00000;
    69.71537    37.00000;
    63.60588    38.00000;
    63.37154    39.00000;
    60.01835    40.00000;
    62.67481    41.00000;
    65.80666    42.00000;
    59.14304    43.00000;
    56.62951    44.00000;
    61.21785    45.00000;
    54.38790    46.00000;
    62.93443    47.00000;
    56.65144    48.00000;
    57.13362    49.00000;
    58.29689    50.00000;
    58.91744    51.00000;
    58.50172    52.00000;
    55.22885    53.00000;
    58.30375    54.00000;
    57.43237    55.00000;
    51.69407    56.00000;
    49.93132    57.00000;
    53.70760    58.00000;
    55.39712    59.00000;
    52.89709    60.00000;
    52.31649    61.00000;
    53.98720    62.00000;
    53.54158    63.00000;
    56.45046    64.00000;
    51.32276    65.00000;
    53.11676    66.00000;
    53.28631    67.00000;
    49.80555    68.00000;
    54.69564    69.00000;
    56.41627    70.00000;
    54.59362    71.00000;
    54.38520    72.00000;
    60.15354    73.00000;
    59.78773    74.00000;
    60.49995    75.00000;
    65.43885    76.00000;
    60.70001    77.00000;
    63.71865    78.00000;
    67.77139    79.00000;
    64.70934    80.00000;
    70.78193    81.00000;
    70.38651    82.00000;
    77.22359    83.00000;
    79.52665    84.00000;
    80.13077    85.00000;
    85.67823    86.00000;
    85.20647    87.00000;
    90.24548    88.00000;
    93.61953    89.00000;
    95.86509    90.00000;
    93.46992    91.00000;
    105.8137    92.00000;
    107.8269    93.00000;
    114.0607    94.00000;
    115.5019    95.00000;
    118.5110    96.00000;
    119.6177    97.00000;
    122.1940    98.00000;
    126.9903    99.00000;
    125.7005   100.00000;
    123.7447   101.00000;
    130.6543   102.00000;
    129.7168   103.00000;
    131.8240   104.00000;
    131.8759   105.00000;
    131.9994    106.0000;
    132.1221    107.0000;
    133.4414    108.0000;
    133.8252    109.0000;
    133.6695    110.0000;
    128.2851    111.0000;
    126.5182    112.0000;
    124.7550    113.0000;
    118.4016    114.0000;
    122.0334    115.0000;
    115.2059    116.0000;
    118.7856    117.0000;
    110.7387    118.0000;
    110.2003    119.0000;
    105.17290    120.0000;
    103.44720    121.0000;
    94.54280    122.0000;
    94.40526    123.0000;
    94.57964    124.0000;
    88.76605    125.0000;
    87.28747    126.0000;
    92.50443    127.0000;
    86.27997    128.0000;
    82.44307    129.0000;
    80.47367    130.0000;
    78.36608    131.0000;
    78.74307    132.0000;
    76.12786    133.0000;
    79.13108    134.0000;
    76.76062    135.0000;
    77.60769    136.0000;
    77.76633    137.0000;
    81.28220    138.0000;
    79.74307    139.0000;
    81.97964    140.0000;
    80.02952    141.0000;
    85.95232    142.0000;
    85.96838    143.0000;
    79.94789    144.0000;
    87.17023    145.0000;
    90.50992    146.0000;
    93.23373    147.0000;
    89.14803    148.0000;
    93.11492    149.0000;
    90.34337    150.0000;
    93.69421    151.0000;
    95.74256    152.0000;
    91.85105    153.0000;
    96.74503    154.0000;
    87.60996    155.0000;
    90.47012    156.0000;
    88.11690    157.0000;
    85.70673    158.0000;
    85.01361    159.0000;
    78.53040    160.0000;
    81.34148    161.0000;
    75.19295    162.0000;
    72.66115    163.0000;
    69.85504    164.0000;
    66.29476    165.0000;
    63.58502    166.0000;
    58.33847    167.0000;
    57.50766    168.0000;
    52.80498    169.0000;
    50.79319    170.0000;
    47.03490    171.0000;
    46.47090    172.0000;
    43.09016    173.0000;
    34.11531    174.0000;
    39.28235    175.0000;
    32.68386    176.0000;
    30.44056    177.0000;
    31.98932    178.0000;
    23.63330    179.0000;
    23.69643    180.0000;
    20.26812    181.0000;
    19.07074    182.0000;
    17.59544    183.0000;
    16.08785    184.0000;
    18.94267    185.0000;
    18.61354    186.0000;
    17.25800    187.0000;
    16.62285    188.0000;
    13.48367    189.0000;
    15.37647    190.0000;
    13.47208    191.0000;
    15.96188    192.0000;
    12.32547    193.0000;
    16.33880    194.0000;
    10.438330    195.0000;
    9.628715    196.0000;
    13.12268    197.0000;
    8.772417    198.0000;
    11.76143    199.0000;
    12.55020    200.0000;
    11.33108    201.0000;
    11.20493    202.0000;
    7.816916    203.0000;
    6.800675    204.0000;
    14.26581    205.0000;
    10.66285    206.0000;
    8.911574    207.0000;
    11.56733    208.0000;
    11.58207    209.0000;
    11.59071    210.0000;
    9.730134    211.0000;
    11.44237    212.0000;
    11.22912    213.0000;
    10.172130    214.0000;
    12.50905    215.0000;
    6.201493    216.0000;
    9.019605    217.0000;
    10.80607    218.0000;
    13.09625    219.0000;
    3.914271    220.0000;
    9.567886    221.0000;
    8.038448    222.0000;
    10.231040    223.0000;
    9.367410    224.0000;
    7.695971    225.0000;
    6.118575    226.0000;
    8.793207    227.0000;
    7.796692    228.0000;
    12.45065    229.0000;
    10.61601    230.0000;
    6.001003    231.0000;
    6.765098    232.0000;
    8.764653    233.0000;
    4.586418    234.0000;
    8.390783    235.0000;
    7.209202    236.0000;
    10.012090    237.0000;
    7.327461    238.0000;
    6.525136    239.0000;
    2.840065    240.0000;
    10.323710    241.0000;
    4.790035    242.0000;
    8.376431    243.0000;
    6.263980    244.0000;
    2.705892    245.0000;
    8.362109    246.0000;
    8.983507    247.0000;
    3.362469    248.0000;
    1.182678    249.0000;
    4.875312    250.0000;
    ]

    parameters = [
    96.0        98.0;
    0.009       0.0105;
    103.0       103.0;
    106.0       105.0;
    18.0        20.0;
    72.0        73.0;
    151.0       150.0;
    18.0        20.0;
    ]
    f(x, beta) = beta[1] * exp(-beta[2]*x) + beta[3] * exp(-(x - beta[4])^2 / beta[5]^2) + beta[6]*exp(-(x - beta[7])^2 / beta[8]^2)

     solution = [
     9.9018328406E+01 ,
     1.0994945399E-02 ,
     1.0188022528E+02 ,
     1.0703095519E+02 ,
     2.3578584029E+01 ,
     7.2045589471E+01 ,
     1.5327010194E+02 ,
     1.9525972636E+01
     ]
    return name, data, parameters, f, solution
end

function DanWood()
    name = "DanWood"
    data = [
    2.138E0        1.309E0;
    3.421E0        1.471E0;
    3.597E0        1.490E0;
    4.340E0        1.565E0;
    4.882E0        1.611E0;
    5.660E0        1.680E0;
    ]
    parameters =[
    1           0.7;
    5           4;
    ]
    f(x, beta) = beta[1] * x^beta[2]
    solution =  [7.6886226176E-01,  3.8604055871E+00]

    return name, data, parameters, f, solution
end

function Misra1b()
    name = "Misra1b"
    data =[
    10.07E0      77.6E0;
    14.73E0     114.9E0;
    17.94E0     141.1E0;
    23.93E0     190.8E0;
    29.61E0     239.9E0;
    35.18E0     289.0E0;
    40.02E0     332.8E0;
    44.82E0     378.4E0;
    50.76E0     434.8E0;
    55.05E0     477.3E0;
    61.01E0     536.8E0;
    66.40E0     593.1E0;
    75.47E0     689.1E0;
    81.78E0     760.0E0;
    ]

    parameters = [
    500         300   ;
    0.0001      0.0002;
    ]

    f(x, beta) =  beta[1]* (1 - (1 / (1 + beta[2] * x / 2)^2))
    solution = [3.3799746163E+02, 3.9039091287E-04]

    return name, data, parameters, f, solution
end

function MGH09()
    name = "MGH09"
    data = [
    1.957000E-01    4.000000E+00;
    1.947000E-01    2.000000E+00;
    1.735000E-01    1.000000E+00;
    1.600000E-01    5.000000E-01;
    8.440000E-02    2.500000E-01;
    6.270000E-02    1.670000E-01;
    4.560000E-02    1.250000E-01;
    3.420000E-02    1.000000E-01;
    3.230000E-02    8.330000E-02;
    2.350000E-02    7.140000E-02;
    2.460000E-02    6.250000E-02;
    ]

    parameters = [
    25          0.25 ;
    39          0.39 ;
    41.5        0.415;
    39          0.39 ;
    ]
    f(x, beta) =  beta[1]*(x^2 + x*beta[2]) / (x^2 + x * beta[3] + beta[4])
    solution = [
    1.9280693458E-01,
    1.9128232873E-01,
    1.2305650693E-01,
    1.3606233068E-01
    ]
    return name, data, parameters, f, solution
end

function Thurber()
    name = "Thurber"
    data = [
    80.574E0      -3.067E0;
    84.248E0      -2.981E0;
    87.264E0      -2.921E0;
    87.195E0      -2.912E0;
    89.076E0      -2.840E0;
    89.608E0      -2.797E0;
    89.868E0      -2.702E0;
    90.101E0      -2.699E0;
    92.405E0      -2.633E0;
    95.854E0      -2.481E0;
    100.696E0      -2.363E0;
    101.060E0      -2.322E0;
    401.672E0      -1.501E0;
    390.724E0      -1.460E0;
    567.534E0      -1.274E0;
    635.316E0      -1.212E0;
    733.054E0      -1.100E0;
    759.087E0      -1.046E0;
    894.206E0      -0.915E0;
    990.785E0      -0.714E0;
    1090.109E0      -0.566E0;
    1080.914E0      -0.545E0;
    1122.643E0      -0.400E0;
    1178.351E0      -0.309E0;
    1260.531E0      -0.109E0;
    1273.514E0      -0.103E0;
    1288.339E0       0.010E0;
    1327.543E0       0.119E0;
    1353.863E0       0.377E0;
    1414.509E0       0.790E0;
    1425.208E0       0.963E0;
    1421.384E0       1.006E0;
    1442.962E0       1.115E0;
    1464.350E0       1.572E0;
    1468.705E0       1.841E0;
    1447.894E0       2.047E0;
    1457.628E0       2.200E0;
    ]

    parameters = [
    1000        1300     ;
    1000        1500     ;
    400         500      ;
    40          75       ;
    0.7         1        ;
    0.3         0.4      ;
    0.03        0.05     ;
    ]
    f(x, beta) =  (beta[1] + beta[2]*x + beta[3]*x^2 + beta[4]*x^3) / (1 + beta[5]*x + beta[6]*x^2 + beta[7]*x^3)
    solution = [
    1.2881396800E+03 ,
     1.4910792535E+03,
     5.8323836877E+02,
     7.5416644291E+01,
     9.6629502864E-01,
     3.9797285797E-01,
     4.9727297349E-02
     ]
    return name, data, parameters, f, solution
end

function BoxBOD()
    name = "BoxBOD"
    data = [
    109             1;
    149             2;
    149             3;
    191             5;
    213             7;
    224            10;
    ]

    parameters = [
    1           100;
    1          0.75;
    ]
    f(x, beta) =  beta[1]* (1 - exp(-beta[2]*x))
    solution = [2.1380940889E+02,  5.4723748542E-01]

    return name, data, parameters, f, solution
end

function Rat42()
    name = "Rat42"
    data = [
    8.930E0        9.000E0;
    10.800E0       14.000E0;
    18.590E0       21.000E0;
    22.330E0       28.000E0;
    39.350E0       42.000E0;
    56.110E0       57.000E0;
    61.730E0       63.000E0;
    64.620E0       70.000E0;
    67.080E0       79.000E0;
    ]

    parameters = [
    100         75;
    1          2.5;
    0.1        0.07;
    ]

    f(x, beta) = beta[1] / (1 + exp(beta[2] - beta[3]*x))
    solution = [ 7.2462237576E+01 , 2.6180768402E+00, 6.7359200066E-02]

    return name, data, parameters, f, solution
end

function MGH10()
    name = "MGH10"
    data = [
    3.478000E+04    5.000000E+01;
    2.861000E+04    5.500000E+01;
    2.365000E+04    6.000000E+01;
    1.963000E+04    6.500000E+01;
    1.637000E+04    7.000000E+01;
    1.372000E+04    7.500000E+01;
    1.154000E+04    8.000000E+01;
    9.744000E+03    8.500000E+01;
    8.261000E+03    9.000000E+01;
    7.030000E+03    9.500000E+01;
    6.005000E+03    1.000000E+02;
    5.147000E+03    1.050000E+02;
    4.427000E+03    1.100000E+02;
    3.820000E+03    1.150000E+02;
    3.307000E+03    1.200000E+02;
    2.872000E+03    1.250000E+02;
    ]

    parameters = [
    2         0.02 ;
    400000      4000 ;
    25000       250    ;
    ]

    f(x, beta) =  beta[1] * exp(beta[2] / (x + beta[3]))
    solution = [5.6096364710E-03, 6.1813463463E+03, 3.4522363462E+02    ]
    return name, data, parameters, f, solution
end

function Eckerle4()
    name = "Eckerle4"
    data = [
    0.0001575E0    400.000000E0;
    0.0001699E0    405.000000E0;
    0.0002350E0    410.000000E0;
    0.0003102E0    415.000000E0;
    0.0004917E0    420.000000E0;
    0.0008710E0    425.000000E0;
    0.0017418E0    430.000000E0;
    0.0046400E0    435.000000E0;
    0.0065895E0    436.500000E0;
    0.0097302E0    438.000000E0;
    0.0149002E0    439.500000E0;
    0.0237310E0    441.000000E0;
    0.0401683E0    442.500000E0;
    0.0712559E0    444.000000E0;
    0.1264458E0    445.500000E0;
    0.2073413E0    447.000000E0;
    0.2902366E0    448.500000E0;
    0.3445623E0    450.000000E0;
    0.3698049E0    451.500000E0;
    0.3668534E0    453.000000E0;
    0.3106727E0    454.500000E0;
    0.2078154E0    456.000000E0;
    0.1164354E0    457.500000E0;
    0.0616764E0    459.000000E0;
    0.0337200E0    460.500000E0;
    0.0194023E0    462.000000E0;
    0.0117831E0    463.500000E0;
    0.0074357E0    465.000000E0;
    0.0022732E0    470.000000E0;
    0.0008800E0    475.000000E0;
    0.0004579E0    480.000000E0;
    0.0002345E0    485.000000E0;
    0.0001586E0    490.000000E0;
    0.0001143E0    495.000000E0;
    0.0000710E0    500.000000E0;
    ]

    parameters = [
    1           1.5 ;
    10           5   ;
    500         450   ;
    ]
    f(x, beta) = beta[1] / beta[2] * exp(-(x - beta[3])^2 / (2*beta[2]^2))
    solution = [1.5543827178E+00    ,   4.0888321754E+00    ,4.5154121844E+02]

    return name, data, parameters, f, solution
end

function Rat43()
    name = "Rat43"
    data =[
    16.08E0     1.0E0;
    33.83E0     2.0E0;
    65.80E0     3.0E0;
    97.20E0     4.0E0;
    191.55E0     5.0E0;
    326.20E0     6.0E0;
    386.87E0     7.0E0;
    520.53E0     8.0E0;
    590.03E0     9.0E0;
    651.92E0    10.0E0;
    724.93E0    11.0E0;
    699.56E0    12.0E0;
    689.96E0    13.0E0;
    637.56E0    14.0E0;
    717.41E0    15.0E0;
    ]

    parameters = [
    100         700;
    10           5;
    1           0.75;
    1           1.3 ;
    ]
    f(x, beta) = beta[1] / (1 + exp(beta[2] - beta[3] * x))^(1/beta[4])
    solution = [
    6.9964151270E+02,
    5.2771253025E+00,
    7.5962938329E-01,
    1.2792483859E+00
  ]

    return name, data, parameters, f, solution
end

function Bennett5()

    name = "Bennet5"
    data = [
    -34.834702E0      7.447168E0;
    -34.393200E0      8.102586E0;
    -34.152901E0      8.452547E0;
    -33.979099E0      8.711278E0;
    -33.845901E0      8.916774E0;
    -33.732899E0      9.087155E0;
    -33.640301E0      9.232590E0;
    -33.559200E0      9.359535E0;
    -33.486801E0      9.472166E0;
    -33.423100E0      9.573384E0;
    -33.365101E0      9.665293E0;
    -33.313000E0      9.749461E0;
    -33.260899E0      9.827092E0;
    -33.217400E0      9.899128E0;
    -33.176899E0      9.966321E0;
    -33.139198E0     10.029280E0;
    -33.101601E0     10.088510E0;
    -33.066799E0     10.144430E0;
    -33.035000E0     10.197380E0;
    -33.003101E0     10.247670E0;
    -32.971298E0     10.295560E0;
    -32.942299E0     10.341250E0;
    -32.916302E0     10.384950E0;
    -32.890202E0     10.426820E0;
    -32.864101E0     10.467000E0;
    -32.841000E0     10.505640E0;
    -32.817799E0     10.542830E0;
    -32.797501E0     10.578690E0;
    -32.774300E0     10.613310E0;
    -32.757000E0     10.646780E0;
    -32.733799E0     10.679150E0;
    -32.716400E0     10.710520E0;
    -32.699100E0     10.740920E0;
    -32.678799E0     10.770440E0;
    -32.661400E0     10.799100E0;
    -32.644001E0     10.826970E0;
    -32.626701E0     10.854080E0;
    -32.612202E0     10.880470E0;
    -32.597698E0     10.906190E0;
    -32.583199E0     10.931260E0;
    -32.568699E0     10.955720E0;
    -32.554298E0     10.979590E0;
    -32.539799E0     11.002910E0;
    -32.525299E0     11.025700E0;
    -32.510799E0     11.047980E0;
    -32.499199E0     11.069770E0;
    -32.487598E0     11.091100E0;
    -32.473202E0     11.111980E0;
    -32.461601E0     11.132440E0;
    -32.435501E0     11.152480E0;
    -32.435501E0     11.172130E0;
    -32.426800E0     11.191410E0;
    -32.412300E0     11.210310E0;
    -32.400799E0     11.228870E0;
    -32.392101E0     11.247090E0;
    -32.380501E0     11.264980E0;
    -32.366001E0     11.282560E0;
    -32.357300E0     11.299840E0;
    -32.348598E0     11.316820E0;
    -32.339901E0     11.333520E0;
    -32.328400E0     11.349940E0;
    -32.319698E0     11.366100E0;
    -32.311001E0     11.382000E0;
    -32.299400E0     11.397660E0;
    -32.290699E0     11.413070E0;
    -32.282001E0     11.428240E0;
    -32.273300E0     11.443200E0;
    -32.264599E0     11.457930E0;
    -32.256001E0     11.472440E0;
    -32.247299E0     11.486750E0;
    -32.238602E0     11.500860E0;
    -32.229900E0     11.514770E0;
    -32.224098E0     11.528490E0;
    -32.215401E0     11.542020E0;
    -32.203800E0     11.555380E0;
    -32.198002E0     11.568550E0;
    -32.189400E0     11.581560E0;
    -32.183601E0     11.594420E0;
    -32.174900E0     11.607121E0;
    -32.169102E0     11.619640E0;
    -32.163300E0     11.632000E0;
    -32.154598E0     11.644210E0;
    -32.145901E0     11.656280E0;
    -32.140099E0     11.668200E0;
    -32.131401E0     11.679980E0;
    -32.125599E0     11.691620E0;
    -32.119801E0     11.703130E0;
    -32.111198E0     11.714510E0;
    -32.105400E0     11.725760E0;
    -32.096699E0     11.736880E0;
    -32.090900E0     11.747890E0;
    -32.088001E0     11.758780E0;
    -32.079300E0     11.769550E0;
    -32.073502E0     11.780200E0;
    -32.067699E0     11.790730E0;
    -32.061901E0     11.801160E0;
    -32.056099E0     11.811480E0;
    -32.050301E0     11.821700E0;
    -32.044498E0     11.831810E0;
    -32.038799E0     11.841820E0;
    -32.033001E0     11.851730E0;
    -32.027199E0     11.861550E0;
    -32.024300E0     11.871270E0;
    -32.018501E0     11.880890E0;
    -32.012699E0     11.890420E0;
    -32.004002E0     11.899870E0;
    -32.001099E0     11.909220E0;
    -31.995300E0     11.918490E0;
    -31.989500E0     11.927680E0;
    -31.983700E0     11.936780E0;
    -31.977900E0     11.945790E0;
    -31.972099E0     11.954730E0;
    -31.969299E0     11.963590E0;
    -31.963501E0     11.972370E0;
    -31.957701E0     11.981070E0;
    -31.951900E0     11.989700E0;
    -31.946100E0     11.998260E0;
    -31.940300E0     12.006740E0;
    -31.937401E0     12.015150E0;
    -31.931601E0     12.023490E0;
    -31.925800E0     12.031760E0;
    -31.922899E0     12.039970E0;
    -31.917101E0     12.048100E0;
    -31.911301E0     12.056170E0;
    -31.908400E0     12.064180E0;
    -31.902599E0     12.072120E0;
    -31.896900E0     12.080010E0;
    -31.893999E0     12.087820E0;
    -31.888201E0     12.095580E0;
    -31.885300E0     12.103280E0;
    -31.882401E0     12.110920E0;
    -31.876600E0     12.118500E0;
    -31.873699E0     12.126030E0;
    -31.867901E0     12.133500E0;
    -31.862101E0     12.140910E0;
    -31.859200E0     12.148270E0;
    -31.856300E0     12.155570E0;
    -31.850500E0     12.162830E0;
    -31.844700E0     12.170030E0;
    -31.841801E0     12.177170E0;
    -31.838900E0     12.184270E0;
    -31.833099E0     12.191320E0;
    -31.830200E0     12.198320E0;
    -31.827299E0     12.205270E0;
    -31.821600E0     12.212170E0;
    -31.818701E0     12.219030E0;
    -31.812901E0     12.225840E0;
    -31.809999E0     12.232600E0;
    -31.807100E0     12.239320E0;
    -31.801300E0     12.245990E0;
    -31.798401E0     12.252620E0;
    -31.795500E0     12.259200E0;
    -31.789700E0     12.265750E0;
    -31.786800E0     12.272240E0;
    ]

    parameters = [
    -2000       -1500;
    50          45;
    0.8         0.85;
    ]
    f(x, beta) = beta[1] *(beta[2] + x)^(-1/beta[3])
    solution = [    -2.5235058043E+03,4.6736564644E+01,9.3218483193E-01]
    return name, data, parameters, f, solution
end

function ff!(fcur, x, f, data)
    for i in 1:size(data, 1)
        fcur[i] = data[i, 1] - f(data[i, 2], x)
    end
end

tests =  [misra1a, Chwirut2, Chwitrut1, Lanczos3, Gauss1, Gauss2, DanWood, Misra1b, MGH09, Thurber, BoxBOD, Rat42, MGH10, Eckerle4, Rat43, Bennett5]

for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
    n = 0
    N = 0
    for test in tests
        global name, data, parameters, f, solution = test()
      #  @show solution
        for j in 1:size(parameters, 2)
            global nls = LeastSquaresProblem(x = parameters[:, j], f! = (fcur, x) -> ff!(fcur, x, f, data), output_length = size(data, 1))
            global r = optimize!(nls, optimizer(LeastSquaresOptim.QR()), x_tol = 1e-50, f_tol = 1e-36, g_tol = 1e-50)
            n += norm(r.minimizer - solution) <= 1e-3
            N += 1
            @test !isnan(mean(r.minimizer) )
        end
    end
    println("strd $optimizer_abbr $n / $N")
end


text:
include( "testingUtils.jl" )
include( "unitTestPredictEnsemble.jl" )
include( "unitTestargminProb.jl" )
include( "unitTestargmaxProb.jl" )
include( "unitTestBootStrapCorEstimation.jl" )
include( "unitTestBootstrapEstimation.jl" )
include( "unitTestDirichletPosterior.jl" )
include( "unitTestLossF.jl" )
include( "unitTestDirectSolution.jl" )


text:
function test_channel_pooling_layer(backend::Backend, pooling::PoolingFunction, tensor_dim::Int, n_input, T, eps)
  println("-- Testing ChannelPooling($(typeof(pooling))) on $(typeof(backend)){$T}...")

  dims = [abs(rand(Int, tensor_dim)) % 7 + 7 for i = 1:n_input]
  op_dim = max(abs(rand(Int)) % tensor_dim, 1)
  pad = (2,2)
  kernel = 3
  stride = 2

  println("    > Setup (pool along dimension $op_dim for $tensor_dim-D tensors)")

  layer = ChannelPoolingLayer(kernel=kernel, stride=stride, pad=pad, pooling=pooling,
      tops=Array(Symbol,n_input), bottoms=Array(Symbol,n_input), channel_dim=op_dim)

  input = [rand(T, dim...) for dim in dims]
  inputs = Blob[make_blob(backend, x) for x in input]
  diffs = Blob[make_blob(backend, x) for x in input]

  state = setup(backend, layer, inputs, diffs)

  println("    > Forward")
  forward(backend, state, inputs)

  payloads = Array(Any, n_input)
  for i = 1:n_input
    expected_output, payloads[i] = channel_pooling_forward(state, i, input[i], op_dim)
    got_output = to_array(state.blobs[i])
    @test all(-eps .< expected_output-got_output .< eps)
  end

  println("    > Backward")
  top_diff = [rand(T, size(state.blobs[i])) for i = 1:n_input]
  for i = 1:n_input
    copy!(state.blobs_diff[i], top_diff[i])
  end
  backward(backend, state, inputs, diffs)

  for i = 1:n_input
    expected_output = channel_pooling_backward(state, i, input[i], top_diff[i], payloads[i], op_dim)
    got_output = to_array(diffs[i])
    @test all(-eps .< expected_output - got_output .< eps)
  end

  shutdown(backend, state)
end

function channel_pooling_forward(state, i, input::Array, op_dim)
  dim_pre, dim_pool, dim_post = split_dims(input, op_dim)
  dim_pooled = size(state.blobs[i], op_dim)

  output = zeros(eltype(input), size(state.blobs[i]))
  if isa(state.layer.pooling, Pooling.Max)
    mask = similar(output, Int)
  end

  canonical_input = reshape(input, dim_pre, dim_pool, dim_post)
  canonical_output = reshape(output, dim_pre, dim_pooled, dim_post)
  if isa(state.layer.pooling, Pooling.Max)
    canonical_mask = reshape(mask, dim_pre, dim_pooled, dim_post)
  end

  for n = 1:dim_post
    for pc = 1:dim_pooled
      cstart = (pc-1)*state.layer.stride - state.layer.pad[1] + 1
      cend = min(cstart + state.layer.kernel - 1, size(input, op_dim))
      cstart = max(1, cstart)

      region = canonical_input[:,cstart:cend,n]
      if isa(state.layer.pooling, Pooling.Max)
        maxval, maxidx = findmax(region, 2)
        canonical_output[:,pc,n] = maxval
        canonical_mask[:,pc,n] = maxidx
      elseif isa(state.layer.pooling, Pooling.Mean)
        canonical_output[:,pc,n] = sum(region, 2) / state.layer.kernel
      else
        error("Unknown pooling $(state.layer.pooling)")
      end
    end
  end

  if isa(state.layer.pooling, Pooling.Max)
    return (output, mask)
  else
    return (output, nothing)
  end
end

function channel_pooling_backward(state, i, input::Array, diff::Array, payload::Any, op_dim)
  dim_pre, dim_pool, dim_post = split_dims(input, op_dim)
  dim_pooled = size(state.blobs[i], op_dim)

  gradient = zeros(eltype(input), size(input))
  canonical_input = reshape(input, dim_pre, dim_pool, dim_post)
  canonical_gradient = reshape(gradient, dim_pre, dim_pool, dim_post)
  canonical_diff = reshape(diff, dim_pre, dim_pooled, dim_post)
  if isa(state.layer.pooling, Pooling.Max)
    canonical_mask = reshape(payload, dim_pre, dim_pooled, dim_post)
  end
  for n = 1:dim_post
    for pc = 1:dim_pooled
      cstart = (pc-1)*state.layer.stride - state.layer.pad[1] + 1
      cend = min(cstart + state.layer.kernel - 1, size(input, op_dim))
      cstart = max(1, cstart)

      if isa(state.layer.pooling, Pooling.Max)
        region = view(canonical_gradient,1:dim_pre,cstart:cend,n)
        maxidx = canonical_mask[:,pc,n]
        region[vec(maxidx)] += vec(canonical_diff[:,pc,n])
      elseif isa(state.layer.pooling, Pooling.Mean)
        for c = cstart:cend
          canonical_gradient[:,c,n] += canonical_diff[:,pc,n] / state.layer.kernel
        end
      else
        error("Unknown pooling $(state.layer.pooling)")
      end
    end
  end
  return gradient
end

function test_channel_pooling_layer(backend::Backend, pooling::PoolingFunction, n_input, T, eps)
  for i = 2:6
    test_channel_pooling_layer(backend, pooling, i, n_input, T, eps)
  end
end
function test_channel_pooling_layer(backend::Backend, n_input, T, eps)
  test_channel_pooling_layer(backend, Pooling.Max(), n_input, T, eps)
  test_channel_pooling_layer(backend, Pooling.Mean(), n_input, T, eps)
end

function test_channel_pooling_layer(backend::Backend)
  test_channel_pooling_layer(backend, 1, Float32, 1e-4)
  test_channel_pooling_layer(backend, 3, Float64, 1e-8)
end

if test_cpu
  test_channel_pooling_layer(backend_cpu)
end
if test_gpu
  test_channel_pooling_layer(backend_gpu)
end


text:
!x=√parse(Int,x)%1
f(a,s="$a",l=length(s))=0∈1:~-l.|>i->s[l-i]<'1'||!s[i+1:end]+!s[1:i]

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
"""
    getdim(f::AbstractFunction)

Returns the dimension of the function `f`, i.e. the number of outputs.
"""
function getdim(::AbstractFunction) end

"""
    getfunction(f::AbstractFunction)

Returns any wrapped function inside the function `f`.
"""
function getfunction end