text:
using IndirectCoStateInit
using Test

@testset "IndirectCoStateInit.jl" begin
    # Write your tests here.
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DhcpOptions <: SwaggerModel
    dnsServers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: dnsServers

    function DhcpOptions(;dnsServers=nothing)
        o = new()
        validate_property(DhcpOptions, Symbol("dnsServers"), dnsServers)
        setfield!(o, Symbol("dnsServers"), dnsServers)
        o
    end
end # type DhcpOptions

const _property_map_DhcpOptions = Dict{Symbol,Symbol}(Symbol("dnsServers")=>Symbol("dnsServers"))
const _property_types_DhcpOptions = Dict{Symbol,String}(Symbol("dnsServers")=>"Vector{String}")
Base.propertynames(::Type{ DhcpOptions }) = collect(keys(_property_map_DhcpOptions))
Swagger.property_type(::Type{ DhcpOptions }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DhcpOptions[name]))}
Swagger.field_name(::Type{ DhcpOptions }, property_name::Symbol) =  _property_map_DhcpOptions[property_name]

function check_required(o::DhcpOptions)
    true
end

function validate_property(::Type{ DhcpOptions }, name::Symbol, val)
end

ext: jl


text:
module Evap
using ModelingToolkit

# @register Base.ifelse(cond, a, b)

function petpt(msalb, srad, tmax, tmin, xhlai, eo)
    td = (0.6*tmax)+(0.4*tmin)
    if (xhlai <= 0.0)
        albedo = msalb
    else
        albedo = 0.23-((0.23-msalb)*exp(-((0.75*xhlai))))
    end
    slang = srad*23.923
    eeq = (slang*(0.000204-(0.000183*albedo)))*(td+29.0)
    if (tmax > 35.0)
        eo = eeq*(((tmax-35.0)*0.05)+1.1)
    else
        eo = eeq*1.1
    end
    if (tmax < 5.0)
        eo = (eeq*0.01)*exp((0.18*(tmax+20.0)))
    else
        eo = eo
    end
    eo = max(eo, 0.0001)
end

op_ifelse(op, a, b) = (op) * a + (1-op) * b
# op_ifelse(op::Operation, a, b) = Operation(ifelse, [op, a,b])

function petpt′(msalb, srad, tmax, tmin, xhlai, eo)
    td = (0.6*tmax)+(0.4*tmin)
    albedo = op_ifelse(xhlai <= 0, msalb, 0.23-((0.23-msalb)*exp(-((0.75*xhlai)))))
    slang = srad*23.923
    eeq = (slang*(0.000204-(0.000183*albedo)))*(td+29.0)
    eo = op_ifelse(tmax > 35, eeq*(((tmax-35.0)*0.05)+1.1), 1.1eeq)
    eo = op_ifelse(tmax < 5, (eeq*0.01)*exp((0.18*(tmax+20.0))), eo)
    eo = max(eo, 0.0001)
end

function petasce(canht, doy, msalb, meevp, srad, tdew, tmax, tmin, windht, windrun, xhlai, xlat, xelev, eo)
    tavg = (tmax+tmin)/2.0
    patm = 101.3*(((293.0-(0.0065*xelev))/293.0)^5.26)
    psycon = 0.000665*patm
    udelta = (2503.0*exp(((17.27*tavg)/(tavg+237.3))))/((tavg+237.3)^2.0)
    emax = 0.6108*exp(((17.27*tmax)/(tmax+237.3)))
    emin = 0.6108*exp(((17.27*tmin)/(tmin+237.3)))
    es = (emax+emin)/2.0
    ea = 0.6108*exp(((17.27*tdew)/(tdew+237.3)))
    rhmin = max(20.0, min(80.0, ((ea/emax)*100.0)))
    if (xhlai <= 0.0)
        albedo = msalb
    else
        albedo = 0.23
    end
    rns = (1.0-albedo)*srad
    pie = π #3.14159265359
    dr = 1.0+(0.033*cos((((2.0*pie)/365.0)*doy)))
    ldelta = 0.409*sin(((((2.0*pie)/365.0)*doy)-1.39))
    ws = acos(-(((1.0*tan(((xlat*pie)/180.0)))*tan(ldelta))))
    ra1 = (ws*sin(((xlat*pie)/180.0)))*sin(ldelta)
    ra2 = (cos(((xlat*pie)/180.0))*cos(ldelta))*sin(ws)
    ra = (((24.0/pie)*4.92)*dr)*(ra1+ra2)
    rso = (0.75+(2e-05*xelev))*ra
    ratio = srad/rso
    if (ratio < 0.3)
        ratio = 0.3
    else
        ratio = ratio
    end
    if (ratio > 1.0)
        ratio = 1.0
    else
        ratio = ratio
    end
    fcd = (1.35*ratio)-0.35
    tk4 = (((tmax+273.16)^4.0)+((tmin+273.16)^4.0))/2.0
    rnl = ((4.901e-09*fcd)*(0.34-(0.14*sqrt(ea))))*tk4
    rn = rns-rnl
    g = 0.0
    windsp = (((windrun*1000.0)/24.0)/60.0)/60.0
    wind2m = windsp*(4.87/log(((67.8*windht)-5.42)))
    cn = 0.0
    cd = 0.0
    if (meevp == "A")
        cd = 0.38
    else
        cd = cd
    end
    if (meevp == "A")
        cn = 1600.0
    else
        cn = cn
    end
    if (meevp == "G")
        cd = 0.34
    else
        cd = cd
    end
    if (meevp == "G")
        cn = 900.0
    else
        cn = cn
    end
    refet = ((0.408*udelta)*(rn-g))+(((psycon*(cn/(tavg+273.0)))*wind2m)*(es-ea))
    refet = refet/(udelta+(psycon*(1.0+(cd*wind2m))))
    refet = max(0.0001, refet)
    skc = 0.8
    kcbmin = 0.3
    kcbmax = 1.2
    if (xhlai <= 0.0)
        kcb = 0.0
    else
        kcb = max(0.0, (kcbmin+((kcbmax-kcbmin)*(1.0-exp(-(((1.0*skc)*xhlai)))))))
    end
    wnd = max(1.0, min(wind2m, 6.0))
    cht = max(0.001, canht)
    kcmax = 0.0
    if (meevp == "A")
        kcmax = max(1.0, (kcb+0.05))
    else
        kcmax = kcmax
    end
    if (meevp == "G")
        kcmax = max((1.2+(((0.04*(wnd-2.0))-(0.004*(rhmin-45.0)))*((cht/3.0)^0.3))), (kcb+0.05))
    else
        kcmax = kcmax
    end
    if (kcb <= kcbmin)
        fc = 0.0
    else
        fc = ((kcb-kcbmin)/(kcmax-kcbmin))^(1.0+(0.5*canht))
    end
    fw = 1.0
    few = min((1.0-fc), fw)
    ke = max(0.0, min((1.0*(kcmax-kcb)), (few*kcmax)))
    eo = (kcb+ke)*refet
    eo = max(eo, 0.0001)
end

function petasce′(canht, doy, msalb, meevp, srad, tdew, tmax, tmin, windht, windrun, xhlai, xlat, xelev, eo)
    tavg = (tmax+tmin)/2.0
    patm = 101.3*(((293.0-(0.0065*xelev))/293.0)^5.26)
    psycon = 0.000665*patm
    udelta = (2503.0*exp(((17.27*tavg)/(tavg+237.3))))/((tavg+237.3)^2.0)
    emax = 0.6108*exp(((17.27*tmax)/(tmax+237.3)))
    emin = 0.6108*exp(((17.27*tmin)/(tmin+237.3)))
    es = (emax+emin)/2.0
    ea = 0.6108*exp(((17.27*tdew)/(tdew+237.3)))
    rhmin = max(20.0, min(80.0, ((ea/emax)*100.0)))
    albedo = op_ifelse(xhlai <= 0.0, msalb, 0.23)
    rns = (1.0-albedo)*srad
    pie = π #3.14159265359
    dr = 1.0+(0.033*cos((((2.0*pie)/365.0)*doy)))
    ldelta = 0.409*sin(((((2.0*pie)/365.0)*doy)-1.39))
    ws = acos(-(((1.0*tan(((xlat*pie)/180.0)))*tan(ldelta))))
    ra1 = (ws*sin(((xlat*pie)/180.0)))*sin(ldelta)
    ra2 = (cos(((xlat*pie)/180.0))*cos(ldelta))*sin(ws)
    ra = (((24.0/pie)*4.92)*dr)*(ra1+ra2)
    rso = (0.75+(2e-05*xelev))*ra
    ratio = srad/rso
    ratio = op_ifelse(ratio < 0.3, 0.3, ratio)
    ratio = op_ifelse(ratio > 1.0, 1.0, ratio)
    fcd = (1.35*ratio)-0.35
    tk4 = (((tmax+273.16)^4.0)+((tmin+273.16)^4.0))/2.0
    rnl = ((4.901e-09*fcd)*(0.34-(0.14*sqrt(ea))))*tk4
    rn = rns-rnl
    g = 0.0
    windsp = (((windrun*1000.0)/24.0)/60.0)/60.0
    wind2m = windsp*(4.87/log(((67.8*windht)-5.42)))
    cn = 0.0
    cd = 0.0
    if (meevp == "A")
        cd = 0.38
    else
        cd = cd
    end
    if (meevp == "A")
        cn = 1600.0
    else
        cn = cn
    end
    if (meevp == "G")
        cd = 0.34
    else
        cd = cd
    end
    if (meevp == "G")
        cn = 900.0
    else
        cn = cn
    end
    refet = ((0.408*udelta)*(rn-g))+(((psycon*(cn/(tavg+273.0)))*wind2m)*(es-ea))
    refet = refet/(udelta+(psycon*(1.0+(cd*wind2m))))
    refet = max(0.0001, refet)
    skc = 0.8
    kcbmin = 0.3
    kcbmax = 1.2
    kcb = op_ifelse(xhlai <= 0.0,
                    0.0,
                    max(0.0,
                        (kcbmin+
                         ((kcbmax-
                           kcbmin)*
                          (1.0-
                           exp(-(((1.0*skc)*
                                  xhlai))))))))
    wnd = max(1.0, min(wind2m, 6.0))
    cht = max(0.001, canht)
    kcmax = 0.0
    if (meevp == "A")
        kcmax = max(1.0, (kcb+0.05))
    else
        kcmax = kcmax
    end
    if (meevp == "G")
        kcmax = max((1.2+(((0.04*(wnd-2.0))-(0.004*(rhmin-45.0)))*((cht/3.0)^0.3))), (kcb+0.05))
    else
        kcmax = kcmax
    end
    fc = op_ifelse(kcb <= kcbmin, 0.0,((kcb-kcbmin)/(kcbmax-kcbmin))^(1.0+(0.5*canht)))
    fw = 1.0
    few = min((1.0-fc), fw)
    ke = max(0.0, min((1.0*(kcmax-kcb)), (few*kcmax)))
    eo = (kcb+ke)*refet
    eo = max(eo, 0.0001)
end
end

using ModelingToolkit

@variables msalb srad tmax tmin xhlai eo
ptexpr = Evap.petpt′(msalb, srad, tmax, tmin, xhlai, eo)
@show vals = (msalb = 1, srad = 1 , tmax = 6, tmin=0, xhlai=1, eo=1)
@show Symbolic.apply(ptexpr, vals)

valsASCE = (canht = 1, doy=15,
            msalb=5, meevp=0.2, srad=1,
            tdew=24, tmax=6, tmin=0,
            windht=10, windrun=2,
            xhlai=1, xlat=0.36, xelev=200,
            eo=1)

meevp = "A"

@variables(canht, doy, msalb, srad, tdew, tmax, tmin,
           windht, windrun, xhlai, xlat, xelev, eo)

asceex = Evap.petasce′(canht, doy, msalb, meevp, srad, tdew,
                  tmax, tmin, windht, windrun, xhlai,
                       xlat, xelev, eo)

@show asceval = Symbolic.apply(asceex, valsASCE)

# using autodiff

# function petptfunc(x)
#     vals = (msalb = 1, srad = 1 , tmax = 6, tmin=0, xhlai=x[1], eo=1);
#     Evap.petpt(vals...)
# end

using ForwardDiff

ptgrad = map(0:0.05:1) do x
    xarr = [1,x,6,0,1,1]
    (input=xarr, value=Evap.petpt(xarr...), gradient=ForwardDiff.gradient(t -> Evap.petpt(t...), xarr))
end

function gradsweep(m, dim::Int, x, perturbations)
    mgrad = map(perturbations) do δ
        xarr = copy(x)
        xarr[dim] += δ
        (input=tuple(xarr...),
         value=m(xarr...),
         gradient=ForwardDiff.gradient(t -> m(t...), xarr)[dim])
    end
    return mgrad
end

x₀= [0.0,1,6,0,1,1]
@show gradsweep(Evap.petpt, 1, x₀, 0:0.05:1)

using Statistics
using LinearAlgebra

function linearity(m, dim::Int, x, perturbations)
    y = gradsweep(m, dim, x, perturbations)
    yvals = map(x->x.value, y)
    grads = map(x->x.gradient, y)
    return (dim=dim,ȳ = mean(yvals), ∇̄=mean(grads), σ_∇=std(grads)/mean(grads))
end

function linearity(m, x₀)
    ForwardDiff.hessian(x->m(x...), x₀) |> diag
end

lintests = map(1:6) do i
    linearity(Evap.petpt, i, x₀, 0:0.05:1)
end
@show lintests
@show linearity(Evap.petpt, x₀)

valsASCE = (canht = 1, doy=15,
            msalb=1, meevp=0.2, srad=1,
            tdew=24, tmax=6, tmin=0,
            windht=10, windrun=2,
            xhlai=1, xlat=0.36, xelev=200,
            eo=1)
@show Evap.petasce′(valsASCE...)
ascehess = linearity(Evap.petasce′, collect(valsASCE))
@show ascehess
# lintests2 = map(1:6) do i
#     linearity(Evap.petpt, i, [0.0,2,6,0,1,1], 0:0.05:1)
# end
# @show lintests2

# lintests3 = map(1:6) do i
#     linearity(Evap.petpt, i, [1.0,1,1,1,1,1], 0:0.05:1)
# end
# @show lintests3
curvatur = map(1:50) do i
    x = collect(valsASCE) .+ randn(Float64, 14)
    y = try
        linearity(Evap.petasce′, x) ./ x
    catch ex ;
        @info "failure"
    end
end |>
    y-> filter(x->x!=nothing, y) |>
    mean

@show curvatur

@show expand_derivatives(
    Differential(msalb)(
        Symbolic.apply(ptexpr,
                       (msalb=msalb,
                        xhlai=1,
                        tmin=0,
                        tmax=6,
                        srad=1,
                        eo=1))))
Dmsalb = Differential(msalb)
@show dpetdmsalb = expand_derivatives(
    Dmsalb(Symbolic.apply(ptexpr,
                          (msalb=msalb,
                           xhlai=1,
                           tmin=0,
                           tmax=tmax,
                           srad=1,
                           eo=1))))
@show Symbolic.apply(dpetdmsalb,
                     (msalb=1, tmax=4))


text:
module EvolutionaryModelingTools

using RecordedArrays
using RecordedArrays.ArrayInterface: indices
using Random
using MacroTools
using Requires

export @cfunc, @ufunc, @reaction, @quickloop, @reaction_eq
export Reaction, gillespie, gillespie!, indices

include("tools.jl")
include("sample.jl")
include("model.jl")
include("scalar.jl")

function __init__()
    @require LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" begin
        DEFAULT_LOOP_MODE[1] = true
    end
end

end # module


text:
module HurdleDMR

using Lasso, StatsBase, StatsModels, Missings, LambertW,
    SparseArrays, Distributed, SharedArrays, LinearAlgebra, Logging
using GLM: FPVector, FP, Link
using LoggingExtras: MinLevelLogger
using Tables: columntable
using StatsModels: TupleTerm, ColumnTable, missing_omit, schema, apply_schema, termvars
import StatsModels: modelcols
import Tables

export dmr, dmrpaths, hdmr, hdmrpaths, fit, coef, srproj, srprojX, @~, mcdmr, posindic,
  DCR, DMR, HDMR, DMRCoefs, DMRPaths, HDMRCoefs, HDMRPaths, @model, Model,
  CIR, predict, coeffwd, coefbwd, MinCVKfold,
  hasintercept, ncategories, nobs, ncoefs, ncovars, ncovarszero, ncovarspos, ncoefszero, ncoefspos,
  Hurdle, PositivePoisson, LogProductLogLink, logpdf_exact, logpdf_approx,
  InclusionRepetition, TableCountsRegressionModel

##############################################
# hurdle glm model involves
#  1. a choice model for the 0 or positive choice (e.g. binomial with logit link)
#  2. a positive (truncated) count model for positive counts (e.g. poisson with logit link)
#############################################

include("positive_poisson.jl")
include("hurdle.jl")
include("inrep.jl")
include("model.jl")
include("dmr.jl")
include("hdmr.jl")
include("srproj.jl")
include("invreg.jl")
include("multicounts.jl")

end


text:
"""
    DistributedOptimiser(optimiser::AbstractOptimiser)

Wrap the `optimiser` in a `DistributedOptimiser`. Before updating the
parameters, this averages the gradients across the processes using
non-blocking Allreduce
"""
struct DistributedOptimiser{B,O<:AbstractOptimiser}
    optimiser::O
    # TODO: Check if non-blocking calls are better. For small scale problems blocking calls are faster
    #       but that could be because we are unable to use Iallreduce with CuArrays
    function DistributedOptimiser(optimiser::O; blocking_communication::Bool=true) where {O<:AbstractOptimiser}
        return new{blocking_communication,O}(optimiser)
    end
end

function getproperty(opt::DistributedOptimiser, name::Symbol)
    if name == :optimiser || name == :optimizer
        return getfield(opt, :optimiser)
    else
        return getproperty(opt.optimiser, name)
    end
end

function setproperty!(opt::DistributedOptimiser, name::Symbol, x::Any)
    if name == :optimiser || name == :optimizer
        error("Cannot update `optimiser` property of `DistributedOptimiser`")
    else
        setproperty!(opt.optimiser, name, x)
    end
end

_get_request_type(ps::Params) = _get_request_type(first(ps))
_get_request_type(::CuArray) = JuliaTaskRequest
_get_request_type(::AbstractArray) = Request

function update!(opt::DistributedOptimiser{false}, xs::Params, gs)
    request_count = 0
    requests = Vector{_get_request_type(xs)}(undef, length(xs))
    # Add the gradients across all processes
    s = total_workers()
    for x in xs
        g = gs[x]
        g === nothing && continue
        # Average out the gradients
        g ./= s
        _, request = Iallreduce!(g, +, MPI.COMM_WORLD)
        request_count += 1
        requests[request_count] = request
    end
    # Wait for all the non-blocking operations to be completed
    Waitall!(requests[1:request_count])
    # Update the parameters
    update!(opt.optimiser, xs, gs)
    return
end

function update!(opt::DistributedOptimiser{true}, xs::Params, gs)
    # Add the gradients across all processes
    s = total_workers()
    for x in xs
        g = gs[x]
        g === nothing && continue
        # Average out the gradients
        g ./= s
        Allreduce!(g, +, MPI.COMM_WORLD)
    end
    # Update the parameters
    update!(opt.optimiser, xs, gs)
    return
end


text:
# Functions to access fields nested within higher-level types

"""
    networks(stationxml) -> ::Vector{Network}

Return the networks contained within a `FDSNStationXML` object.
"""
networks(sxml::FDSNStationXML) = sxml.network

"""
    stations(stationxml) -> ::Vector{Station}
    stations(network) -> ::Vector{Station}

Return a flattened list of all stations within the entire FDSNStationXML object, or within
a `network`.
"""
stations(sxml::FDSNStationXML) = collect(Iterators.flatten(stations(net) for net in sxml.network))
stations(net::Network) = net.station

"""
    channels(stationxml) -> ::Vector{Channel}
    channels(network) -> ::Vector{Channel}
    channels(station) -> ::Vector{Channel}

Return a flattened list of channels contained within a FDSNStationXML object, a whole
`network`, or a single `station`.
"""
channels(sxml::FDSNStationXML) = collect(Iterators.flatten(channels(net) for net in sxml.network))
channels(net::Network) = collect(Iterators.flatten(channels(sta) for sta in net.station))
channels(sta::Station) = sta.channel

Base.getproperty(net::AbstractArray{Network}, f::Symbol) = getfield.(net, f)
Base.getproperty(sta::AbstractArray{Station}, f::Symbol) = getfield.(sta, f)
Base.getproperty(cha::AbstractArray{Channel}, f::Symbol) = getfield.(cha, f)

"""
    channel_codes(network) -> ::Vector{String}
    channel_codes(stationxml) -> ::Vector{String}

Return a list of the channel codes of all the channels within
a `network` or `stationxml` document.
"""
function channel_codes(network::Network)
    cha_codes = String[]
    net = network.code
    for station in stations(network)
        sta = station.code
        for channel in channels(station)
            cha = channel.code
            loc = coalesce(channel.location_code, "")
            push!(cha_codes, "$net.$sta.$loc.$cha")
        end
    end
    cha_codes
end

channel_codes(sxml::FDSNStationXML) = reduce(vcat, channel_codes.(networks(sxml)))


text:
# Integrations with DFTK.jl

# -----------------------------------------------------------------------------
# Integration with AtomsBase
# -----------------------------------------------------------------------------

# ! all functions in this section are copied directly from experimental DFTK code
# ! this code will be deleted once it is live in DFTK

function parse_system(system::AbstractSystem{D}) where {D}
    if !all(periodicity(system))
        error("DFTK only supports calculations with periodic boundary conditions.")
    end

    # Parse abstract system and return data required to construct model
    mtx = austrip.(hcat(bounding_box(system)...))
    T = eltype(mtx)
    lattice = zeros(T, 3, 3)
    lattice[1:D, 1:D] .= mtx

    # Cache for instantiated pseudopotentials
    # (such that the respective objects are indistinguishable)
    cached_pseudos = Dict{String,Any}()
    atoms = map(system) do atom
        if hasproperty(atom, :potential)
            potential = atom.potential
        elseif hasproperty(atom, :pseudopotential)
            pspkey = atom.pseudopotential
            potential = get!(cached_pseudos, pspkey) do
                ElementPsp(AtomsBase.atomic_symbol(atom); psp = load_psp(pspkey))
            end
        else
            potential = ElementCoulomb(AtomsBase.atomic_symbol(atom))
        end

        coordinate = zeros(T, 3)
        coordinate[1:D] = inv(lattice[1:D, 1:D]) * T.(austrip.(position(atom)))
        potential => Vec3{T}(coordinate)
    end

    oldatoms = oldatoms_from_new(atoms)

    (; lattice, atoms = oldatoms)
end
function oldatoms_from_new(atomic_potentials)
    potentials = first.(atomic_potentials)
    potential_groups = [findall(Ref(pot) .== potentials) for pot in Set(potentials)]
    [first(atomic_potentials[first(group)]) => last.(atomic_potentials[group]) for group in potential_groups]
end
function DFTK.model_LDA(system::AbstractSystem; kwargs...)
    parsed = parse_system(system)
    model_LDA(parsed.lattice, parsed.atoms; kwargs...)
end

# -----------------------------------------------------------------------------
# Integration with InteratomicPotentials
# -----------------------------------------------------------------------------

# ! This integration should ultimately move to the DFTK package itself

@kwdef struct DFTKPotential <: ArbitraryPotential
    Ecut::Real
    kgrid::AbstractVector{<:Integer}
    n_bands::Union{Integer,Nothing} = nothing
    tol::Union{AbstractFloat,Nothing} = nothing
    damping::Union{AbstractFloat,Nothing} = nothing
    mixing::Union{Mixing,Nothing} = nothing
    previous_scfres::Ref{Any} = Ref{Any}()
    potential_energy_cache::Dict{Float64,Float64} = Dict{Float64,Float64}()
end
function DFTKPotential(Ecut::Real, kgrid::AbstractVector{<:Integer}; kwargs...)
    DFTKPotential(; Ecut = Ecut, kgrid = kgrid, kwargs...)
end
function DFTKPotential(Ecut::Unitful.Energy, kgrid::AbstractVector{<:Integer}; kwargs...)
    DFTKPotential(; Ecut = austrip(Ecut), kgrid = kgrid, kwargs...)
end

function calculate_scf(system::AbstractSystem, potential::DFTKPotential)
    model = model_LDA(system)
    basis = PlaneWaveBasis(model; Ecut = potential.Ecut, kgrid = potential.kgrid)

    args = (f => getfield(potential, f) for f ∈ (:n_bands, :tol, :damping, :mixing) if !isnothing(getfield(potential, f)))
    extra_args = isassigned(potential.previous_scfres) ? (ψ = potential.previous_scfres[].ψ, ρ = potential.previous_scfres[].ρ) : (;)
    scfres = self_consistent_field(basis; args..., extra_args...)
    potential.previous_scfres[] = scfres
end

function InteratomicPotentials.potential_energy(system::AbstractSystem, potential::DFTKPotential)
    calculate_scf(system, potential).energies.total
end

function InteratomicPotentials.potential_energy(system::DynamicSystem, potential::DFTKPotential)
    get!(potential.potential_energy_cache, austrip(system.time)) do
        InteratomicPotentials.potential_energy(system.system, potential)
    end
end

function InteratomicPotentials.force(system::AbstractSystem, potential::DFTKPotential)
    # TODO: support multiple species
    compute_forces_cart(calculate_scf(system, potential))[1]
end

function InteratomicPotentials.force(system::DynamicSystem, potential::DFTKPotential)
    scf = calculate_scf(system, potential)
    potential.potential_energy_cache[austrip(system.time)] = scf.energies.total
    # TODO: support multiple species
    compute_forces_cart(scf)[1]
end

# -----------------------------------------------------------------------------
# Experimental functions
# -----------------------------------------------------------------------------

# ! This function is experimental and will eventually be removed

function analyze_convergence(system::AbstractSystem, potential::DFTKPotential, cutoffs::AbstractVector{<:Unitful.Energy})
    energies = zeros(Float64, length(cutoffs))
    for i ∈ 1:length(cutoffs)
        parameters = DFTKPotential(
            Ecut = cutoffs[i],
            kgrid = potential.kgrid,
            n_bands = potential.n_bands,
            tol = potential.tol,
            damping = potential.damping,
            mixing = potential.mixing
        )
        @info "Ecut: $(cutoffs[i])"
        energies[i] = InteratomicPotentials.potential_energy(system, parameters)
    end

    plot(
        title = "DFTK Analysis",
        xlab = "Ecut",
        ylab = "Total Energy",
        legend = false,
        cutoffs,
        energies * ENERGY_UNIT
    )
end


text:
@testset "basic bandstructures" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, range = 1/√3))
    b = bandstructure(h, resolution = 13)
    @test length(bands(b)) == 1

    h = LatticePresets.honeycomb() |>
        hamiltonian(onsite(0.5, sublats = :A) + onsite(-0.5, sublats = :B) +
                    hopping(-1, range = 1/√3))
    b = bandstructure(h, resolution = 13)
    @test length(bands(b)) == 2

    h = LatticePresets.square() |> hamiltonian(hopping(1)) |> unitcell(2)
    b = bandstructure(h, resolution = 13)
    @test length(bands(b)) == 4
end

@testset "functional bandstructures" begin
    hc = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, sublats = :A=>:B, plusadjoint = true))
    matrix = similarmatrix(hc, LinearAlgebraPackage())
    hf((x,)) = bloch!(matrix, hc, (x, -x))
    mesh = marchingmesh(range(0, 1, length = 13))
    b = bandstructure(hf, mesh)
    @test length(bands(b)) == 2

    hc2 = LatticePresets.honeycomb() |> hamiltonian(hopping(-1))
    hp2 = parametric(hc2, @hopping!((t; s) -> s*t))
    matrix2 = similarmatrix(hc2, LinearAlgebraPackage())
    hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x))
    mesh2 = marchingmesh(range(0, 1, length = 13), range(0, 1, length = 13))
    b = bandstructure(hf2, mesh2)
    @test length(bands(b)) == 2
end

@testset "bandstructures cuts & transforms" begin
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2) + hopping(-1, range = 1/√3))
    mesh1D = marchingmesh(range(0, 2π, length = 13))
    b = bandstructure(h, mesh1D, cut = φ -> (φ, -φ), transform = inv)
    b´ = transform!(inv, bandstructure(h, mesh1D, cut = φ -> (φ, -φ)))
    @test length(bands(b)) == length(bands(b´)) == 2
    @test vertices(bands(b)[1]) == vertices(bands(b´)[1])
    h´ = unitcell(h)
    s1 = spectrum(h´, transform = inv)
    s2 = transform!(inv,spectrum(h´))
    @test energies(s1) == energies(s2)
end

@testset "parametric bandstructures" begin
    ph = LatticePresets.linear() |> hamiltonian(hopping(-I), orbitals = Val(2)) |> unitcell(2) |>
         parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k)-> t - k*I))
    mesh2D = marchingmesh(range(0, 1, length = 13), range(0, 2π, length = 13))
    mesh1D = marchingmesh(range(0, 2π, length = 13))
    b = bandstructure(ph, mesh2D)
    @test length(bands(b)) == 4
    b = bandstructure(ph, mesh2D, cut = (k, φ) -> (k + φ/2π, φ))
    @test length(bands(b)) == 4
    b = bandstructure(ph, mesh2D, cut = (k, φ) -> (k + φ/2π, φ))
    @test length(bands(b)) == 4
    b = bandstructure(ph, mesh1D, cut = k -> (k, 2π*k))
    @test length(bands(b)) == 4
    b = bandstructure(ph, mesh1D, cut = φ -> (φ, -φ))
    @test length(bands(b)) == 4
end


text:
struct AtomicVars
    # number of species
    nspecies::Int64
    # number of atoms for each species
    natoms::Vector{Int64} # of length (maxspecies)
    # maximum number of atoms over all the species
    natmmax::Int64
    # total number of atoms
    natmtot::Int64
    # index to atoms and species
    idxas::Array{Int64,2} # (maxatoms,maxspecies)
    # inverse atoms and species indices
    idxis::Vector{Int64} # (maxatoms*maxspecies)
    idxia::Vector{Int64} # (maxatoms*maxspecies)
    # molecule is .true. is the system is an isolated molecule
    molecule::Bool
    # primcell is .true. if primitive unit cell is to be found automatically
    primcell::Bool
    # atomic positions in lattice coordinates
    atposl::Array{Float64,3} #(3,maxatoms,maxspecies)
    # atomic positions in Cartesian coordinates
    atposc::Array{Float64,3} # (3,maxatoms,maxspecies)
    # magnitude of random displacements added to the atomic positions
    rndatposc::Float64
end

function AtomicVars(
    nspecies::Int64,
    natoms_::Vector{Int64},
    atposl_::Array{Float64,3},
    lattice_vars::LatticeVars
)

    epslat = lattice_vars.epslat
    avec = lattice_vars.avec

    natoms = copy(natoms_)

    # FIXME: Should not have this limitation
    maxspecies = 8 # maximum allowed species
    maxatoms = 200 # maximum allowed atoms per species

    atposl = copy(atposl_)
    atposc = zeros(3,maxatoms,maxspecies)

    for is in 1:nspecies
        for ia in 1:natoms[is]
            # map atomic lattice coordinates to [0,1)
            @views r3frac!( epslat, atposl[:,ia,is] )
            # determine atomic Cartesian coordinates
            @views r3mv!( avec,atposl[:,ia,is], atposc[:,ia,is])
        end
    end

    molecule = false
    rndatposc = 0.0
    primcell = false

    idxas = zeros(Int64, maxatoms, maxspecies)
    # inverse atoms and species indices
    idxis = zeros(Int64, maxatoms*maxspecies)
    idxia = zeros(Int64, maxatoms*maxspecies)
    natmmax = 0
    ias = 0
    for is in 1:nspecies
        for ia in 1:natoms[is]
            ias = ias + 1
            idxas[ia,is] = ias
            idxis[ias] = is
            idxia[ias] = ia
        end
        # maximum number of atoms over all species
        natmmax = max(natmmax, natoms[is])
    end
    # total number of atoms
    natmtot = ias

    return AtomicVars(
        nspecies, natoms, natmmax, natmtot,
        idxas, idxis, idxia,
        molecule, primcell, atposl, atposc, rndatposc
    )

end


text:
include("../src/pValues.jl")
using Test

test_dict = Dict(1 => Float64[1,2,3,4])

#These tests need updating
@test FindpValues(test_dict, 2.5, 1) == Dict(1 => 0.5)
@test FindpValues(test_dict, 6, 1) == Dict(1 => 0)
@test FindpValues(test_dict, 0, 1) == Dict(1 => 1)

test_dict = Dict(1 => Float64[1,2,3,4], 2 => Float64[5,6,7,8])
@test FindpValues(test_dict, 2.5, 1) == Dict(1 => 0.5, 2 => 1)

#This test fails but output is correct
@test FindpValues(test_dict, 2.5, 40) == Dict(1 => missing, 2 => missing)

meanArr = [1,2,3,4]
mean_val = 1
threshold = 1

@test GlobalpValues(meanArr, mean_val, 1.5) == 0.75

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
@testset "lz4hc" begin
    test_in = """
        Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor
        incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud
        exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat
        """
    test_size = convert(Cint, length(test_in))
    bound = CodecLz4.LZ4_compressBound(test_size)
    bufsize = bound
    half_bufsize = floor(Cint, bufsize/2)
    buffer = Vector{UInt8}(undef, bufsize)
    dec_buffer = Vector{UInt8}(undef, test_size)

    @testset "compress" begin
        result = CodecLz4.LZ4_compress_HC(pointer(test_in), pointer(buffer), test_size, bufsize, 13)

        result = CodecLz4.LZ4_decompress_safe(pointer(buffer), pointer(dec_buffer), result, test_size)
        @test unsafe_string(pointer(dec_buffer), test_size) == test_in

        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_compress_HC(pointer(test_in), pointer(buffer), test_size, half_bufsize, 12)
    end

    @testset "streams" begin
        comp_str = CodecLz4.LZ4_createStreamHC()
        try
            @test CodecLz4.LZ4_resetStreamHC(comp_str, 1) === nothing

            result = CodecLz4.LZ4_compress_HC_continue(comp_str, pointer(test_in), pointer(buffer), test_size, bufsize)
            @test result > 0
            result = CodecLz4.LZ4_decompress_safe(pointer(buffer), pointer(dec_buffer), result, test_size)

            @test result == test_size
            @test unsafe_string(pointer(dec_buffer), result) == test_in
        finally
            @test CodecLz4.LZ4_freeStreamHC(comp_str) == 0
        end
    end

    @testset "errors" begin
        stream = Ref{Ptr{CodecLz4.LZ4_streamHC_t}}(C_NULL)

        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_resetStreamHC(stream[])
        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_compress_HC_continue(stream[], pointer(test_in), pointer(buffer), test_size, bufsize)
        @test CodecLz4.LZ4_freeStreamHC(stream[]) == 0
    end
end


text:
using Wigxjpf
using Test

@testset "Wigxjpf.jl" begin
    Wigxjpf.wig_table_init(400, 9)
    Wigxjpf.wig_temp_init(400)

    @testset "Evaluate Wigner 3-j symbols" begin
        @test begin
            Wigxjpf.wig3jj(3, 3, 2, 3, -1, -2) ≈ -0.316227766016838
        end
        @test begin
            Wigxjpf.wig3jj(2, 4, 8, 2, -4, 2) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig3jj(2, 4, 2, 2, 2, -4) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig3jj(2, 4, 6, 1, -1, 0) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig3jj(2, 2, 2, 2, 0, -2) ≈ -0.408248290463863
        end
        @test begin
            Wigxjpf.wig3jj(6, 6, 6, -2, -2, 4) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig3jj(7, 5, 4, 7, -3, -4) ≈ -0.235702260395516
        end
        @test begin
            Wigxjpf.wig3jj(15, 15, 0, 3, -3, 0) ≈ 0.25
        end
        @test begin
            Wigxjpf.wig3jj(16, 11, 9, 4, -7, 3) ≈ -0.120622637445762
        end
        @test begin
            Wigxjpf.wig3jj(16, 14, 12, 6, 0, -6) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig3jj(30, 60, 80, 4, 4, -8) ≈ -0.0190815797991916
        end
        @test begin
            Wigxjpf.wig3jj(60, 60, 60, 0, 30, -30) ≈ -0.0152857154897839
        end
        @test begin
            Wigxjpf.wig3jj(286, 200, 120, -20, 120, -100) ≈ 4.16710792865494e-10
        end
        @test begin
            Wigxjpf.wig3jj(320, 200, 120, -20, 120, -100) ≈ 3.81124616116626e-21
        end
        @test begin
            Wigxjpf.wig3jj(400, 400, 400, -20, 120, -100) ≈ 0.000749392731398951
        end
        @test begin
            Wigxjpf.wig3jj(140, 150, 160, 40, -80, 40) ≈ -0.00845794815539587
        end
    end

    @testset "Evaluate Wigner 6-j symbols" begin
        @test begin
    Wigxjpf.wig6jj(4, 4, 2, 4, 2, 2) ≈ -0.1
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 0, 4, 4, 2) ≈ 0.258198889747161
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 2, 4, 4, 2) ≈ -0.223606797749979
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 4, 4, 4, 2) ≈ 0.152752523165195
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 0, 4, 4, 4) ≈ -0.258198889747161
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 0, 4, 4, 2) ≈ 0.258198889747161
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 2, 4, 4, 4) ≈ 0.074535599249993
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 2, 4, 4, 2) ≈ -0.223606797749979
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 4, 4, 4, 4) ≈ 0.152752523165195
        end
        @test begin
            Wigxjpf.wig6jj(2, 2, 4, 4, 4, 2) ≈ 0.152752523165195
        end
        @test begin
            Wigxjpf.wig6jj(4, 4, 2, 4, 2, 4) ≈ 0.152752523165195
        end
        @test begin
            Wigxjpf.wig6jj(7, 6, 3, 2, 3, 6) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig6jj(8, 8, 2, 8, 6, 2) ≈ -0.0277777777777778
        end
            @test begin
            Wigxjpf.wig6jj(12, 12, 8, 9, 7, 11) ≈ 0.0378232878232878
        end
        @test begin
            Wigxjpf.wig6jj(13, 12, 3, 6, 7, 12) ≈ 0.0
        end
        @test begin
            Wigxjpf.wig6jj(16, 16, 16, 16, 14, 14) ≈ 0.020221595293569
        end
        @test begin
            Wigxjpf.wig6jj(16, 16, 16, 16, 16, 16) ≈ -0.0126520807231535
        end
        @test begin
            Wigxjpf.wig6jj(40, 40, 40, 40, 40, 40) ≈ -0.00502940645686796
        end
        @test begin
            Wigxjpf.wig6jj(256, 240, 144, 224, 96, 160) ≈ 1.19377022891155e-17
            end
        @test begin
            Wigxjpf.wig6jj(460, 160, 300, 380, 460, 240) ≈ 2.42701001304468e-24
        end
        @test begin
            Wigxjpf.wig6jj(400, 400, 400, 400, 400, 400) ≈ 0.000155903212413242
        end
        @test begin
            Wigxjpf.wig6jj(160, 140, 120, 120, 140, 160) ≈ 0.000271672365647525
        end
    end
    @testset "Evaluate Wigner 9-j symbols" begin
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 2, 4, 6, 0, 4, 4) ≈ 0.03614031611621
        end
        @test begin
            Wigxjpf.wig9jj(1, 3, 4, 3, 1, 4, 2, 2, 0) ≈ 0.0645497224367903
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 4, 6, 10) ≈ -0.00141609931936472
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 4, 8, 10) ≈ 0.00793146728313409
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 6, 4, 10) ≈ 0.00759958641837765
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 6, 6, 10) ≈ 0.00851133416822978
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 8, 2, 10) ≈ 0.0108588135723727
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 8, 4, 10) ≈ 0.00335110262166189
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 8, 6, 10) ≈ -0.0050301051199454
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 6, 4, 8, 8, 8, 10) ≈ 0.00012246297616039
        end
        @test begin
            Wigxjpf.wig9jj(2, 4, 6, 4, 6, 8, 6, 8, 10) ≈ -0.000922146636432351
        end
        @test begin
            Wigxjpf.wig9jj(2, 6, 4, 3, 5, 4, 1, 1, 0) ≈ -0.0281718084909506
        end
        @test begin
            Wigxjpf.wig9jj(6, 8, 4, 7, 7, 4, 1, 1, 2) ≈ 0.0230021853114118
        end
        @test begin
            Wigxjpf.wig9jj(2, 6, 8, 1, 7, 6, 1, 1, 2) ≈ -0.00595238095238095
        end
        @test begin
            Wigxjpf.wig9jj(6, 7, 7, 5, 6, 7, 1, 1, 2) ≈ -0.00215400852516591
        end
        @test begin
            Wigxjpf.wig9jj(17, 19, 14, 25, 16, 17, 16, 21, 19) ≈
            0.000281298301912545
        end
        @test begin
        Wigxjpf.wig9jj(30, 30, 30, 30, 6, 30, 30, 36, 20) ≈
            -7.78324615309539e-05
        end
        @test begin
        Wigxjpf.wig9jj(90, 60, 40, 90, 30, 70, 180, 90, 90) ≈
            -6.51130967998823e-07
        end
        @test begin
            Wigxjpf.wig9jj(120, 140, 260, 100, 140, 240, 120, 100, 80) ≈
            5.35253924139905e-16
        end
        @test begin
            Wigxjpf.wig9jj(200, 160, 100, 100, 200, 140, 120, 100, 200) ≈
            1.05597798065761e-07
        end
        @test begin
            Wigxjpf.wig9jj(10, 1, 9, 10, 1, 11, 18, 2, 20) ≈ 0.00257786760370727
        end
        @test begin
            Wigxjpf.wig9jj(30, 30, 60, 30, 6, 30, 30, 36, 60) ≈
            0.000168384777511567
        end
        @test begin
            Wigxjpf.wig9jj(40, 20, 60, 60, 60, 120, 20, 40, 60) ≈
            0.000268744961031981
        end
        @test begin
            Wigxjpf.wig9jj(60, 40, 20, 60, 20, 40, 120, 60, 60) ≈
            0.000268744961031981
        end
        @test begin
            Wigxjpf.wig9jj(90, 60, 40, 90, 30, 70, 180, 90, 90) ≈
            -6.51130967998823e-07
        end
        @test begin
            Wigxjpf.wig9jj(100, 102, 104, 106, 108, 110, 112, 114, 116) ≈
            4.18840599233205e-06
        end
        @test begin
            Wigxjpf.wig9jj(148, 134, 120, 164, 120, 140, 20, 40, 60) ≈
            4.94488048916531e-06
        end
        @test begin
            Wigxjpf.wig9jj(168, 180, 122, 120, 128, 80, 104, 80, 60) ≈
            -7.94917060961424e-08
        end
        @test begin
            Wigxjpf.wig9jj(188, 134, 172, 122, 146, 104, 140, 120, 160) ≈
            -1.86808574631986e-07
        end
        @test begin
            Wigxjpf.wig9jj(120, 140, 230, 100, 140, 220, 120, 100, 80) ≈
            -1.36799650122649e-06
        end
        @test begin
            Wigxjpf.wig9jj(200, 160, 100, 100, 200, 140, 120, 100, 200) ≈
            1.05597798065761e-07
        end
        @test begin
            Wigxjpf.wig9jj(40, 40, 80, 40, 40, 80, 40, 40, 80) ≈
            0.000111633558388972
        end
    end

    Wigxjpf.wig_temp_free()
    Wigxjpf.wig_table_free()
end


text:
__precompile__()

module DiffEqWebBase
  using OrdinaryDiffEq, DiffEqBase

  import DiffEqBase: build_solution

  immutable QuickODEProblem{uType,tType,isinplace} <: AbstractODEProblem{uType,tType,isinplace}
    f::Function
    u0::uType
    tspan::Tuple{tType,tType}
    callback::C
  end

  function QuickODEProblem(f,u0,tspan)
    callback = CallbackSet()
    QuickODEProblem{typeof(u0),eltype(tspan),true,typeof(callback)}(f,u0,tspan,callback)
  end

  immutable QuickODESolution{uType,tType,rateType,P,A} <: AbstractODESolution
    u::uType
    t::tType
    k::rateType
    prob::P
    alg::A
    interp::OrdinaryDiffEq.InterpolationData
    dense::Bool
    tslocation::Int
  end
  (sol::QuickODESolution)(t) = sol.interp(t)

  function build_solution{uType,tType,isinplace}(
          prob::QuickODEProblem{uType,tType,isinplace},
          alg,t,u;dense=false,
          k=[],interp = (tvals) -> nothing,kwargs...)
    QuickODESolution(u,t,k,prob,alg,interp,dense,0)
  end

  type QuickSDEProblem{uType,tType,isinplace,NoiseClass,F3,C,ND} <: AbstractSDEProblem{uType,tType,isinplace,NoiseClass,ND}
    f::Function
    g::Function
    u0::uType
    tspan::Tuple{tType,tType}
    noise::DiffEqBase.NoiseProcess{NoiseClass,F3}
    callback::C
    noise_rate_prototype::ND
  end

  function QuickSDEProblem(f,g,u0,tspan)
    noise = INPLACE_WHITE_NOISE
    callback = CallbackSet()
    noise_rate_prototype = nothing
    QuickSDEProblem{typeof(u0),eltype(tspan),true,true,:White,typeof(noise.noise_func),typeof(callback),typeof(noise_rate_prototype)}(f,g,u0,tspan,noise,callback,noise_rate_prototype)
  end

  export QuickODEProblem, QuickODESolution,QuickSDEProblem
end


text:
# Test of vbroadcast

using NumericExtensions
using Base.Test

#### 2D ####

x = randn(5, 6)
yr = rand(6) + 0.5
yrm = reshape(yr, 1, 6)
yc = rand(5) + 0.5

# badd

@test_approx_eq badd(x, yc, 1) x .+ yc
@test_approx_eq badd(x, yr, 2) x .+ yrm

r = copy(x); badd!(r, yc, 1)
@test_approx_eq r x .+ yc

r = copy(x); badd!(r, yr, 2)
@test_approx_eq r x .+ yrm

# bsubtract

@test_approx_eq bsubtract(x, yc, 1) x .- yc
@test_approx_eq bsubtract(x, yr, 2) x .- yrm

r = copy(x); bsubtract!(r, yc, 1)
@test_approx_eq r x .- yc

r = copy(x); bsubtract!(r, yr, 2)
@test_approx_eq r x .- yrm

# bmultiply

@test_approx_eq bmultiply(x, yc, 1) x .* yc
@test_approx_eq bmultiply(x, yr, 2) x .* yrm

r = copy(x); bmultiply!(r, yc, 1)
@test_approx_eq r x .* yc

r = copy(x); bmultiply!(r, yr, 2)
@test_approx_eq r x .* yrm

# bdivide

@test_approx_eq bdivide(x, yc, 1) x ./ yc
@test_approx_eq bdivide(x, yr, 2) x ./ yrm

r = copy(x); bdivide!(r, yc, 1)
@test_approx_eq r x ./ yc

r = copy(x); bdivide!(r, yr, 2)
@test_approx_eq r x ./ yrm

# #### 3D ####

x = randn(4, 5, 6)

y1 = rand(4) + 0.5; y1m = reshape(y1, 4, 1, 1)
y2 = rand(5) + 0.5; y2m = reshape(y2, 1, 5, 1)
y3 = rand(6) + 0.5; y3m = reshape(y3, 1, 1, 6)

# badd

@test_approx_eq badd(x, y1, 1) x .+ y1m
@test_approx_eq badd(x, y2, 2) x .+ y2m
@test_approx_eq badd(x, y3, 3) x .+ y3m

# bsubtract

@test_approx_eq bsubtract(x, y1, 1) x .- y1m
@test_approx_eq bsubtract(x, y2, 2) x .- y2m
@test_approx_eq bsubtract(x, y3, 3) x .- y3m

# bmultiply

@test_approx_eq bmultiply(x, y1, 1) x .* y1m
@test_approx_eq bmultiply(x, y2, 2) x .* y2m
@test_approx_eq bmultiply(x, y3, 3) x .* y3m

# bdivide

@test_approx_eq bdivide(x, y1, 1) x ./ y1m
@test_approx_eq bdivide(x, y2, 2) x ./ y2m
@test_approx_eq bdivide(x, y3, 3) x ./ y3m


text:
"""

Support for out-of-core sparse arrays backed by a HDF5 file stored on disk. Provides `H5SparseMatrixCSC`.
"""
module H5Sparse

using HDF5, SparseArrays

export H5SparseMatrixCSC

"""
    H5SparseMatrixCSC{Tv, Ti<:Integer, Td<:HDF5.H5DataStore} <: SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}

Out-of-core `AbstractSparseMatrixCSC` backed by a dataset stored on disk, of type `Td<:HDF5.H5DataStore`, e.g., a `HDF5.File`.

```julia
# Conversion from SparseMatrixCSC; writes B to a dataset "A" in the file "foo.h5"
using SparseArrays
B = sprand(10, 10, 0.5)
A = H5SparseMatrixCSC("foo.h5", "A", B)

# or, equivalently
using HDF5
fid = h5open("foo.h5", "cw")
A = H5SparseMatrixCSC(fid, "A", B)

# kwargs are passed on to h5writecsc
A = H5SparseMatrixCSC("foo.h5", "A", B, overwrite=true) # Overwrites any existing dataset with name A

# Construct from an existing file
A = H5SparseMatrixCSC("foo.h5", "A")
A = H5SparseMatrixCSC(fid, "A")

# Construct a view into a subset of the rows and/or columns stored in a file
A = H5SparseMatrixCSC("foo.h5", "A", :, 2:5)

# Colon or UnitRange indexing returns a new H5SparseMatrixCSC that is a view into the specified subset of rows and/or columns
A[:, 1:10]
A[1:4, :]

# Integer indexing returns the requested element
A[1, 1]

# Concatenate with a SparseMatrixCSC to the right; useful for constructing large matrices in an iterative fashion
# Returns a new H5SparseMatrixCSC spanning all columns of the resulting matrix
C = sprand(10, 5, 0.5)
A = hcat(A, C)      # A is now of size (10, 15)

# Load a H5SparseMatrixCSC matrix into memory using SparseArrays.sparse
# Only the columns covered by this particular H5SparseMatrixCSC will be read from disk
# However, A must cover all rows of the underlying matrix
sparse(A)           # SparseMatrixCSC
sparse(A[:, 1:4])   # The first 4 columns of A as a new SparseMatrixCSC
sparse(A[1:4, :])   # Results in an error (not implemented)

# If converting to a dense matrix, Matrix(sparse(A)) is likely orders of magnitude faster than calling Matrix(A) directly
# (since Matrix(A) does not take advantage of A being sparse)
Matrix(A)           # Matrix (slow)
Matrix(sparse(A))   # Matrix (fast)
```

"""

struct H5SparseMatrixCSC{Tv, Ti<:Integer, Td<:HDF5.H5DataStore} <: SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}
    fid::Td                 # Backing storage
    name::String            # Dataset name, i.e., data is stored in fid[name]
    rows::UnitRange{Int}    # Subset of rows stored in fid[name] accessible via this instance
    cols::UnitRange{Int}    # Subset of columns stored in fid[name] accessible via this instance
    function H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString, rows::UnitRange{Int}, cols::UnitRange{Int}) where {Tv,Ti<:Integer}
        name in keys(fid) || throw(ArgumentError("$name is not in $fid"))
        g = fid[name]
        g isa HDF5.Group || throw(ArgumentError("fid[name] is $g, but must be a HDF5.Group"))
        "m" in keys(g) || throw(ArgumentError("m is not in $g"))
        "n" in keys(g) || throw(ArgumentError("n is not in $g"))
        "colptr" in keys(g) || throw(ArgumentError("colptr is not in $g"))
        "rowval" in keys(g) || throw(ArgumentError("rowval is not in $g"))
        "nzval" in keys(g) || throw(ArgumentError("nzval is not in $g"))
        length(size(g["colptr"])) == 1 || throw(ArgumentError("expected colptr to have 1 dimension, but got $(size(g["colptr"]))"))
        length(size(g["rowval"])) == 1 || throw(ArgumentError("expected rowval to have 1 dimension, but got $(size(g["rowval"]))"))
        length(size(g["nzval"])) == 1 || throw(ArgumentError("expected nzval to have 1 dimension, but got $(size(g["nzval"]))"))
        eltype(g["colptr"]) == eltype(g["rowval"]) || throw(ArgumentError("colptr has eltype $(g["colptr"])), but rowval has eltype $(g["rowval"]))"))
        m, n = g["m"][], g["n"][]
        (iszero(m) && first(rows) == 0) || 0 < first(rows) || throw(ArgumentError("first row is $(first(rows)), but m is $m"))
        (iszero(m) && last(rows) == 0)  || 0 < last(rows) || throw(ArgumentError("last row is $(last(rows)), but m is $m"))
        first(rows) <= m || throw(ArgumentError("first row is $(first(rows)), but m is $m"))
        last(rows) <= m || throw(ArgumentError("last row is $(last(rows)), but m is $m"))
        first(rows) <= last(rows) || throw(ArgumentError("first row is $(first(rows)), but last row is $(last(rows))"))
        (iszero(n) && first(cols) == 0) || 0 < first(cols) || throw(ArgumentError("first column is $(first(cols)), but n is $n"))
        (iszero(n) && last(cols) == 0)  || 0 < last(cols) || throw(ArgumentError("last column is $(last(cols)), but n is $n"))
        first(cols) <= n || throw(ArgumentError("first column is $(first(cols)), but n is $n"))
        last(cols) <= n || throw(ArgumentError("last column is $(last(cols)), but n is $n"))
        first(cols) <= last(cols) || throw(ArgumentError("first columns is $(first(cols)), but last column is $(last(cols))"))
        new{eltype(g["nzval"]),eltype(g["rowval"]),typeof(fid)}(fid, String(name), rows, cols)
    end
end
H5SparseMatrixCSC(filename::AbstractString, args...; kwargs...) = H5SparseMatrixCSC(h5open(filename, "cw"), args...; kwargs...)
function H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString, B::SparseMatrixCSC; kwargs...)
    h5writecsc(fid, name, B; kwargs...)
    H5SparseMatrixCSC(
        fid, name,
        size(B, 1) > 0 ? (1:size(B, 1)) : (0:0),
        size(B, 2) > 0 ? (1:size(B, 2)) : (0:0),
        )
end
function H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString)
    m, n = h5size(fid, name)
    H5SparseMatrixCSC(fid, name, 1:m, 1:n)
end
function H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString, ::Colon, cols)
    m, n = h5size(fid, name)
    H5SparseMatrixCSC(fid, name, 1:m, cols)
end
function H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString, rows, ::Colon)
    m, n = h5size(fid, name)
    H5SparseMatrixCSC(fid, name, rows, 1:n)
end
H5SparseMatrixCSC(fid::HDF5.H5DataStore, name::AbstractString, ::Colon, ::Colon) = H5SparseMatrixCSC(fid, name)

Base.show(io::IOContext, A::H5SparseMatrixCSC) = show(io.io, A)

function Base.show(io::IO, A::H5SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    print(io, "H5SparseMatrixCSC{$Tv,$Ti}($(A.fid.filename), $(A.name), $(A.rows), $(A.cols))")
end

function Base.display(A::H5SparseMatrixCSC)
    show(stdout, A)
end

Base.size(A::H5SparseMatrixCSC) = (last(A.rows)-first(A.rows)+1, last(A.cols)-first(A.cols)+1)
function Base.eltype(::H5SparseMatrixCSC{Tv}) where Tv
    Tv
end

struct H5View{Tv} <: AbstractVector{Tv}
    x::HDF5.Dataset # Data
    i::Int          # Lower index
    j::Int          # Upper index
end
H5View(x::HDF5.Dataset, i::Integer, j::Integer) = H5View{eltype(x)}(x, i, j)
H5View(x::HDF5.Dataset) = H5View(x, 1, length(x))

function Base.show(io::IO, a::H5View{Tv}) where Tv
    print(io, "H5View{$Tv}($(a.x), $(x.i), $(x.j))")
end

function Base.display(a::H5View)
    show(stdout, a)
end

Base.length(a::H5View) = a.j-a.i+1
Base.size(a::H5View) = (length(a),)
function Base.eltype(::H5View{Tv}) where Tv
    Tv
end
function Base.getindex(a::H5View{Tv}, k::Integer)::Tv where Tv
    @boundscheck 0 < k <= length(a)
    a.x[a.i+k-1]
end

SparseArrays.getcolptr(A::H5SparseMatrixCSC) = H5View(A.fid[A.name]["colptr"])
SparseArrays.rowvals(A::H5SparseMatrixCSC) = H5View(A.fid[A.name]["rowval"])
SparseArrays.nonzeros(A::H5SparseMatrixCSC) = H5View(A.fid[A.name]["nzval"])

function Base.getindex(A::H5SparseMatrixCSC{Tv}, row::Integer, col::Integer)::Tv where Tv
    m, n = size(A)
    @boundscheck 0 < row <= m || throw(BoundsError(A, (row, col)))
    @boundscheck 0 < col <= n || throw(BoundsError(A, (row, col)))
    row = first(A.rows) + row - 1
    col = first(A.cols) + col - 1
    g = A.fid[A.name]
    i = g["colptr"][col]
    j = g["colptr"][col+1] - 1
    v = H5View(g["rowval"], i, j)
    k = searchsortedfirst(v, row) + i - 1
    if k > j || g["rowval"][k] != row
        return zero(Tv)
    else
        return g["nzval"][k]
    end
end

Base.getindex(A::H5SparseMatrixCSC, ::Colon, cols::UnitRange{<:Integer}) = getindex(A, 1:size(A, 1), cols)
Base.getindex(A::H5SparseMatrixCSC, rows::UnitRange{<:Integer}, ::Colon) = getindex(A, rows, 1:size(A, 2))
Base.getindex(A::H5SparseMatrixCSC, ::Colon, ::Colon) = getindex(A, 1:size(A, 1), 1:size(A, 2))
function Base.getindex(A::H5SparseMatrixCSC{Tv,Ti}, rows::UnitRange{<:Integer}, cols::UnitRange{<:Integer})::H5SparseMatrixCSC{Tv,Ti} where {Tv,Ti}
    m, n = size(A)
    @boundscheck 0 < first(rows) <= m || throw(BoundsError(A, (rows, cols)))
    @boundscheck 0 < last(rows) <= m || throw(BoundsError(A, (rows, cols)))
    @boundscheck first(rows) <= last(rows) || throw(ArgumentError("first row is $(first(rows)), but last row is $(last(rows))"))
    @boundscheck 0 < first(cols) <= n || throw(BoundsError(A, (rows, cols)))
    @boundscheck 0 < last(cols) <= n || throw(BoundsError(A, (rows, cols)))
    @boundscheck first(cols) <= last(cols) || throw(ArgumentError("first column is $(first(cols)), but last column is $(last(cols))"))
    H5SparseMatrixCSC(A.fid, A.name, (first(A.rows)+first(rows)-1):(first(A.rows)+last(rows)-1), (first(A.cols)+first(cols)-1):(first(A.cols)+last(cols)-1))
end

function Base.setindex!(::H5SparseMatrixCSC, ::Any, ::Integer, ::Integer)
    error("H5SparseMatrixCSC data is immutable")
end

"""
    append!(A::H5SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

Appends `B` to the right of `A`.
"""
function Base.hcat(A::H5SparseMatrixCSC{Tv,Ti}, Bs::SparseMatrixCSC{Tv,Ti}...) where {Tv,Ti}
    ncols = 0
    for B in Bs
        size(A, 1) == size(B, 1) || throw(DimensionMismatch("A has dimensions $(size(A)), but B has dimensions $(size(B))"))
        h5appendcsc(A.fid, A.name, B)
        ncols += size(B, 2)
    end
    H5SparseMatrixCSC(A.fid, A.name, A.rows, first(A.cols):(last(A.cols)+ncols))
end

function SparseArrays.sparse(A::H5SparseMatrixCSC{Tv,Ti})::SparseMatrixCSC{Tv,Ti} where {Tv,Ti}
    m, n = size(A)
    if m != h5size(A.fid, A.name)[1]
        error("not implemented")
    end
    g = A.fid[A.name]
    colptr::Vector{Ti} = g["colptr"][first(A.cols):(last(A.cols)+1)]
    i = colptr[1]
    j = colptr[end] - 1
    rowval::Vector{Ti} = j >= i ? g["rowval"][i:j] : Vector{Ti}()
    nzval::Vector{Tv} = j >= i ? g["nzval"][i:j] : Vector{Tv}()
    colptr .-= i-1
    SparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval)
end

function h5size(fid::HDF5.H5DataStore, name::AbstractString)
    g = fid[name]
    g["m"][], g["n"][]
end

"""
    h5writecsc(fid::HDF5.H5DataStore, name::AbstractString, B::SparseMatrixCSC; kwargs...)

Write `B` to `fid[name]`, overwriting any existing dataset with name `name` if `overwrite=true`.
This is a low-level routine without error checking; user should typically use `H5SparseMatrixCSC`
instead.
"""
function h5writecsc(fid::HDF5.H5DataStore, name::AbstractString, B::SparseMatrixCSC; kwargs...)
    m, n = size(B)
    h5writecsc(fid, name, m, n, SparseArrays.getcolptr(B), rowvals(B), nonzeros(B); kwargs...)
end
function h5writecsc(fid::HDF5.H5DataStore, name::AbstractString, m::Integer, n::Integer, colptr::AbstractVector{<:Integer}, rowval::AbstractVector{<:Integer}, nzval::AbstractVector; overwrite=false, colptr_chunk=length(colptr) > 0 ? HDF5.heuristic_chunk(colptr) : [1], val_chunk=length(rowval) > 0 ? HDF5.heuristic_chunk(rowval) : [1], blosc=5, kwargs...)
    if name in keys(fid)
        if overwrite
            delete_object(fid, name)
        else
            throw(ArgumentError("$name already exists in $fid"))
        end
    end
    g = create_group(fid, name)
    g["m"] = m
    g["n"] = n
    ds = create_dataset(
        g, "colptr",
        eltype(colptr),
        ((length(colptr),), (-1,)),
        chunk=colptr_chunk,
        blosc=blosc,
        kwargs...,
    )
    if length(colptr) > 0
        ds[1:length(colptr)] = colptr
    end
    ds = create_dataset(
        g, "rowval",
        eltype(rowval),
        ((length(rowval),), (-1,)),
        chunk=val_chunk,
        blosc=blosc,
        kwargs...,
    )
    if length(rowval) > 0
        ds[1:length(rowval)] = rowval
    end
    ds = create_dataset(
        g, "nzval",
        eltype(nzval),
        ((length(nzval),), (-1,)),
        chunk=val_chunk,
        blosc=blosc,
        kwargs...,
    )
    if length(nzval) > 0
        ds[1:length(nzval)] = nzval
    end
    return
end

"""
    h5appendcsc(fid::HDF5.H5DataStore, name::AbstractString, B::SparseMatrixCSC)

Append `B` to the right of the `SparseMatrixCSC` stored in `fid[name]`. This is a low-level routine
without error checking; user should typically use `append` instead.
"""
function h5appendcsc(fid::HDF5.H5DataStore, name::AbstractString, B::SparseMatrixCSC)
    m, n = size(B)
    h5appendcsc(fid, name, m, n, SparseArrays.getcolptr(B), rowvals(B), nonzeros(B))
end
function h5appendcsc(fid::HDF5.H5DataStore, name::AbstractString, m::Integer, n::Integer, colptr::AbstractVector{<:Integer}, rowval::AbstractVector{<:Integer}, nzval::AbstractVector)
    g = fid[name]
    # colptr
    i = size(g["colptr"], 1)
    j = i + length(colptr) - 1
    HDF5.set_extent_dims(g["colptr"], (j,))
    offset = size(g["rowval"], 1)
    g["colptr"][i:j]  = colptr .+ offset
    # rowval
    i = size(g["rowval"], 1) + 1
    j = i + length(rowval) - 1
    HDF5.set_extent_dims(g["rowval"], (j,))
    g["rowval"][i:j]  = rowval
    # nzval
    i = size(g["nzval"], 1) + 1
    j = i + length(nzval) - 1
    HDF5.set_extent_dims(g["nzval"], (j,))
    g["nzval"][i:j]  = nzval
    # n
    old_n = g["n"][]
    delete_object(g, "n")
    g["n"] = old_n + n
    return
end

"""

Return `true` if `fid[name]` is a valid `H5SparseMatrixCSC` dataset, and `false` otherwise.
"""
function h5isvalidcsc(fid::HDF5.H5DataStore, name::AbstractString)
    name in keys(fid) || return false
    g = fid[name]
    g isa HDF5.Group || return false
    "m" in keys(g) || return false
    "n" in keys(g) || return false
    "colptr" in keys(g) || return false
    "rowval" in keys(g) || return false
    "nzval" in keys(g) || return false
    length(size(g["colptr"])) == 1 || return false
    length(size(g["rowval"])) == 1 || return false
    length(size(g["nzval"])) == 1 || return false
    eltype(g["colptr"]) == eltype(g["rowval"]) || return false
    true
end

"""

Return `true` if `h5open(filename)[name]` is a valid `H5SparseMatrixCSC` dataset, and `false`
otherwise.
"""
function h5isvalidcsc(filename::AbstractString, name::AbstractString)
    HDF5.ishdf5(filename) || return false
    h5open(filename) do fid
        return h5isvalidcsc(fid, name)
    end
end

end


text:
import HiGHS

@testset "FEASIBILITY_SENSE zeros objective" begin
    model = DiffOpt.diff_optimizer(HiGHS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)

    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol=ATOL rtol=RTOL

    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
end

@testset "Forward or reverse without optimizing throws" begin
    model = DiffOpt.diff_optimizer(HiGHS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    # do not optimize, just try to differentiate
    @test_throws ErrorException DiffOpt.forward_differentiate!(model)
    @test_throws ErrorException DiffOpt.reverse_differentiate!(model)
    # impossible constraint
    MOI.add_constraint(model, x, MOI.LessThan(0.5))
    MOI.optimize!(model)
    @test_throws ErrorException DiffOpt.forward_differentiate!(model)
    @test_throws ErrorException DiffOpt.reverse_differentiate!(model)
end


text:
using Test, Cassette, LinearAlgebra
using Cassette: @context, @pass, @overdub, overdub, hasmetadata, metadata, hasmetameta,
                metameta, untag, tag, enabletagging, untagtype, istagged, istaggedtype,
                Tagged, fallback, canoverdub, similarcontext

println("running unit tests")
@time @testset "unit tests" begin include("unittests.jl") end

println("running misc. tests (w/o tagging)")
@time @testset "misc. tests (w/o tagging)" begin include("misctests.jl") end

println("running misc. tests (w/ tagging)")
@time @testset "misc. tests (w/ tagging)" begin include("misctaggingtests.jl") end


text:
struct Common
    runid::String
    feeder::Space
    fictive_nest::Space
    track::Track
    pellets::Union{Missing, Vector{Space}}
    nest::Union{Missing, Space}
    pickup::Union{Missing, Space}
    dropoff::Space
end

Base.getproperty(x::Common, k::Symbol) = k === :homing ? x.track.homing :
                                         k === :searching ? x.track.searching :
                                         k === :center_of_search ? x.track.center_of_search :
                                         k === :turning_point ? x.track.turning_point :
                                         getfield(x, k)

using Unitful, Statistics, StaticArrays, Rotations, CoordinateTransformations, LinearAlgebra

function getmissingfeeder(guess, displacement, nest, dropoff, nest2feeder)
    c = calculatec(guess, displacement, nest, dropoff)
    u = LinearAlgebra.normalize(nest - c)
    return nest + nest2feeder*u
end

getfeeder(x, metadata, nest, dropoff, nest2feeder) = haskey(x, :initialfeeder) ? x[:initialfeeder] :
                                           haskey(x, :feeder) ? x[:feeder] :
                                           getmissingfeeder(x[:guess], metadata[:expected_dropoff], nest, dropoff, nest2feeder)

getpickup(data, ::Missing, _) = missing
getpickup(data, _, feeder) = haskey(data, :rightdowninitial) ? mean(data[k] for k in (:rightdowninitial, :leftdowninitial, :rightupinitial, :leftupinitial)) :
                             haskey(data, :pickup) ? data[:pickup] :
                             get(data, :initialfeeder, feeder)

# getdropoff(data, ::Nothing) = nothing
getdropoff(data) = haskey(data, :rightdownfinal) ? mean(data[k] for k in (:rightdownfinal, :leftdownfinal, :rightupfinal, :leftupfinal)) :
                      haskey(data, :dropoff) ? data[:dropoff] :
                      data[:feeder]

getnest2feeder(x, metadata) = haskey(x, :nestbefore) ? norm(x[:nestbefore] - x[:feederbefore]) :
                    get(metadata, :nest2feeder, missing)

# function ayse_fix!(x)
#     if !haskey(x, :feeder) && !haskey(x, :guess)
#         x[:feeder] = first(x[:track]).xy
#     end
# end

function getdata(x, metadata)
    # ayse_fix!(x)
    nest = get(x, :nest, missing)
    track = Track(x[:track], metadata[:turning_point])
    _pellet = get(x, :pellet, StructVector{SpaceTime}(undef, 0))
    pellet = _pellet.xy
    dropoff = getdropoff(x)
    nest2feeder = getnest2feeder(x, metadata)
    feeder = getfeeder(x, metadata, nest, dropoff, nest2feeder)
    pickup = getpickup(x, nest, feeder)
    feeder, nest, track, pellet, pickup, dropoff, nest2feeder
end

function common(x, metadata)
    feeder, nest, track, pellet, pickup, dropoff, nest2feeder = getdata(x, metadata)
    fictive_nest = getfictive_nest(x, metadata, pickup, nest, dropoff, nest2feeder)
    Common(metadata[:runid], feeder, fictive_nest, track, pellet, nest, pickup, dropoff)
end

getfictive_nest(x, metadata, pickup::Missing, nest::Space, dropoff::Missing, _) = nest
getfictive_nest(x, metadata, pickup::Space, nest::Space, dropoff::Space, _) = nest + dropoff - pickup
function getfictive_nest(x, metadata, pickup::Missing, nest::Missing, dropoff::Space, nest2feeder)
    south = x[:south]
    north = x[:north]
    v = north - south
    azimuth = getazimuth(x, metadata)
    α = atan(v[2], v[1]) - azimuth# - π
    u = Space(cos(α), sin(α))
    return dropoff + u*nest2feeder
end

function calculatec(guess, displacement, nest, dropoff)
    v = dropoff - nest
    α = π/2 - atan(v[2], v[1])
    t = LinearMap(Angle2d(α)) ∘ Translation(-nest)
    ṫ = inv(t)
    ab = norm(displacement)
    bc², ac² = displacement.^2
    cy = (ab^2 + ac² - bc²)/2ab
    Δ = sqrt(ac² - cy^2)
    c = [ṫ([i*Δ, cy]) for i in (-1, 1)]
    l = norm.(c .- Ref(guess))
    _, i = findmin(l)
    c[i]
end

function anglebetween(north, south, nest, feeder)
    v = LinearAlgebra.normalize(north - south)
    u = LinearAlgebra.normalize(nest - feeder)
    acos(v ⋅ u)
    # atan(v[2], v[1]) - atan(u[2], u[1])
end
getazimuth(x, metadata) = haskey(x, :southbefore) ? anglebetween(x[:northbefore], x[:southbefore], x[:nestbefore], x[:feederbefore]) :
                haskey(metadata, :azimuth) ?  deg2rad(metadata[:azimuth]) :
                missing

######################### END ######################

struct Standardized
    runid::String
    nest2feeder::Float64
    fictive_nest::Space
    track::Track
    pellets::Union{Missing, Vector{Space}}
    pickup::Union{Missing, Space}
    dropoff::Space
    nest::Union{Missing, Space}
end

Base.getproperty(x::Standardized, k::Symbol) = k === :homing ? x.track.homing :
                                         k === :searching ? x.track.searching :
                                         k === :center_of_search ? x.track.center_of_search :
                                         k === :turning_point ? x.track.turning_point :
                                         k == :feeder ? Space(0.0, -x.nest2feeder) :
                                         getfield(x, k)

function Standardized(x)
    t = get_transform(x.nest, x.fictive_nest, x.feeder)
    fictive_nest = t(x.fictive_nest)
    map!(t, x.track.coords, x.track.coords)
    for i in eachindex(x.track.rawcoords)
      xyt = x.track.rawcoords[i]
      x.track.rawcoords[i] = SpaceTime(t(xyt.xy), xyt.t)
    end
    map!(t, x.pellets, x.pellets)
    pickup = t(x.pickup)
    dropoff = t(x.dropoff)
    nest2feeder = abs(last(t(x.feeder)))
    Standardized(x.runid, nest2feeder, fictive_nest, x.track, x.pellets, pickup, dropoff, t(x.nest))
end

get_center(::Missing, fictive_nest) = fictive_nest
get_center(nest, _) = nest
function get_rotation(v)
    α = atan(v[2], v[1]) + π/2
    LinearMap(Angle2d(-α))
end

function get_transform(nest, fictive_nest, feeder)
    c = get_center(nest, fictive_nest)
    trans = Translation(-c)
    rot = get_rotation(trans(feeder))
    passmissing(rot ∘ trans)
end


text:
using Revise
using NNFEM

init_gmsh()
pts = [
    addPoint(0,0, 1.0),
    addPoint(8,0, 1.0),
    addPoint(8,4, 1.0),
    addPoint(0,4, 1.0),
    addPoint(3.13, 0.0, 0.8),
    addPoint(4.0, 0.5, 0.8)
]
rectangle = addPlaneSurface(
    [addCurveLoop([
        addLine(1,2),
        addLine(2,3),
        addLine(3,4),
        addLine(4,1)
    ])]
)
# meshsize("")
line = addLine(5,6)
embedLine([line], rectangle)
filename = finalize_gmsh()
cp(filename, "test.msh")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
# Residual block from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export ResidualBlock, ResidualBlock3D

"""
    RB = ResidualBlock(n_in, n_hidden; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false)
    RB = ResidualBlock3D(n_in, n_hidden; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false)

or

    RB = ResidualBlock(W1, W2, W3, b1, b2; p1=1, p2=1, s1=1, s2=1, fan=false)
    RB = ResidualBlock3D(W1, W2, W3, b1, b2; p1=1, p2=1, s1=1, s2=1, fan=false)

 Create a (non-invertible) residual block, consisting of three convolutional layers and activation functions.
 The first convolution is a downsampling operation with a stride equal to the kernel dimension. The last
 convolution is the corresponding transpose operation and upsamples the data to either its original dimensions
 or to twice the number of input channels (for `fan=true`). The first and second layer contain a bias term.

 *Input*:

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 - `fan`: bool to indicate whether the ouput has twice the number of input channels. For `fan=false`, the last
    activation function is a gated linear unit (thereby bringing the output back to the original dimensions).
    For `fan=true`, the last activation is a ReLU, in which case the output has twice the number of channels
    as the input.

or

 - `W1`, `W2`, `W3`: 4D tensors of convolutional weights

 - `b1`, `b2`: bias terms

 *Output*:

 - `RB`: residual block layer

 *Usage:*

 - Forward mode: `Y = RB.forward(X)`

 - Backward mode: `ΔX = RB.backward(ΔY, X)`

 *Trainable parameters:*

 - Convolutional kernel weights `RB.W1`, `RB.W2` and `RB.W3`

 - Bias terms `RB.b1` and `RB.b2`

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct ResidualBlock <: NeuralNetLayer
    W1::Parameter
    W2::Parameter
    W3::Parameter
    b1::Parameter
    b2::Parameter
    fan::Bool
    strides
    pad
end

@Flux.functor ResidualBlock

#######################################################################################################################
#  Constructors

# Constructor
function ResidualBlock(n_in, n_hidden; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false, ndims=2)

    k1 = Tuple(k1 for i=1:ndims)
    k2 = Tuple(k2 for i=1:ndims)
    # Initialize weights
    W1 = Parameter(glorot_uniform(k1..., n_in, n_hidden))
    W2 = Parameter(glorot_uniform(k2..., n_hidden, n_hidden))
    W3 = Parameter(glorot_uniform(k1..., 2*n_in, n_hidden))
    b1 = Parameter(zeros(Float32, n_hidden))
    b2 = Parameter(zeros(Float32, n_hidden))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, (s1, s2), (p1, p2))
end

# Constructor for given weights
function ResidualBlock(W1, W2, W3, b1, b2; p1=1, p2=1, s1=1, s2=1, fan=false, ndims=2)

    # Make weights parameters
    W1 = Parameter(W1)
    W2 = Parameter(W2)
    W3 = Parameter(W3)
    b1 = Parameter(b1)
    b2 = Parameter(b2)

    return ResidualBlock(W1, W2, W3, b1, b2, fan, (s1, s2), (p1, p2))
end

ResidualBlock3D(args...; kw...) = ResidualBlock(args...; kw..., ndims=3)
#######################################################################################################################
# Functions

# Forward
function forward(X1::AbstractArray{Float32, N}, RB::ResidualBlock; save=false) where {N}
    inds =[i!=(N-1) ? 1 : (:) for i=1:N]

    Y1 = conv(X1, RB.W1.data; stride=RB.strides[1], pad=RB.pad[1]) .+ reshape(RB.b1.data, inds...)
    X2 = ReLU(Y1)

    Y2 = X2 + conv(X2, RB.W2.data; stride=RB.strides[2], pad=RB.pad[2]) .+ reshape(RB.b2.data, inds...)
    X3 = ReLU(Y2)

    cdims3 = DCDims(X1, RB.W3.data; nc=2*size(X1, N-1), stride=RB.strides[1], padding=RB.pad[1])
    Y3 = ∇conv_data(X3, RB.W3.data, cdims3)
    RB.fan == true ? (X4 = ReLU(Y3)) : (X4 = GaLU(Y3))

    if save == false
        return X4
    else
        return Y1, Y2, Y3, X2, X3
    end
end

# Backward
function backward(ΔX4::AbstractArray{Float32, N}, X1::AbstractArray{Float32, N},
                  RB::ResidualBlock; set_grad::Bool=true) where {N}
    inds = [i!=(N-1) ? 1 : (:) for i=1:N]
    dims = collect(1:N-1); dims[end] +=1

    # Recompute forward states from input X
    Y1, Y2, Y3, X2, X3 = forward(X1, RB; save=true)

    # Cdims
    cdims2 = DenseConvDims(X2, RB.W2.data; stride=RB.strides[2], padding=RB.pad[2])
    cdims3 = DCDims(X1, RB.W3.data; nc=2*size(X1, N-1), stride=RB.strides[1], padding=RB.pad[1])

    # Backpropagate residual ΔX4 and compute gradients
    RB.fan == true ? (ΔY3 = ReLUgrad(ΔX4, Y3)) : (ΔY3 = GaLUgrad(ΔX4, Y3))
    ΔX3 = conv(ΔY3, RB.W3.data, cdims3)
    ΔW3 = ∇conv_filter(ΔY3, X3, cdims3)

    ΔY2 = ReLUgrad(ΔX3, Y2)
    ΔX2 = ∇conv_data(ΔY2, RB.W2.data, cdims2) + ΔY2
    ΔW2 = ∇conv_filter(X2, ΔY2, cdims2)
    Δb2 = sum(ΔY2, dims=dims)[inds...]

    cdims1 = DenseConvDims(X1, RB.W1.data; stride=RB.strides[1], padding=RB.pad[1])

    ΔY1 = ReLUgrad(ΔX2, Y1)
    ΔX1 = ∇conv_data(ΔY1, RB.W1.data, cdims1)
    ΔW1 = ∇conv_filter(X1, ΔY1, cdims1)
    Δb1 = sum(ΔY1, dims=dims)[inds...]

    # Set gradients
    if set_grad
        RB.W1.grad = ΔW1
        RB.W2.grad = ΔW2
        RB.W3.grad = ΔW3
        RB.b1.grad = Δb1
        RB.b2.grad = Δb2
    else
        Δθ = [Parameter(ΔW1), Parameter(ΔW2), Parameter(ΔW3), Parameter(Δb1), Parameter(Δb2)]
    end

    set_grad ? (return ΔX1) : (return ΔX1, Δθ)
end

## Jacobian-related functions
function jacobian(ΔX1::AbstractArray{Float32, N}, Δθ::Array{Parameter, 1},
                  X1::AbstractArray{Float32, N}, RB::ResidualBlock) where {N}
    inds = [i!=(N-1) ? 1 : (:) for i=1:N]
    # Cdims
    cdims1 = DenseConvDims(X1, RB.W1.data; stride=RB.strides[1], padding=RB.pad[1])

    Y1 = conv(X1, RB.W1.data, cdims1) .+ reshape(RB.b1.data, inds...)
    ΔY1 = conv(ΔX1, RB.W1.data, cdims1) + conv(X1, Δθ[1].data, cdims1) .+ reshape(Δθ[4].data, inds...)
    X2 = ReLU(Y1)
    ΔX2 = ReLUgrad(ΔY1, Y1)

    cdims2 = DenseConvDims(X2, RB.W2.data; stride=RB.strides[2], padding=RB.pad[2])

    Y2 = X2 + conv(X2, RB.W2.data, cdims2) .+ reshape(RB.b2.data, inds...)
    ΔY2 = ΔX2 + conv(ΔX2, RB.W2.data, cdims2) + conv(X2, Δθ[2].data, cdims2) .+ reshape(Δθ[5].data, inds...)
    X3 = ReLU(Y2)
    ΔX3 = ReLUgrad(ΔY2, Y2)

    cdims3 = DCDims(X1, RB.W3.data; nc=2*size(X1, N-1), stride=RB.strides[1], padding=RB.pad[1])
    Y3 = ∇conv_data(X3, RB.W3.data, cdims3)
    ΔY3 = ∇conv_data(ΔX3, RB.W3.data, cdims3) + ∇conv_data(X3, Δθ[3].data, cdims3)
    if RB.fan == true
        X4 = ReLU(Y3)
        ΔX4 = ReLUgrad(ΔY3, Y3)
    else
        ΔX4, X4 = GaLUjacobian(ΔY3, Y3)
    end

    return ΔX4, X4

end

# 2D/3D
function adjointJacobian(ΔX4::AbstractArray{Float32, N}, X1::AbstractArray{Float32, N}, RB::ResidualBlock) where N
    return backward(ΔX4, X1, RB; set_grad=false)
end

## Other utils
# Clear gradients
function clear_grad!(RB::ResidualBlock)
    RB.W1.grad = nothing
    RB.W2.grad = nothing
    RB.W3.grad = nothing
    RB.b1.grad = nothing
    RB.b2.grad = nothing
end

get_params(RB::ResidualBlock) = [RB.W1, RB.W2, RB.W3, RB.b1, RB.b2]


text:
module OHDSICohortExpressions

include("datamodel.jl")
include("expression.jl")
include("translate.jl")

end


text:
using Profile, PProf
using Pkg;Pkg.activate(".")
using PairwiseVelocities
using StaticArrays
using BenchmarkTools
using DelimitedFiles

boxsize = 2000.
rbins = LinRange(0.,50,50)
rbins_c = 0.5*(rbins[2:end] + rbins[1:end-1])

filename = "moments.csv"
logn = -5.
run = 101
snapshot = 15

positions, velocities = read_data(run, snapshot, 10^logn, boxsize)
positions = convert(Array{Float64}, positions)
velocities = convert(Array{Float64}, velocities)

moments = PairwiseVelocities.compute_pairwise_velocity_moments(
    positions,
    velocities,
    rbins,
    boxsize,
)
open(filename; write=true) do f
    write(f, "# r_c v_r sigma_r sigma_t skewness_r skewness_rt kurtosis_r kurtosis_t kurtosis_rt \n")
    writedlm(f,
            zip(rbins_c, moments[2][:],
            moments[3][:], moments[6][:],
            moments[4][:], moments[7][:],
            moments[5][:], moments[8][:], moments[9][:]))
    end

println("Wrote file !")


text:
@testset "240.search-a-2d-matrix-ii.jl" begin
arr = [
    1   4   7   11  15;
    2   5   8   12  19;
    3   6   9   16  22;
    10  13  14  17  24;
    18  21  23  26  30;
]

@test searchMatrix(arr, 5) == true
@test searchMatrix(arr, 20) == false
end


text:
#########
# @pass #
#########

"""
```
Cassette.@pass transform
```

Return a Cassette pass that can be provided to the `Context` constructor's `pass` keyword
argument in order to apply `transform` to the lowered IR representations of all methods
invoked during contextual execution.

`transform` must be a Julia object that is callable with the following signature:

```
transform(::Type{<:Context}, ::Cassette.Reflection)::Union{Expr,CodeInfo}
```

If `isa(transform(...), Expr)`, then the returned `Expr` will be emitted immediately
without any additional processing. Otherwise, if `isa(transform(...), CodeInfo)`, then
the returned `CodeInfo` will undergo the rest of Cassette's overdubbing transformation
before being emitted from the `overdub` generator.

Two special `Expr` heads are available to Cassette pass authors that are not normally valid
in Julia IR. `Expr`s with these heads can be used to interact with the downstream built-in
Cassette passes that consume them.

- `:nooverdub`: Wrap an `Expr` with this head value around the first argument in an
    `Expr(:call)` to tell downstream built-in Cassette passes not to overdub that call. For
    example, `Expr(:call, Expr(:nooverdub, GlobalRef(MyModule, :myfunc)), args...)`.

- `:contextslot`: Cassette will replace any `Expr(:contextslot)` with the actual `SlotNumber`
    corresponding to the context object associated with the execution trace. For example, one
    could construct an IR element that accesses the context's `metadata` field by emitting:
    `Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), Expr(:contextslot), QuoteNode(:metadata))`

Cassette provides a few IR-munging utility functions of interest to pass authors; for details,
see [`insert_statements!`](@ref), [`replace_match!`](@ref), and [`is_ir_element`](@ref).

Note that the `@pass` macro expands to an `eval` call and thus should only be called at
top-level. Furthermore, to avoid world-age issues, `transform` should not be overloaded after
it has been registered with `@pass`.

Note also that `transform` should be "relatively pure." More specifically, Julia's compiler
has license to apply `transform` multiple times, even if only compiling a single method
invocation once. Thus, it is required that `transform` always return a generally "equivalent"
`CodeInfo` for a given context, method body, and signature. If your `transform`
implementation is not naturally "pure" in this sense, then it is still possible to guarantee
this property by memoizing your implementation (i.e. maintaining a cache of previously
computed IR results, instead of recomputing results every time).

See also: [`Context`](@ref), [`overdub`](@ref)
"""
macro pass(transform)
    Pass = gensym("PassType")
    line = Expr(:quote, __source__.line)
    file = Expr(:quote, __source__.file)
    return esc(quote
        import Cassette.__overdub_generator__
        struct $Pass <: $Cassette.AbstractPass end
        (::Type{$Pass})(ctxtype, reflection) = $transform(ctxtype, reflection)
        Core.eval($__module__, $Cassette.overdub_definition($line, $file))
        $Pass()
    end)
end

#############
# utilities #
#############

"""
```
replace_match!(replace, ismatch, x)
```

Return `x` with all subelements `y` replaced with `replace(y)` if `ismatch(y)`. If
`!ismatch(y)`, but `y` is of type `Expr`, `Array`, or `SubArray`, then replace `y`
in `x` with `replace_match!(replace, ismatch, y)`.

Generally, `x` should be of type `Expr`, `Array`, or `SubArray`.

Note that this function modifies `x` (and potentially its subelements) in-place.

See also: [`insert_statements!`](@ref), [`is_ir_element`](@ref)
"""
function replace_match!(replace, ismatch, x)
    if ismatch(x)
        return replace(x)
    elseif isa(x, Array) || isa(x, SubArray)
        for i in eachindex(x)
            x[i] = replace_match!(replace, ismatch, x[i])
        end
    elseif isa(x, Expr)
        replace_match!(replace, ismatch, x.args)
    end
    return x
end

"""
```
insert_statements!(code::Vector, codelocs::Vector, stmtcount, newstmts)
```

For every statement `stmt` at position `i` in `code` for which `stmtcount(stmt, i)` returns
an `Int`, remove `stmt`, and in its place, insert the statements returned by
`newstmts(stmt, i)`. If `stmtcount(stmt, i)` returns `nothing`, leave `stmt` alone.

For every insertion, all downstream `SSAValue`s, label indices, etc. are incremented
appropriately according to number of inserted statements.

Proper usage of this function dictates that following properties hold true:

- `code` is expected to be a valid value for the `code` field of a `CodeInfo` object.
- `codelocs` is expected to be a valid value for the `codelocs` field of a `CodeInfo` object.
- `newstmts(stmt, i)` should return a `Vector` of valid IR statements.
- `stmtcount` and `newstmts` must obey `stmtcount(stmt, i) == length(newstmts(stmt, i))` if
    `isa(stmtcount(stmt, i), Int)`.

To gain a mental model for this function's behavior, consider the following scenario. Let's
say our `code` object contains several statements:

```
code = Any[oldstmt1, oldstmt2, oldstmt3, oldstmt4, oldstmt5, oldstmt6]
codelocs = Int[1, 2, 3, 4, 5, 6]
```

Let's also say that for our `stmtcount` returns `2` for `stmtcount(oldstmt2, 2)`, returns `3`
for `stmtcount(oldstmt5, 5)`, and returns `nothing` for all other inputs. From this setup, we
can think of `code`/`codelocs` being modified in the following manner:

```
newstmts2 = newstmts(oldstmt2, 2)
newstmts5 = newstmts(oldstmt5, 5)
code = Any[oldstmt1,
           newstmts2[1], newstmts2[2],
           oldstmt3, oldstmt4,
           newstmts5[1], newstmts5[2], newstmts5[3],
           oldstmt6]
codelocs = Int[1, 2, 2, 3, 4, 5, 5, 5, 6]
```

See also: [`replace_match!`](@ref), [`is_ir_element`](@ref)
"""
function insert_statements!(code, codelocs, stmtcount, newstmts)
    ssachangemap = fill(0, length(code))
    labelchangemap = fill(0, length(code))
    worklist = Tuple{Int,Int}[]
    for i in 1:length(code)
        stmt = code[i]
        nstmts = stmtcount(stmt, i)
        if nstmts !== nothing
            addedstmts = nstmts - 1
            push!(worklist, (i, addedstmts))
            ssachangemap[i] = addedstmts
            if i < length(code)
                labelchangemap[i + 1] = addedstmts
            end
        end
    end
    Core.Compiler.renumber_ir_elements!(code, ssachangemap, labelchangemap)
    for (i, addedstmts) in worklist
        i += ssachangemap[i] - addedstmts # correct the index for accumulated offsets
        stmts = newstmts(code[i], i)
        @assert length(stmts) == (addedstmts + 1)
        code[i] = stmts[end]
        for j in 1:(length(stmts) - 1) # insert in reverse to maintain the provided ordering
            insert!(code, i, stmts[end - j])
            insert!(codelocs, i, codelocs[i])
        end
    end
end

"""
```
is_ir_element(x, y, code::Vector)
```

Return `true` if `x === y` or if `x` is an `SSAValue` such that
`is_ir_element(code[x.id], y, code)` is `true`.

See also: [`replace_match!`](@ref), [`insert_statements!`](@ref)
"""
function is_ir_element(x, y, code::Vector)
    result = false
    while true # break by default
        if x === y #
            result = true
            break
        elseif isa(x, SSAValue)
            x = code[x.id]
        else
            break
        end
    end
    return result
end


text:
using TrajectoryOptimization
using Ipopt
using MathOptInterface
const MOI = MathOptInterface
const TO = TrajectoryOptimization

if !isdefined(Main,:TEST_TIME)
    TEST_TIME = true
end

# Parallel Park
prob = DubinsCarProblem(:parallel_park)
TO.add_dynamics_constraints!(prob)

nlp = TO.TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)

optimizer = Ipopt.Optimizer()
TO.build_MOI!(nlp, optimizer)
MOI.optimize!(optimizer)
@test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
@test cost(nlp) < 0.0541
@test max_violation(nlp) < 1e-11
# TEST_TIME && @test optimizer.solve_time < 0.1

@test norm(states(nlp)[end] - prob.xf) < 1e-10
@test norm(states(nlp)[1] - prob.x0) < 1e-10

# Cartpole
prob = CartpoleProblem()
# prob = Problems.Cartpole()[1]
TO.add_dynamics_constraints!(prob)

nlp = TO.TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
optimizer = Ipopt.Optimizer()
TO.build_MOI!(nlp, optimizer)
MOI.optimize!(optimizer)
@test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
@test cost(nlp) < 1.50
@test max_violation(nlp) < 1e-11
# TEST_TIME && @test optimizer.solve_time < 1

@test norm(states(nlp)[end] - prob.xf) < 1e-10
@test norm(states(nlp)[1] - prob.x0) < 1e-10

# # Pendulum
# prob = Problems.Pendulum()[1]
# TO.add_dynamics_constraints!(prob)
#
# nlp = TO.TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
# optimizer = Ipopt.Optimizer()
# TO.build_MOI!(nlp, optimizer)
# MOI.optimize!(optimizer)
# @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
# @test cost(nlp) < 0.0185
# @test max_violation(nlp) < 1e-9
# # TEST_TIME && @test optimizer.solve_time < 0.1
#
# @test norm(states(nlp)[end] - prob.xf) < 1e-10
# @test norm(states(nlp)[1] - prob.x0) < 1e-10
#
# # 3 Obstacles
# prob = Problems.DubinsCar(:three_obstacles)[1]
# TO.add_dynamics_constraints!(prob)
#
# nlp = TO.TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
# optimizer = Ipopt.Optimizer()
# TO.build_MOI!(nlp, optimizer)
# MOI.optimize!(optimizer)
# @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
# @test cost(nlp) < 12.1
# @test max_violation(nlp) < 1e-8
# # TEST_TIME && @test optimizer.solve_time < 0.5
#
# @test norm(states(nlp)[end] - prob.xf) < 1e-10
# @test norm(states(nlp)[1] - prob.x0) < 1e-10
#
# # Escape
# prob = Problems.DubinsCar(:escape)[1]
# TO.add_dynamics_constraints!(prob)
#
# nlp = TO.TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
# optimizer = Ipopt.Optimizer()
# TO.build_MOI!(nlp, optimizer)
# MOI.optimize!(optimizer)
# @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
# @test cost(nlp) < 0.333
# @test max_violation(nlp) < 1e-8
# # TEST_TIME && @test optimizer.solve_time < 5
# @test norm(states(nlp)[end] - prob.xf) < 1e-10
# @test norm(states(nlp)[1] - prob.x0) < 1e-10


text:
push!(LOAD_PATH,"../")

using Documenter, OREnvironment, ORInterface

makedocs(sitename="OREnvironment.jl",
format = Documenter.HTML(prettyurls = false)
)

deploydocs(
	repo = "github.com/DavidGarHeredia/OREnvironment.jl.git",
	devbranch = "dev"
)


text:
#This example has been adopted from https://github.com/IntelLabs/ParallelAccelerator.jl/blob/master/examples/black-scholes/black-scholes.jl

function blackscholes_serial(sptprice,
                             strike,
                             rate,
                             volatility,
                             time)
    logterm = log10( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(- rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
end

@inline function cndf2(in)
    out = 0.5f0 + 0.5f0 * erf(0.707106781f0 * in)
    return out
end

function runs(iterations)
    sptprice   = Float32[ 42.0 for i = 1:iterations ]
    initStrike = Float32[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float32[ 0.5 for i = 1:iterations ]
    volatility = Float32[ 0.2 for i = 1:iterations ]
    time       = Float32[ 0.5 for i = 1:iterations ]

    sptprice_gpu = AFArray(sptprice)
    initStrike_gpu = AFArray(initStrike)
    rate_gpu = AFArray(rate)
    volatility_gpu = AFArray(volatility)
    time_gpu = AFArray(time)

    put1 = blackscholes_serial.(sptprice, initStrike, rate, volatility, time)
    put2 = blackscholes_serial.(sptprice_gpu, initStrike_gpu, rate_gpu, volatility_gpu, time_gpu)
    sync(put2)
    @test sum(put1) ≈ sum(put2)
    put2
end

@inferred runs(10^3)


text:
import Base.string
import Base.print
import Base.show
import Millboard

"""
    to_dict(m::Any) :: Dict{String,Any}
    to_string_dict(m::Any; all_output::Bool = false) :: Dict{String,String}
    to_string_dict(m::Any, fields::Array{Symbol,1}; all_output::Bool = false) :: Dict{String,String}

Creates a `Dict` using the fields and the values of `m`.
"""
function to_dict(m::Any) :: Dict{String,Any}
  Dict(string(f) => getfield(m, Symbol(f)) for f in fieldnames(typeof(m)))
end
function to_string_dict(m::Any; all_output::Bool = false) :: Dict{String,String}
  to_string_dict(m, fieldnames(m), all_output = all_output)
end
function to_string_dict(m::Any, fields::Array{Symbol,1}; all_output::Bool = false) :: Dict{String,String}
  response = Dict{String,String}()
  for f in fields
    response[string(f)] = string(getfield(m, Symbol(f)))
  end

  response
end