text:
abstract type AbstractScalarFunctionSpace{dim,T,FE<:FiniteElement,MM,N1,N2,N3,fdim} <: DiscreteFunctionSpace{dim,T,FE} end

struct ScalarFunctionSpace{dim,T<:Real,FE<:FiniteElement,MM,N1,N2,N3,fdim} <: AbstractScalarFunctionSpace{dim,T,FE,MM,N1,N2,N3,fdim}
    N::Matrix{T}
    dNdx::Matrix{Vec{dim,T}}
    dNdξ::Matrix{Vec{dim,T}}
    detJ::ScalarWrapper{T}
    M::Matrix{T}
    dMdξ::Matrix{Vec{dim,T}}
    qr_weights::Vector{T}
    fe::FE
    mesh::PolygonalMesh{dim,N1,N2,N3,T}
    Jinv::ScalarWrapper{Tensor{2,dim,T,MM}}
    update_face_data::Bool
    L::Array{T,3}
    dLdξ::Array{Vec{dim,T},3}
    E :: Matrix{Vector{T}}
    detJf::Vector{T}
    normals::Vector{Vec{dim,T}}
    qr_face_weigths::Vector{T}
    qr_face_points::Vector{Vec{fdim,T}}
end

function ScalarFunctionSpace(mesh::PolygonalMesh, felem::FiniteElement{dim,shape,order,gorder};update_face_data::Bool = true,
    quad_degree::Int = order+1) where {dim, shape, order,gorder}
    quad_rule = QuadratureRule{dim,shape}(DefaultQuad(), quad_degree)
    f_quad_rule = QuadratureRule{dim-1,shape}(DefaultQuad(), quad_degree)
    fs = _scalar_fs(Float64, mesh, quad_rule, f_quad_rule, felem, update_face_data)
end

function _scalar_fs(::Type{T}, mesh::PolygonalMesh{dim,N1,N2,N3,T}, quad_rule::QuadratureRule{dim,shape},
    f_quad_rule::QuadratureRule{fdim,shape},
    felem::FiniteElement{dim,shape,order,1}, update_face_data::Bool) where {dim, fdim,T,shape<:AbstractRefShape,N1,N2,N3,order}
    n_qpoints = length(getpoints(quad_rule))
    n_cells = getncells(mesh)
    n_faces = n_faces_per_cell(mesh)
    n_faceqpoints = length(getpoints(f_quad_rule))
    q_ref_facepoints = getpoints(f_quad_rule)
    q_ref_faceweights = getweights(f_quad_rule)
    face_quad_rule = create_face_quad_rule(f_quad_rule)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getngeombasefunctions(felem)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)

    # Face interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    E    = fill(zeros(T,n_faces) * T(NaN), n_func_basefuncs, n_faceqpoints)

    # Geometry face Interpolation
    n_geom_basefuncs = getngeombasefunctions(felem)
    L =    fill(zero(T)          * T(NaN), n_geom_basefuncs, n_faceqpoints, n_faces)
    dLdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_faceqpoints, n_faces)

    if update_face_data
        for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
            for i in 1:n_geom_basefuncs
                dLdξ[i, qp, face], L[i, qp, face] = gradient(ξ -> geom_value(felem, i, ξ), ξ, :all)
            end
        end
        for qp in 1:n_faceqpoints
            for i in 1:n_func_basefuncs
                E_f = zeros(T, n_faces)
                for j in 1:n_faces
                    #Evaluate shape function on q_point map to edge j
                    E_f[j] = value(felem, i, face_quad_rule[j].points[qp])
                end
                E[i, qp] = E_f
            end
        end
    end

    for (qp, ξ) in enumerate(quad_rule.points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp]  = gradient(ξ -> value(felem, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp], M[i, qp] = gradient(ξ -> geom_value(felem, i, ξ), ξ, :all)
        end
    end

    detJ = ScalarWrapper(T(NaN))
    detJf = fill(T(NaN), n_faces)
    normals = fill(zero(Vec{dim,T}) * T(NaN), n_faces)
    Jinv = ScalarWrapper(zero(Tensor{2,dim,T,2*dim}) )

    MM = Tensors.n_components(Tensors.get_base(typeof(Jinv[])))

    ScalarFunctionSpace{dim,T,typeof(felem),MM,N1,N2,N3,fdim}(N, dNdx, dNdξ, detJ,
    M, dMdξ, getweights(quad_rule), felem, mesh, Jinv,update_face_data, L, dLdξ, E, detJf, normals,
    q_ref_faceweights, q_ref_facepoints)
end

function reinit!(fs::ScalarFunctionSpace{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_geom_basefuncs = getngeombasefunctions(fs.fe)
    n_func_basefuncs = getnbasefunctions(fs.fe)
    @assert length(x) == n_geom_basefuncs
    fecv_J = zero(Tensor{2,dim})
    for j in 1:n_geom_basefuncs
        fecv_J += x[j] ⊗ fs.dMdξ[j, 1]
    end
    detJ = det(fecv_J)
    detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
    fs.detJ[] = detJ
    fs.Jinv[] = inv(fecv_J)

    for j in 1:n_func_basefuncs, i in 1:length(fs.qr_weights)
        fs.dNdx[j, i] = fs.dNdξ[j, i] ⋅ fs.Jinv[]
    end
    if fs.update_face_data
        n_faces = n_faces_per_cell(fs.mesh)
        shape = getshape(fs.fe)
        for l in 1:n_faces
            fef_J = zero(Tensor{2,dim})
            for j in 1:n_geom_basefuncs
                fef_J += x[j] ⊗ fs.dLdξ[j, 1, l]
            end
            weight_norm = weighted_normal(fef_J, l, shape, Val{dim})
            detJ_f = norm(weight_norm)
            fs.normals[l] = weight_norm / detJ_f
            detJ_f > 0.0 || throw(ArgumentError("det(Jf) is not positive: det(Jf) = $(detJ_f)"))
            fs.detJf[l] = detJ_f
        end
    end
end

########### Data Functions
@inline getngeobasefunctions(fs::AbstractScalarFunctionSpace) = size(fs.M, 1)
@inline getnquadpoints(fs::AbstractScalarFunctionSpace) = length(fs.qr_weights)
@inline getnbasefunctions(fs::AbstractScalarFunctionSpace) = size(fs.N,1)
@inline getdetJdV(fs::ScalarFunctionSpace, q_point::Int) = fs.detJ[]*fs.qr_weights[q_point]
@inline shape_value(fs::AbstractScalarFunctionSpace, q_point::Int, base_func::Int) = fs.N[base_func, q_point]
@inline shape_gradient(fs::ScalarFunctionSpace, q_point::Int, base_func::Int) = fs.dNdx[base_func, q_point]
@inline shape_divergence(fs::ScalarFunctionSpace, q_point::Int, base_func::Int) = sum(fs.dNdx[base_func, q_point])
@inline geometric_value(fs::AbstractScalarFunctionSpace, q_point::Int, base_func::Int) = fs.M[base_func, q_point]
@inline getdim(::AbstractScalarFunctionSpace{dim}) where {dim} = dim
@inline reference_coordinate(fs::AbstractScalarFunctionSpace{dim,T},cell::Int, mesh::PolygonalMesh, x::Vec{dim,T}) where {dim,T} = fs.Jinv[]⋅(x-mesh.nodes[mesh.cells[cell].nodes[1]].x)
@inline getfiniteelement(fs::AbstractScalarFunctionSpace) = fs.fe
@inline getnlocaldofs(fs::AbstractScalarFunctionSpace) = getnbasefunctions(fs)
@inline getmesh(fs::AbstractScalarFunctionSpace) = fs.mesh

# Face Data
@inline getnfacegeobasefunctions(fs::AbstractScalarFunctionSpace) = size(fs.L, 1)
@inline getnfacequadpoints(fs::AbstractScalarFunctionSpace) = length(fs.qr_face_weigths)
@inline getfacedetJdS(fs::ScalarFunctionSpace, face::Int, q_point::Int) = fs.detJf[face]*fs.qr_face_weigths[q_point]
@inline shape_value(fs::AbstractScalarFunctionSpace, face::Int, q_point::Int, base_func::Int, orientation::Bool = true) = orientation ? fs.E[base_func, q_point][face] : fs.E[base_func, end - q_point+1][face]
@inline geometric_value(fs::AbstractScalarFunctionSpace, face::Int, q_point::Int, base_func::Int) = fs.L[base_func, q_point, face]

"""
    getnormal(fs::ScalarFunctionSpace, cell::Int, face::Int, qp::Int)
Return the normal at the quadrature point `qp` for the face `face` at
cell `cell` of the `ScalarFunctionSpace` object.
"""
@inline get_normal(fs::ScalarFunctionSpace, face::Int) = fs.normals[face]


text:
using MacroTools
using Compat
using Compat.Test

let
  x = @match :(2+3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == (2, 3)
end

let
  x = @match :(2-3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == (3, 2)
end

let
  x = @match :(2/3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
  end
  @test x == nothing
end

let
  x = @match :(2/3) begin
    (a_+b_) => (a, b)
    (a_-b_) => (b, a)
    _ => :default
  end
  @test x == :default
end

let
  ex = :(mutable struct Foo
           x::Int
           y
         end)
  @capture(ex, mutable struct T_ fields__ end)
  @test T == :Foo
  @test fields == [:(x::Int), :y]

  @capture(ex, mutable struct T_ fields__ end)
  @test T == :Foo
  @test fields == [:(x::Int), :y]
end

let
  ex = :(f(x))
  @capture(ex, f_(xs__))
  @test f == :f
  @test xs == [:x]
end

let
  ex = :(f(x, y, z))
  @capture(ex, f_(x_, xs__))
  @test f == :f
  @test x == :x
  @test xs == [:y, :z]
end

let
  ex = quote
    function foo(a, b)
      return a+b
    end
  end
  @assert @capture(shortdef(ex), f_(args__) = body_)
end

let
  ex = :(a = b)
  @capture(ex, a_ = b_)
  @test (a, b) == (:a, :b)
end

let
  ex = :(f(a = b))
  @capture(ex, f(a_ = b_))
  @test (a, b) == (:a, :b)
  @capture(ex, f(x_))
  @test isexpr(x, :kw)
end

macro nothing_macro()
end
@test @expand(@nothing_macro) === nothing

macro splitcombine(fundef) # should be a no-op
    dict = splitdef(fundef)
    esc(MacroTools.combinedef(dict))
end

# Macros for testing that splitcombine doesn't break
# macrocalls in bodies
macro zeroarg()
   :(1)
end
macro onearg(x)
   :(1+$(esc(x)))
end

let
    # Ideally we'd compare the result against :(function f(x)::Int 10 end),
    # but it fails because of :line and :block differences
    @test longdef(:(f(x)::Int = 10)).head == :function
    @test longdef(:(f(x::T) where U where T = 2)).head == :function
    @test shortdef(:(function f(x)::Int 10 end)).head != :function
    @test map(splitarg, (:(f(a=2, x::Int=nothing, y, args...))).args[2:end]) ==
        [(:a, :Any, false, 2), (:x, :Int, false, :nothing),
         (:y, :Any, false, nothing), (:args, :Any, true, nothing)]
    @test splitarg(:(::Int)) == (nothing, :Int, false, nothing)

    @splitcombine foo(x) = x+2
    @test foo(10) == 12
    @splitcombine add(a, b=2; c=3, d=4)::Float64 = a+b+c+d
    @test add(1; d=10) === 16.0
    @splitcombine fparam(a::T) where {T} = T
    @test fparam([]) == Vector{Any}
    struct Orange end
    @splitcombine (::Orange)(x) = x+2
    @test Orange()(10) == 12
    @splitcombine fwhere(a::T) where T = T
    @test fwhere(10) == Int
    @splitcombine manywhere(x::T, y::Vector{U}) where T <: U where U = (T, U)
    @test manywhere(1, Number[2.0]) == (Int, Number)
    @splitcombine fmacro0() = @zeroarg
    @test fmacro0() == 1
    @splitcombine fmacro1() = @onearg 1
    @test fmacro1() == 2

    struct Foo{A, B}
        a::A
        b::B
    end
    # Parametric outer constructor
    @splitcombine Foo{A}(a::A) where A = Foo{A, A}(a,a)
    @test Foo{Int}(2) == Foo{Int, Int}(2, 2)
end

include("destruct.jl")

ext: jl


text:
using Test, Pkg

@testset "Utilities" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false
    for submodule in ["TicToc", "VariableTemplates"]
        if all_tests || "$submodule" in ARGS || "Utilities" in ARGS
            include_test(submodule)
        end
    end

end


text:
#Input.jl
"""
Handling input for DEB-IBM simulations.
Contains functions for parameter conversions and corrections, as well as dictionaries with AMP values and default simulation settings.
"""

using DocStringExtensions

#### ----- default global dictionaries ---- ####

# parameters for population experiment with D. magna conducted in early 2019
population_dict = Dict(
    "remove_juveniles" => [], # weekdays on which to remove juveniles
    "renew_medium" => [0, 3], # weekdays on which to remove medium
    "frct_medium_renewal" => 0.25, # fraction of medium to replace per renewal
    "N0_juvenile" => 8, # initial juveniles
    "N0_adult"    => 2, # initial adults
    "timestep"    => 24, # number of timesteps per day
    "tmax"        => 42, # max. simulated time (days)
    "food_level"  => 8e7, # food level (# cells/L/day)
    "age_range_init" => 24, # maximum initial age of juveniles (timesteps)
    "stressor_conc" => [0., 0., 0.]
)

# settings of a life-table experiment as conducted in late 2019/early 2020
life_table_dict = Dict(
    "remove_juveniles" => [0, 2, 5], # remove juveniles 3x per week
    "renew_medium"     => [0, 2, 5], # renew media on the same days
    "frct_medium_renewal" => 1.00, # 100% of media are renewed
    "N0_juvenile" => 1, # simulation starts with 1 juvenila
    "N0_adult"    => 0,
    "timestep"    => 24, # number of timesteps per day
    "tmax"        => 21, # test lasts for 21 days
    "food_level"  => 5., # food level in mgC//L3/day
    "age_range_init" => 24, # maximum initial age of juveniles (timesteps)
    "stressor_conc" => [0., 0., 0.]
)

# reference scenario for ecorelevance project
ecorelevance_reference_dict = Dict(
    "remove_juveniles" => [],
    "renew_medium" => [0,3], # renew media on mondays and thursdays
    "frct_medium_renewal" => 0.25, # renew 25% media
    "N0_juvenile" => 3, # start with 3 juveniles
    "N0_adult" => 1, # and one adult
    "age_range_init" => 48,
    "timestep" => 24,
    "tmax" => 56,
    "food_level" => 2.5/1000, # food level in mgC/cm^3/day
    "stressor_1_conc" => 0, # concentration of stressor 1
    "stressor_2_conc" => 0, # concentration of stressor 2
    "stressor_3_conc" => 0  # concentration of stressor 3
)

"""
Conversion of AMP parameters for use in IBM.
See user manual of the generic DEB-IBM implementation for explanations.
$(SIGNATURES)
"""
function convert_parameters(fld)
    # maintenance rate follows from volume-specific somatic structure
    # and volume-specific costs for structure
    fld["k_M_rate"]    = fld["p_M"]  / fld["E_G"]
     # compound parameters following from scaling by assimilation rate
    fld["g"]           = fld["E_G"]  * fld["v_rate"] / (fld["kap"] * fld["p_Am"])
    fld["U_Hb"]        = fld["E_Hb"] / fld["p_Am"]
    fld["U_Hp"]        = fld["E_Hp"] / fld["p_Am"]
    # functional response parameters derived from primary parameters
    fld["J_XAm_rate"]  = fld["p_Am"] / fld["y_EX"]
    fld["K"]           = fld["J_XAm_rate"] / fld["F_m"]
    return fld
end

"""
Estimate embryonal parameters using bisection method.
"""
function estimate_egg_weight(
        timestep::Int64,
        v_rate::Float64,
        U_Hb::Float64,
        kap::Float64,
        g::Float64,
        k_J_rate::Float64,
        k_M_rate::Float64
)
    embryo_timestep = timestep * 1000
    lower_bound = 0
    upper_bound = 1
    maxiter = 100
    egg_weight = 0.5
    U_H_embryo = 0
    U_E_embryo = 0
    e_scaled_embryo = 0
    e_ref = 1
    L_embryo = 0.00001

    for i in 1:maxiter
        egg_weight = .5 * (lower_bound + upper_bound)
        L_embryo = 0.00001
        U_E_embryo = egg_weight
        U_H_embryo = 0
        e_scaled_embryo = v_rate * (U_E_embryo/L_embryo)^3

        while (U_H_embryo < U_Hb) & (e_scaled_embryo >1)
             # calculate reserve density based on energy conductance, scaled reserves and structural volume
            e_scaled_embryo = v_rate * (U_E_embryo/L_embryo^3)
            # calculate corresponding mobilization flux
            S_C_embryo      = L_embryo^2 * (g * e_scaled_embryo / ( g + e_scaled_embryo)) * (1 + (L_embryo / (g * (v_rate / (g * k_M_rate)))))
            # calculate change in state varaibles
            dU_E_embryo     = -1 * S_C_embryo
            dU_H_embryo     = (1-kap)*S_C_embryo - k_J_rate*U_H_embryo
            dL_embryo       = ((1/3) * (((v_rate / (g*(L_embryo^2))) * S_C_embryo) - k_M_rate*L_embryo))
            # update state variables
            U_E_embryo = U_E_embryo + dU_E_embryo / embryo_timestep
            U_H_embryo = U_H_embryo + dU_H_embryo / embryo_timestep
            L_embryo   = L_embryo   + dL_embryo   / embryo_timestep
            e_scaled_embryo = v_rate * (U_E_embryo / L_embryo^3)
        end

        # in order to accept an estimate:
        # 1. estimated reserve density of embryo has to be close to one
        # (to fulfil the criterion that embryos are mostly reserve)
        # 2. estimated maturity has to be larger or equal to maturity at birth (by definition)
        if (e_scaled_embryo < 0.01+e_ref) & (e_scaled_embryo > -0.01+e_ref) & (U_H_embryo >= U_Hb)
            @goto escape_label # jump to @label escape_label
        end
        # if estimate has not been accpeted:
        # adjust upper/lower bound based on result of previous estimation
        if U_H_embryo > U_Hb
            upper_bound = egg_weight
        else
            lower_bound = egg_weight
        end
        #if i == maxiter
        #    error("Embryo submodel is not converging. Smaller timesteps could help.")
        #end
    end
    @label escape_label

    return egg_weight, L_embryo, U_E_embryo, U_H_embryo
end

"""
Estimation of intitial state variables of an Embryo, determining the amount of energy
that has to be invested per offspring.
$(SIGNATURES)
"""
function estimate_egg_weight(fleas_dict::Dict)
    # estimate embryonal parameters
    fleas_dict["egg_weight"],
    fleas_dict["L_embryo"],
    fleas_dict["U_E_embryo"],
    fleas_dict["U_H_embryo"] = estimate_egg_weight(
        24,
        fleas_dict["v_rate"],
        fleas_dict["U_Hb"],
        fleas_dict["kap"],
        fleas_dict["g"],
        fleas_dict["k_J_rate"],
        fleas_dict["k_M_rate"])
    return fleas_dict
end

#### ---- default flea dictionaries ---- ####
"""
Parameter values for Daphnia magna according to AMP database.
Last checked: Jan 2020
Toxicity parameters from Daphnia magna models (Pereira et al. 2019; Karel Vlaeminck, unpublished).
"""
DM_AMP = Dict(
    "species" => "DM_AMP",
    "cv"      => 0.1,

    # AMP parameters
    "p_Am" => 313.169, #
    "F_m"  => 30., # filtration rate (cm^3/day*cm^2)
    "E_Hb"       => 0.05464,
    "E_Hp"       => 1.09,
    "p_M"        => 1200,
    "E_G"        => 4400,
    "k_J_rate"   => 0.2537,
    "v_rate"     => 0.1858,
    "kap"        => 0.5809,
    "kap_R"      => 0.95,
    "shape_factor"       => 0.264,

    "sG"  => -0.3,
    "h_a" => 0.0002794,
    "time_between_molts" => 2.324,

    # embryonal parameters
    "L0"  => 0.0001,
    "egg_weight" => 3.771781921386719E-4,
    "L_embryo"   => 0.023340451224020977,
    "U_E_embryo" => 6.039673323323045E-5,

    "starvation_hazard_rate" => 0.,
    "crit_mass" => 0.,

    # Toxicity parameters
    # Toxicity parameters
    "pmoas"  => ["Growth", "Maintenance", "Assimilation", "Reproduction"],
    "k_d"    => [ 0., 0., 0.],
    "h_max"  => [ 0., 0., 0.],
    "h_beta" => [ 0., 0., 0.],
    "h_ED50" => [ 0., 0., 0.],
    "S_max"  => [ 0., 0., 0.],
    "S_beta" => [ 0., 0., 0.],
    "S_ED50" => [ 0., 0., 0.]
)

"""
Parameter values for Daphnia magna, fed with mixture of Pseudokirchneriella subcapitata and Chlamydomonas reinhardtii, simulated as a single food source.
Food densities in # cells/L.
"""
DM1_2019 = Dict(
    "species" => "DM1_2019",
    "cv"      => 0.1,

    # corrected parameters

    "p_Am" => 260.4, # (CI = (260.3, 309.5))
    "F_m"  => 10.0, # (CI = (8, 11))
    "y_EX" => 3.7e-6, # (CI = (3e-6, 4e-6))
    "mortality_constant" => 0.41, # (CI = (0.3, 0.63))

    # AMP parameters

    "E_Hb"       => 0.05464,
    "E_Hp"       => 1.09,
    "p_M"        => 1200,
    "E_G"        => 4400,
    "k_J_rate"   => 0.2537,
    "v_rate"     => 0.1858,
    "kap"        => 0.5809,
    "kap_R"      => 0.95,
    "shape_factor"       => 0.264,

    "sG"  => -0.3,
    "h_a" => 0.0002794,
    "time_between_molts" => 2.324,

    # embryonal parameters
    "L0"  => 0.0001,
    "egg_weight" => 3.771781921386719E-4,
    "L_embryo"   => 0.023340451224020977,
    "U_E_embryo" => 6.039673323323045E-5,

    "starvation_hazard_rate" => 0.,
    "crit_mass" => 0.,

    "toxicity_particle" => Dict(
        "stressors" => ["Cu", "Ni", "Zn"],
        "pMoAs" => ["Growth", "Maintenance", "Assimilation", "Reproduction"],
        # one k_d for Cu, Ni and Zn
        "k_d" => [0., 0., 0.],
        # one parameter for every stressor (rows) and every pMoA (columns)
        "S_max" => hcat(
            [0., 0., 0., 0.],
            [0., 0., 0., 0.],
            [0., 0., 0., 0.]
        ),
        "S_ED50" => fill(Inf, 4, 3),
        "S_beta" => fill(2., 4, 3),
        "h_max" => [0., 0., 0.],
        "h_ED50" => [Inf, Inf, Inf],
        "h_beta" => [2., 2., 2.]
    )
)

DL_AMP = Dict(
    "species" => "DL_AMP",
    "cv"      => 0.1,

    # AMP parameters
    "E_Hb"       => 0.01746,
    "E_Hp"       => 1.401,
    "p_M"        => 5686.83,
    "E_G"        => 4400.59,
    "k_J_rate"   => 0.002,
    "v_rate"     => 0.020455,
    "kap"        => 0.5568,
    "kap_R"      => 0.95,
    "shape_factor"       => 0.18946,
    "mortality_constant" => 0,
    "crit_mass"          => 0.4,
    "starvation_hazard_rate" => 0.9,
    "sG"  => 0.001,
    "h_a" => 8.603e-5,
    "L0"  => 0.0001,

    # food-dependent parameters
    "p_Am" => 558.735, # J/d*cm^2
    "y_EX" => 1e-3, # yield of reserve (J) on food (# cells/cm^3)
    "F_m"  => 270.833, # filtration rate (cm^3/day*cm^2)

    # embryonal parameters
    "egg_weight" => 3.771781921386719E-4,
    "L_embryo"   => 0.023340451224020977,
    "U_E_embryo" => 6.039673323323045E-5,

    # reproductive timing
    "time_between_molts" => 2.324,

    # Toxicity parameters
    # compound 1
    "compound_1" => "Copper",
    "pmoa_1"     => "Growth",
    "h_b_1"      => 0/24, # background hazard rate
    "k_d_1"      => 0, # dominant rate constant
    "k_k_1"      => 0, # ?
    "m_n_1"      => 0, # ?
    "s_s_1"      => 0,
    "s_d_1"      => 0, # ?
    "A_eq_1" => 0,
    "B_eq_1" => 0,
    "C_eq_1" => 0,

    # compound 2
    "compound_2" => "Nickel",
    "pmoa_2"     => "Maintenance",
    "h_b_2"      => 0,
    "k_d_2"      => 0,
    "k_k_2"      => 0,
    "m_n_2"      => 0,
    "s_s_2"      => 0,
    "s_d_2"      => 0,
    "A_eq_2"     => 0,
    "B_eq_2"     => 0,
    "C_eq_2"     => 0,

    # compound 3
    "compound_3" => "Zinc",
    "pmoa_3"     => "Assimilation",
    "h_b_3"      => 0,
    "k_d_3"      => 0,
    "k_k_3"      => 0,
    "m_n_3"      => 0,
    "s_s_3"      => 0,
    "s_d_3"      => 0,
    "A_eq_3"     => 0,
    "B_eq_3"     => 0,
    "C_eq_3"     => 0
)

DL1 = Dict(
    "species" => "DL_AMP",
    "cv"      => 0.1,

    # AMP parameters
    "E_Hb"       => 0.01746,
    "E_Hp"       => 1.401,
    "p_M"        => 5686.83,
    "E_G"        => 4400.59,
    "k_J_rate"   => 0.002,
    "v_rate"     => 0.020455,
    "kap"        => 0.331941,
    "kap_R"      => 0.95,
    "shape_factor"       => 0.18946,
    "mortality_constant" => 0.93567,
    "sG"  => 0.001,
    "h_a" => 8.603e-5,
    "L0"  => 0.0001,

    # food-dependent parameters
    "p_Am" => 651.572, # 2020-03-24 J/d*cm^2
    "y_EX" => 71312.7, # yield of reserve (J) on food (# cells/cm^3)
    "F_m"  => 606.58, # filtration rate (cm^3/day*cm^2)

    # embryonal parameters
    "egg_weight" => 3.771781921386719E-4,
    "L_embryo"   => 0.023340451224020977,
    "U_E_embryo" => 6.039673323323045E-5,

    # other
    "time_between_molts" => 2.324,
    "crit_mass"          => 0.442002,
    "starvation_hazard_rate" => 0.784624,

    # Toxicity parameters
    # compound 1
    "compound_1" => "Copper",
    "pmoa_1"     => "Growth",
    "h_b_1"      => 0/24, # background hazard rate
    "k_d_1"      => 0, # dominant rate constant
    "k_k_1"      => 0, # ?
    "m_n_1"      => 0, # ?
    "s_s_1"      => 0,
    "s_d_1"      => 0, # ?
    "A_eq_1" => 0,
    "B_eq_1" => 0,
    "C_eq_1" => 0,

    # compound 2
    "compound_2" => "Nickel",
    "pmoa_2"     => "Maintenance",
    "h_b_2"      => 0,
    "k_d_2"      => 0,
    "k_k_2"      => 0,
    "m_n_2"      => 0,
    "s_s_2"      => 0,
    "s_d_2"      => 0,
    "A_eq_2"     => 0,
    "B_eq_2"     => 0,
    "C_eq_2"     => 0,

    # compound 3
    "compound_3" => "Zinc",
    "pmoa_3"     => "Assimilation",
    "h_b_3"      => 0,
    "k_d_3"      => 0,
    "k_k_3"      => 0,
    "m_n_3"      => 0,
    "s_s_3"      => 0,
    "s_d_3"      => 0,
    "A_eq_3"     => 0,
    "B_eq_3"     => 0,
    "C_eq_3"     => 0
)

DL1_ECOREL = Dict(
    "species" => "DL_AMP",
    "cv"      => 0.1,

    # AMP parameters
    "E_Hb"       => 0.01746,
    "E_Hp"       => 1.401,
    "p_M"        => 5686.83,
    "E_G"        => 4400.59,
    "k_J_rate"   => 0.002,
    "v_rate"     => 0.020455,
    "kap"        => 0.331941,
    "kap_R"      => 0.95,
    "shape_factor"       => 0.18946,
    "sG"  => 0.001,
    "h_a" => 8.603e-5,
    "L0"  => 0.0001,

    # food-dependent parameters
    "p_Am" => 651.572, # 2020-03-24 J/d*cm^2
    "y_EX" => 71312.7, # yield of reserve (J) on food (# cells/cm^3)
    "F_m"  => 606.58, # filtration rate (cm^3/day*cm^2)

    # embryonal parameters
    "egg_weight" => 3.771781921386719E-4,
    "L_embryo"   => 0.023340451224020977,
    "U_E_embryo" => 6.039673323323045E-5,

    # reproductive timing
    "time_between_molts" => 2.324,

    # starvation
    "mortality_constant" => 0.93567,
    "starvation_hazard_rate" => 0.784624,
    "crit_mass" => 0.442002,

    # Toxicity parameters
    # compound 1
    "compound_1" => "Copper",
    "pmoa_1"     => "Growth",
    "h_b_1"      => 0/24, # background hazard rate
    "k_d_1"      => 0, # dominant rate constant
    "k_k_1"      => 0, # ?
    "m_n_1"      => 0, # ?
    "s_s_1"      => 0,
    "s_d_1"      => 0, # ?
    "A_eq_1" => 0,
    "B_eq_1" => 0,
    "C_eq_1" => 0,

    # compound 2
    "compound_2" => "Nickel",
    "pmoa_2"     => "Maintenance",
    "h_b_2"      => 0,
    "k_d_2"      => 0,
    "k_k_2"      => 0,
    "m_n_2"      => 0,
    "s_s_2"      => 0,
    "s_d_2"      => 0,
    "A_eq_2"     => 0,
    "B_eq_2"     => 0,
    "C_eq_2"     => 0,

    # compound 3
    # Zinc toxicity parameters based on mean µg dissolved Zn/L
    "compound_3" => "Zinc",
    "pmoa_3"     => "Maintenance",
    "k_d_3"      => 0.727032,
    "m_n_3"      => 344.759,
    "k_k_3"      => 2.67,
    "s_s_3"      => 0.01,
    "A_eq_3"     => 298.85,
    "B_eq_3"     => 2.11,
    "C_eq_3"     => 9.40,

    # TODO: remove these parameters entirely from all code
    "h_b_3" => 0,
    "s_d_3" => 0
)


text:
typealias NormMetric Union{Euclidean, Chebyshev, Cityblock, Minkowski, WeightedEuclidean, WeightedCityblock, WeightedMinkowski, Mahalanobis}

immutable HyperSphere{N, T <: AbstractFloat}
    center::SVector{N, T}
    r::T
end

HyperSphere{N, T1, T2}(center::SVector{N, T1}, r::T2) = HyperSphere(center, convert(T1, r))

@inline function intersects{T <: AbstractFloat, N, M <: Metric}(m::M,
                                                             s1::HyperSphere{N, T},
                                                             s2::HyperSphere{N, T})
    evaluate(m, s1.center, s2.center) <= s1.r + s2.r
end

@inline function encloses{T <: AbstractFloat, N, M <: Metric}(m::M,
                                                           s1::HyperSphere{N, T},
                                                           s2::HyperSphere{N, T})
    evaluate(m, s1.center, s2.center) + s1.r <= s2.r
end

@inline function interpolate{V <: AbstractVector, M <: NormMetric}(::M,
                                                                  c1::V,
                                                                  c2::V,
                                                                  x,
                                                                  d,
                                                                  ab)
    alpha = x / d
    @assert length(c1) == length(c2)
    @inbounds for i in eachindex(ab.center)
        ab.center[i] = (1 - alpha) .* c1[i] + alpha .* c2[i]
    end
    return ab.center, true
end

@inline function interpolate{V <: AbstractVector, M <: Metric}(::M,
                                                               c1::V,
                                                               ::V,
                                                               ::Any,
                                                               ::Any,
                                                               ::Any)
    return c1, false
end

function create_bsphere{V}(data::Vector{V}, metric::Metric, indices::Vector{Int}, low, high, ab)
    n_dim = size(data,1)
    n_points = high - low + 1
    # First find center of all points
    fill!(ab.center, 0.0)
    for i in low:high
        for j in 1:length(ab.center)
            ab.center[j] += data[indices[i]][j]
        end
    end
    scale!(ab.center, 1 / n_points)

    # Then find r
    r = zero(get_T(eltype(V)))
    for i in low:high
        r = max(r, evaluate(metric, data[indices[i]], ab.center))
    end
    r += eps(get_T(eltype(V)))
    return HyperSphere(SVector{length(V), eltype(V)}(ab.center), r)
end

# Creates a bounding sphere from two other spheres
function create_bsphere{N, T <: AbstractFloat}(m::Metric,
                                               s1::HyperSphere{N, T},
                                               s2::HyperSphere{N, T},
                                               ab)
    if encloses(m, s1, s2)
        return HyperSphere(s2.center, s2.r)
    elseif encloses(m, s2, s1)
        return HyperSphere(s1.center, s1.r)
    end

    # Compute the distance x along a geodesic from s1.center to s2.center
    # where the new center should be placed (note that 0 <= x <= d because
    # neither s1 nor s2 contains the other)
    dist = evaluate(m, s1.center, s2.center)
    x = 0.5 * (s2.r - s1.r + dist)
    center, is_exact_center = interpolate(m, s1.center, s2.center, x, dist, ab)
    if is_exact_center
        rad = 0.5 * (s2.r + s1.r + dist)
    else
        rad = max(s1.r + evaluate(m, s1.center, center), s2.r + evaluate(m, s2.center, center))
    end

    return HyperSphere(SVector{N, T}(center), rad)
end


text:
function col2type(A,splitSep)

    r,c,v = find(A);
    cType, cVal = SplitStr(c,splitSep);
    return Assoc(r,Array{Union{AbstractString,Number}}(cType),Array{Union{AbstractString,Number}}(cVal));

end


text:
ComponentSol{T,N,C} = DiffEqBase.AbstractODESolution{T,N,C} where C<:AbstractVector{<:ComponentArray}

# Plotting stuff
function DiffEqBase.getsyms(sol::ComponentSol)
    if DiffEqBase.has_syms(sol.prob.f)
        return sol.prob.f.syms
    else
        return labels(sol.u[1])
    end
end

# A little bit of type piracy. Should probably make this a PR to DiffEqBase
DiffEqBase.cleansyms(syms::AbstractArray{<:String}) = DiffEqBase.cleansyms.(syms)
DiffEqBase.cleansyms(syms::String) = syms

# DiffEqBase.interpret_vars(var::AbstractString, sol::ComponentSol, syms) =
#     DiffEqBase.interpret_vars(_to_sol_index(var, syms), sol, syms)
DiffEqBase.interpret_vars(var::Union{AbstractString, Symbol}, sol::ComponentSol, syms) =
    DiffEqBase.interpret_vars(_to_sol_index(var, syms), sol, syms)
# This doesn't work with mixed strings and integers
DiffEqBase.interpret_vars(vars::Tuple{Vararg{<:Union{<:AbstractString, Symbol}}}, sol::ComponentSol, syms) =
    DiffEqBase.interpret_vars(map(var->_to_sol_index(var, syms), vars), sol, syms)
DiffEqBase.interpret_vars(vars::AbstractArray{<:Union{<:AbstractString, Symbol}}, sol::ComponentSol, syms) =
    DiffEqBase.interpret_vars(collect(Iterators.flatten(map(var->_to_sol_index(var, syms), vars))), sol, syms)

_to_sol_index(var::Symbol, syms) = _to_sol_index(string(var), syms)
_to_sol_index(var::AbstractString, syms) = findall(startswith.(syms, var))
_to_sol_index(var, syms) = var

function (p::DiffEqBase.DefaultLinSolve)(x,A::ComponentMatrix,b,update_matrix=false;tol=nothing, kwargs...)
    if p.iterable isa Vector && eltype(p.iterable) <: LinearAlgebra.BlasInt # `iterable` here is the pivoting vector
        F = LU{eltype(A)}(A, p.iterable, zero(LinearAlgebra.BlasInt))
        ldiv!(x, F, b)
        return nothing
    end

    if update_matrix
        blasvendor = BLAS.vendor()
        if ArrayInterface.can_setindex(x) && (size(A,1) <= 100 || ((blasvendor === :openblas || blasvendor === :openblas64) && size(A,1) <= 500))
            p.A = RecursiveFactorization.lu!(A)
        else
            p.A = lu!(A)
        end
    end

    x .= b
    ldiv!(p.A,x)
    return nothing
end


text:
function nchoosek(n, k; with_replacement=true)
    perms = []
    all_branches = collect(1:n)
    for i in 2:k
        next_branches = collect(1:n)
        if !with_replacement
            deleteat!(next_branches, i)
        end
        all_branches = add_branches(all_branches, next_branches)
    end
    return all_branches
end

function add_branches(roots, options)
    new_branches = []
    for r in roots
        for o in options
            push!(new_branches, [r o])
        end
    end
    return new_branches
end

function calc_permutations(n, k)
    perms = []
    current_branches = collect(1:n)
    if k==1; return current_branches; end

    all_branches = nothing
    for i in 2:k
        tmp = copy(current_branches)
        all_branches = nothing
        for b in tmp
            next_branches = collect(b[end]:n)
            b = add_branches([b], next_branches)
            if all_branches == nothing
                all_branches = b
            else
                all_branches = vcat(all_branches, b)
            end
        end
        current_branches = all_branches
    end
    return all_branches
end

export nchoosek, calc_permutations


text:
# This script looks at the coefficient plot shown in the first example of [Robust causal interference
# using direct acyclic graphs: the R package `dagitty`](https://academic.oup.com/ije/article/45/6/1887/2907796)
# by J Textor, B van der Zanden, M Gilthorpe, M Lisliewicz and G Ellison.
# [R Dagitty](https://doi.org/10.1093/ije/dyw341)

using StatisticalRethinking, StatsPlots

ProjDir = @__DIR__

for i in 1:5
  include(scm_path("..", "examples", "sports", "m$i.jl"))
end

models=[m1s, m2s, m3s, m4s, m5s];
pars=[:bC, :bFL, :bPGP, :bTM]

df = DataFrame!(CSV.File(scm_path("..", "examples", "sports", "sports.csv"), delim=','));

scatter(df[:, :TM], df[:, :WUE])
savefig("$(ProjDir)/sports_01.png")

r1 = plotcoef(models, pars, "$(ProjDir)/sports_02.png", "Normal (Particles) estimates")
display(r1)

println()
r2 = plotcoef(models, pars, "$(ProjDir)/sports_03.png", "Quap estimates", quap)
display(r2)


text:
#########################################################################
using AutoPreallocation, BenchmarkTools

foo() = ones(1, 2096) * ones(2096, 1024) * ones(1024,1)
@btime foo()

const foo_res, foo_record = record_allocations(foo)

@btime avoid_allocations($foo_record, foo)

#########################################################################

using Test
twos(dims) = 2*ones(dims)

const twos_res, twos_record = record_allocations(twos, (3,6))
@test_throws TypeError avoid_allocations(twos_record, twos, (3,6,9))

# If the type is the same and only size differs AutoPreallocation right now won't even
# error, it will just silently return the wrong result
avoid_allocations(twos_record, twos, (30,60))

###############################################################################

function bar(x)
    ys = fill(10.0, 30)
    for ii in 1:100
        push!(ys, x)
    end
    return ys
end

@btime bar(3.14);

const _, bar_record = record_allocations(bar, 3.14);
reinitialize!(bar_record);
@btime avoid_allocations($bar_record, bar, 42.0);

reinitialize!(bar_record);
@btime avoid_allocations($bar_record, bar, 24601);

############################################

function mat(x)
    out = zeros(2,2)
    out[1, 1] = x
    out[2, 2] = x
    return out
end

const mat1, mat_record = record_allocations(mat, 1);
(mat1,)

mat2 = avoid_allocations(mat_record, mat, 2);
(mat1, mat2)  # Notice mat1 has changed

mat3 = avoid_allocations(mat_record, mat, 3);
(mat1, mat2, mat3)  # Notice: mat1 and mat2 have changed


text:
using Revise
using PyPlot
using Infiltrator
using LinearAlgebra
using Bem2d

"""
    bemsolve(g, rho, lambda, mu)

Solve BEM with particular integral approach with gravitysquareparticular
"""
function bemsolve(els, nels, g, rho, lambda, mu, nu, x, y)
    # Build BEM operator, TH
    idx = getidxdict(els)
    T_pU_qall, _ = PUTC(slip2dispstress, els, idx["B"], 1:1:els.endidx, mu, nu)
    _, H_pT_qall = PUTC(slip2dispstress, els, [idx["R"] ; idx["T"]; idx["L"]], 1:1:els.endidx, mu, nu)
    TH = [T_pU_qall ; H_pT_qall]

    # Modified displacment BCs for bottom
    Upart, _ = gravityparticularfunctions(els.xcenter[idx["B"]], els.ycenter[idx["B"]], g, rho, lambda, mu)
    bcs = zeros(2 * els.endidx)
    bcs[1:2:2*nels] = -Upart[:, 1]
    bcs[2:2:2*nels] = -Upart[:, 2]

    # Modified traction BCs for top
    bcs[4*nels+2:2:6*nels] = -rho * g * (els.ycenter[idx["T"]] .* els.ynormal[idx["T"]])
    Ueff = TH \ bcs

    # Forward solution on grid
    Ucomp, _ = constdispstress(slip2dispstress, x, y, els, [idx["B"] ; idx["R"] ; idx["T"] ; idx["L"]],
                               Ueff[1:2:end], Ueff[2:2:end], mu, nu)
    Uint, _ = gravityparticularfunctions(x, y, g, rho, lambda, mu)
    Utotal = Ucomp .+ Uint
    Umag = sqrt.(Utotal[:, 1].^2 + Utotal[:, 2].^2)
    return Umag
end

"""
    gravitysquareparticular()

Experiments with gravity body force.
"""
function gravitysquareparticular()
    close("all")
    fontsize = 20
    mu = 3e10
    lambda = 3e10
    nu = 0.25
    rho = 2700
    g = 9.81
    nels = 100
    npts = 200
    L = 1e4
    offset = 100

    ##### Below ground box
    top = 0
    x, y = obsgrid(-L+offset, top-2*L+offset, L-offset, top-offset, npts)
    els = Elements(Int(1e5))
    addelsez!(els, discretizedline(-L, top-2*L, L, top-2*L, nels)..., "B") # Bottom
    addelsez!(els, discretizedline(L, top-2*L, L, top, nels)... , "R") # Right hand side
    addelsez!(els, discretizedline(L, top, -L, top, nels)..., "T") # Top
    addelsez!(els, discretizedline(-L, top, -L, top-2*L, nels)..., "L") # Left hand side

    Umag = bemsolve(els, nels, g, rho, lambda, mu, nu, x, y)
    figure(figsize=(15, 5))
    subplot(1, 3, 1)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts), reshape(Umag, npts, npts), 50)
    colorbar()
    gca().set_aspect("equal")
    title("y(top) = 0")

    ##### Now with top of box not at zero
    top = 2*L
    x, y = obsgrid(-L+offset, top-2*L+offset, L-offset, top-offset, npts)
    els = Elements(Int(1e5))
    addelsez!(els, discretizedline(-L, top-2*L, L, top-2*L, nels)..., "B") # Bottom
    addelsez!(els, discretizedline(L, top-2*L, L, top, nels)... , "R") # Right hand side
    addelsez!(els, discretizedline(L, top, -L, top, nels)..., "T") # Top
    addelsez!(els, discretizedline(-L, top, -L, top-2*L, nels)..., "L") # Left hand side

    Umagoffset = bemsolve(els, nels, g, rho, lambda, mu, nu, x, y)
    subplot(1, 3, 2)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts), reshape(Umagoffset, npts, npts), 50)
    colorbar()
    gca().set_aspect("equal")
    title("y(top) = yoffset")

    subplot(1, 3, 3)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts), reshape(abs.(Umagoffset - Umag), npts, npts), 50)
    colorbar()
    gca().set_aspect("equal")
    title("residual")

end
gravitysquareparticular()


text:
"""
Compute the cost function for a model and loss function given a vector of
parameters, a dataset and a vector of frequency counts as weights.
"""
function cost_function{T <: Model}(
    ::Type{T},
    loss::Loss,
    λ::Real,
    parameters::Vector,
    inputs::Vector,
    weights::Vector,
)
    # Construct a model object that allows stack allocation of model parameters.
    m = T(parameters)

    # We iterate over all observations in the dataset.
    n = length(inputs)

    # Iterate over the full dataset accumating the total cost.
    total_cost = 0.0
    for i in 1:n
        @inbounds x2, t2, x1, t1, y = inputs[i]
        @inbounds w = weights[i]
        p = predict(m, x2, t2, x1, t1)
        total_cost += w * loss_value(loss, p, y)
    end

    # Add L2 regularization.
    total_cost += λ * norm(parameters)^2

    # Return the total cost.
    return total_cost
end

"""
Generate a closure that computes the cost function for a model given only a
vector of parameters as an argument.
"""
function cost_function_closure{T <: Model}(
    ::Type{T},
    loss::Loss,
    λ::Real,
    inputs::Vector,
    weights::Vector,
)
    return parameters -> cost_function(T, loss, λ, parameters, inputs, weights)
end

"""
Compute the gradient of a cost function for a model and loss function given a
vector of parameters, a dataset and a vector of frequency counts as weights.
"""
function cost_function_gradient!{T <: Model}(
    gr::Vector,
    ::Type{T},
    loss::Loss,
    λ::Real,
    parameters::Vector,
    inputs::Vector,
    weights::Vector,
)
    # Set the gradient to zero because we will accumulate the gradient by
    # mutating summation.
    fill!(gr, 0.0)

    # Construct a model object that allows stack allocation of model parameters.
    m = T(parameters)

    # We iterate over all observations in the dataset.
    n = length(inputs)
    K = length(parameters)

    # Iterate over the full dataset accumating the gradient.
    for i in 1:n
        @inbounds x2, t2, x1, t1, y = inputs[i]
        @inbounds w = weights[i]
        p = predict(m, x2, t2, x1, t1)
        gradient_tuple = gradient_component(m, x2, t2, x1, t1)
        σ = w * loss_derivative(loss, p, y)
        for j in 1:K
            @inbounds gr[j] += σ * gradient_tuple[j]
        end
    end

    # Add L2 regularization.
    for j in 1:K
        @inbounds gr[j] += λ * 2.0 * parameters[j]
    end

    # Return nothing. The gradient is already stored in the mutated argument.
    return
end

"""
Generate a closure that computes the gradient of the cost function for a model
given only a vector of parameters and a mutable vector that will be updated to
contain the gradient as arguments.
"""
function cost_function_gradient_closure{T <: Model}(
    ::Type{T},
    loss::Loss,
    λ::Real,
    inputs::Vector,
    weights::Vector,
)
    return (parameters, gr) -> cost_function_gradient!(
        gr,
        T,
        loss,
        λ,
        parameters,
        inputs,
        weights,
    )
end

"""
Generate the closures we use for model fitting.
"""
function make_closures{T <: Model}(
    ::Type{T},
    loss::Loss,
    λ::Real,
    inputs::Vector,
    weights::Vector,
)
    return (
        cost_function_closure(T, loss, λ, inputs, weights),
        cost_function_gradient_closure(T, loss, λ, inputs, weights),
    )
end

"""
Fit a model given inputs and weights.

Set loss to the type of loss function that will be minimized when fitting the
model's parameters.

Set λ to the regularization parameter to increase the L2 penalty on the model
parameters in the unconstrained space.

Set Θ₀ to any vector of initial parameters for the model being fit. Defaults to
initial_parameters(T) for a model of type T.

Set method to any method provided by the Optim package. Defaults to
Optim.Newton(), which is a straight-forward implementation of Newton's method.
"""
function fit{T <: Model}(
    ::Type{T},
    inputs::Vector,
    weights::Vector,
    loss::Loss = LogLoss(),
    λ::Real = 0.0,
    Θ₀::Vector = initial_parameters(T),
    method::Optim.Optimizer = Optim.Newton(),
)
    cost, cost_gradient! = make_closures(T, loss, λ, inputs, weights)

    opts = Optim.OptimizationOptions(
        f_tol = 0.0,
        show_trace = false,
        autodiff = true,
    )

    # TODO: Provide analytic Hessians to Optim.optimize.
    res = Optim.optimize(cost, cost_gradient!, Θ₀, method, opts)

    # Require convergence.
    if !Optim.converged(res)
        error("Failed to converge")
    end

    return res
end


text:
include("../required.jl")

## Get layers
# Define coordinates range
coords = (left=-145.0, right=-50.0, bottom=20.0, top=75.0)

# Test loading variables
lc_vars = SimpleSDMPredictor(Copernicus, LandCover, 1:10; resolution=10.0, coords...)

# Load worldclim variables to compare
wc_vars = SimpleSDMPredictor(WorldClim, BioClim, 1:19; resolution=10.0, coords...);

## Temporary dimensions fix
# Check which line is extra
wc_test = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords]
wc_vars[1].grid
wc_test.grid # last line seems extra
isequal(wc_test.grid[1:(end - 1), :], wc_vars[1].grid) # true

# Test clip fix
wc_clip = wc_test[top=coords.top - stride(wc_test, 2)]
isequal(wc_clip.grid, wc_vars[1].grid) # true
wc_vars[1].top
wc_clip.top # not exactly the same
isapprox(wc_vars[1].top, wc_clip.top) # approximately the same

# Test clip from world layer instead
wc_test2 = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords][top=coords.top -
                                                                 stride(wc_test, 2)] # same limits
wc_test2.top # not helping

# Fix landcover layers
lc_vars = [l[top=coords.top - stride(l, 2)] for l in lc_vars]
isequal(latitudes(wc_vars[1]), latitudes(lc_vars[1])) # false
isapprox(latitudes(wc_vars[1]), latitudes(lc_vars[1])) # true so should be ok

## Export to CSV
# Combine environmental data
env_vars = [wc_vars; lc_vars]
# Create env matrix
env_mat = mapreduce(x -> vec(x.grid), hcat, env_vars)
replace!(x -> isnothing(x) ? NaN : x, env_mat)
# Create env dataframe
env_df = DataFrame(env_mat, :auto)
rename!(env_df, vcat(Symbol.("wc", 1:size(wc_vars, 1)), Symbol.("lc", 1:size(lc_vars, 1))))
insertcols!(env_df, 1, :site => 1:nrow(env_df))

# Get sites latitudes
lats = repeat(collect(latitudes(wc_vars[1])); outer=size(wc_vars[1].grid, 2))
# Get sites longitudes
lons = repeat(collect(longitudes(wc_vars[1])); inner=size(wc_vars[1].grid, 1))
# Create spa matrix
spa_mat = [lats lons]
# Create spa dataframe
spa_df = DataFrame(; site=eachindex(lats), lat=lats, lon=lons)

## Check difference with previously saved data

# Copy data
run(
    `cp ./data/proc/distributions_env_full.csv ./data/proc/distributions_env_full_backup.csv`,
)
run(
    `cp ./data/proc/distributions_spa_full.csv ./data/proc/distributions_spa_full_backup.csv`,
)
run(`cp ./data/proc/distributions_spa_qc.csv ./data/proc/distributions_spa_qc_backup.csv`)

# Export new files
CSV.write(joinpath("data", "proc", "distributions_env_full.csv"), env_df; delim="\t")
CSV.write(joinpath("data", "proc", "distributions_spa_full.csv"), spa_df; delim="\t")
# Git diff says spa files are different
# Sizes are different too
# Looking at the diff, it seems to be a float approximation difference, which should be fine

# Let's make sure it's just that
spa_df_backup = CSV.read(
    joinpath("data", "proc", "distributions_spa_full_backup.csv"), DataFrame
)
spa_df_backup
spa_df
isequal(spa_df_backup, spa_df) # false
isapprox(spa_df_backup, spa_df) # true
isequal(spa_df_backup.site, spa_df.site) # true
isapprox(spa_df_backup.lat, spa_df.lat) # true
isapprox(spa_df_backup.lon, spa_df.lon) # true
# Yep, should just be approximation error, so its fine

# Make sure env files are equal too, even if diff says they are
env_df_backup = CSV.read(
    joinpath("data", "proc", "distributions_env_full_backup.csv"), DataFrame
)
env_df_backup
env_df
isequal(env_df_backup, env_df) # false
isapprox(env_df_backup, env_df) # false ?!
isapprox(NaN, NaN) # oh probably just because of NaN
# whatever let's trust Git

## Export QC sites coordinates (for smaller scale analyses)
coords_qc = (left=-80.0, right=-55.0, bottom=45.0, top=63.0)
# Get site indices
spa_qc = filter(
    x ->
        (coords_qc.left - stride(wc_vars[1]; dims=2) <= x.lon < coords_qc.right) &&
            (coords_qc.bottom - stride(wc_vars[1]; dims=1) <= x.lat < coords_qc.top),
    spa_df,
)

# Test dimensions for Quebec coordinates
wc_vars_qc = SimpleSDMPredictor(WorldClim, BioClim, 1; coords_qc...)
wc_test_qc = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords_qc]
isequal(size(wc_vars_qc), size(wc_test_qc)) # false
isapprox(longitudes(wc_vars_qc), unique(spa_qc.lon)) # true
isapprox(longitudes(wc_test_qc), unique(spa_qc.lon)) # false
isequal(length(wc_vars_qc.grid), nrow(spa_qc)) # equal
# So dimensions in the exported file should be ok
CSV.write(joinpath("data", "proc", "distributions_spa_qc.csv"), spa_qc; delim="\t")
# Git diff says there's difference
spa_qc_backup = CSV.read(
    joinpath("data", "proc", "distributions_spa_qc_backup.csv"), DataFrame
)
isapprox(spa_qc_backup, spa_qc) # true
# Again just an approximation error

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"""
    open-loop policy
"""
mutable struct OpenLoop{T} <: Policy{T}
    u::Vector{Vector{T}} # nominal controls
    idx::Int
    cnt::Int
    N_sample::Int
end

function open_loop_policy(u; N_sample = 1)
    OpenLoop(u, 0, N_sample, N_sample)
end

function policy(p::OpenLoop, traj, t)
    # reset
    if t == 1
        p.idx = 0
        p.cnt = p.N_sample
    end

    if p.cnt == p.N_sample
        p.idx += 1
        p.cnt = 0
    end

    p.cnt += 1

    return p.u[p.idx] ./ p.N_sample
end


text:
using Test
using StaticCompiler
using Libdl
using LinearAlgebra
using LoopVectorization
using ManualMemory
using Distributed
using StaticTools
using StrideArraysCore
addprocs(1)
@everywhere using StaticCompiler, StrideArraysCore

const GROUP = get(ENV, "GROUP", "All")

@static if GROUP == "Core" || GROUP == "All"
    include("testcore.jl")
end

@static if GROUP == "Integration" || GROUP == "All"
    include("testintegration.jl")
end


text:
using Test
using ReinforcementLearning
using ReinforcementLearning.Utils
using ReinforcementLearningEnvironments

using StatsBase
using Random
using Flux

@testset "ReinforcementLearning" begin
    include("Utils/utils.jl")
    include("components/components.jl")
    include("glue/glue.jl")
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
@gsp " " :-  # reset session
@gsp "unset key" "set multi layout 1,2 title 'Multiplot title'" :-
@gsp :- 1 x y Z "w l palette" "set view 55, 65" :-
@gsp :- 2 x y Z "w pm3d" "set view 55, 65" "set key off" :-
@gsp :- x y Z "w l lc 'white'"
save(term="pngcairo size 800,400", output="plt3d_ex2_2.png")


text:
using Base64
struct Viewer
    fn::String
end
"""
visualize gif result at Plots pane
See https://discourse.julialang.org/t/base-show-and-plotpane/17835/3
"""
function Base.show(io::IO, b::Viewer)
    mime = "image/png"
    suffix = ";base64,"
    str = stringmime(mime, read(b.fn))
    str = string("<img src=\"data:", mime, suffix, str, "\">")
    Juno.render(Juno.PlotPane(), HTML(str))
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
##############################################################################
#
# All Document types share a common metadata profile using DocumentMetadata
#
##############################################################################

mutable struct DocumentMetadata
    language
    title::String
    author::String
    timestamp::String
end
DocumentMetadata() = DocumentMetadata(
    Languages.English(),
    "Untitled Document",
    "Unknown Author",
    "Unknown Time"
)

##############################################################################
#
# The abstract Document type
#
##############################################################################

abstract type AbstractDocument; end

##############################################################################
#
# FileDocument type and constructors
#
##############################################################################

mutable struct FileDocument <: AbstractDocument
    filename::String
    metadata::DocumentMetadata
end

function FileDocument(f::AbstractString)
    d = FileDocument(String(f), DocumentMetadata())
    d.metadata.title = f
    return d
end

##############################################################################
#
# StringDocument type and constructors
#
##############################################################################

mutable struct StringDocument{T<:AbstractString} <: AbstractDocument
    text::T
    metadata::DocumentMetadata
end

StringDocument(txt::AbstractString) = StringDocument(txt, DocumentMetadata())

##############################################################################
#
# TokenDocument type and constructors
#
##############################################################################

mutable struct TokenDocument{T<:AbstractString} <: AbstractDocument
    tokens::Vector{T}
    metadata::DocumentMetadata
end
function TokenDocument(txt::AbstractString, dm::DocumentMetadata)
    TokenDocument(tokenize(dm.language, String(txt)), dm)
end
function TokenDocument(tkns::Vector{T}) where T <: AbstractString
    TokenDocument(tkns, DocumentMetadata())
end
TokenDocument(txt::AbstractString) = TokenDocument(String(txt), DocumentMetadata())

##############################################################################
#
# NGramDocument type and constructors
#
##############################################################################

mutable struct NGramDocument{T<:AbstractString} <: AbstractDocument
    ngrams::Dict{T,Int}
    n::Union{Int,Vector{Int}}
    metadata::DocumentMetadata
end
function NGramDocument(txt::AbstractString, dm::DocumentMetadata, n::Integer...=1)
    NGramDocument(ngramize(dm.language, tokenize(dm.language, String(txt)), n...), (length(n) == 1) ? Int(first(n)) : Int[n...], dm)
end
function NGramDocument(txt::AbstractString, n::Integer...=1)
    NGramDocument(txt, DocumentMetadata(), n...)
end
function NGramDocument(ng::Dict{T, Int}, n::Integer...=1) where T <: AbstractString
    NGramDocument(merge(Dict{AbstractString,Int}(), ng), (length(n) == 1) ? Int(first(n)) : Int[n...], DocumentMetadata())
end

##############################################################################
#
# text() / text!(): Access to document text as a string
#
##############################################################################

function text(fd::FileDocument)
    !isfile(fd.filename) && error("Can't find file: $(fd.filename)")
    read(fd.filename, String)
end

text(sd::StringDocument) = sd.text
function text(td::TokenDocument)
    @warn("TokenDocument's can only approximate the original text")
    join(td.tokens, " ")
end
function text(ngd::NGramDocument)
    error("The text of an NGramDocument cannot be reconstructed")
end

text!(sd::StringDocument, new_text::AbstractString) = (sd.text = new_text)
function text!(d::AbstractDocument, new_text::AbstractString)
    error("The text of a $(typeof(d)) cannot be edited")
end

##############################################################################
#
# tokens() / tokens!(): Access to document text as a token array
#
##############################################################################

tokens(d::(Union{FileDocument, StringDocument})) = tokenize(language(d), text(d))
tokens(d::TokenDocument) = d.tokens
function tokens(d::NGramDocument)
    error("The tokens of an NGramDocument cannot be reconstructed")
end

tokens!(d::TokenDocument, new_tokens::Vector{T}) where {T <: AbstractString} = (d.tokens = new_tokens)
function tokens!(d::AbstractDocument, new_tokens::Vector{T}) where T <: AbstractString
    error("The tokens of a $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# ngrams() / ngrams!(): Access to document text as n-gram counts
#
##############################################################################

function ngrams(d::NGramDocument, n::Integer)
    error("The n-gram complexity of an NGramDocument cannot be increased")
end
ngrams(d::AbstractDocument, n::Integer...) = ngramize(language(d), tokens(d), n...)
ngrams(d::NGramDocument) = d.ngrams
ngrams(d::AbstractDocument) = ngrams(d, 1)

ngrams!(d::NGramDocument, new_ngrams::Dict{AbstractString, Int}) = (d.ngrams = new_ngrams)
function ngrams!(d::AbstractDocument, new_ngrams::Dict)
    error("The n-grams of $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

function Base.length(d::NGramDocument)
    error("NGramDocument's do not have a well-defined length")
end

Base.length(d::AbstractDocument) = length(text(d))

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

ngram_complexity(ngd::NGramDocument) = ngd.n

function ngram_complexity(d::AbstractDocument)
    error("$(typeof(d))'s have no n-gram complexity")
end

##############################################################################
#
# Union type that refers to a generic, non-abstract document type
#
##############################################################################

const GenericDocument = Union{
    FileDocument,
    StringDocument,
    TokenDocument,
    NGramDocument
}

##############################################################################
#
# Easier Document() constructor that decides types based on inputs
#
##############################################################################

Document(str::AbstractString) = isfile(str) ? FileDocument(str) : StringDocument(str)
Document(tkns::Vector{T}) where {T <: AbstractString} = TokenDocument(tkns)
Document(ng::Dict{String, Int}) = NGramDocument(ng)

##############################################################################
#
# Conversion rules
#
##############################################################################

function Base.convert(::Type{StringDocument}, d::FileDocument)
    StringDocument(text(d), d.metadata)
end
function Base.convert(::Type{TokenDocument}, d::(Union{FileDocument, StringDocument}))
    TokenDocument(tokens(d), d.metadata)
end
function Base.convert(::Type{NGramDocument},
            d::(Union{FileDocument, StringDocument, TokenDocument}))
    NGramDocument(ngrams(d), 1, d.metadata)
end
Base.convert(::Type{TokenDocument}, d::TokenDocument) = d
Base.convert(::Type{NGramDocument}, d::NGramDocument) = d

##############################################################################
#
# getindex() methods: StringDocument("This is text and that is not")["is"]
#
##############################################################################

Base.getindex(d::AbstractDocument, term::AbstractString) = ngrams(d)[term]


text:
using Distributions

# FIXME implement actual distribution using Distributions.MixtureModel
"""
A mixture of 2 Cauchy distributions.
Random values are further constrained to `[0.0, 1.0]` range either by
truncating the initial unconstrained value or by generating new random value
until it fits the range.
"""
struct BimodalCauchy
    a::Cauchy
    b::Cauchy
    mix_prob::Float64
    # When sampling it is common to truncate in either or both ends.
    clampBelow0::Bool
    clampAbove1::Bool

    BimodalCauchy(location1, scale1, location2, scale2;
                  mix_prob::Number=0.5, clampBelow0::Bool=true, clampAbove1::Bool=true) =
        new(Cauchy(location1, scale1), Cauchy(location2, scale2),
            mix_prob, clampBelow0, clampAbove1)
end

function Random.rand(distr::BimodalCauchy)
    while true
        value = rand() < distr.mix_prob ? rand(distr.a) : rand(distr.b)
        if value >= 1.0
            distr.clampAbove1 && return 1.0
        elseif value <= 0.0
            distr.clampBelow0 && return 0.0
        else
            return value
        end
    end
end

function Random.rand!(distr::BimodalCauchy, A::AbstractArray)
    for i in eachindex(A)
        @inbounds A[i] = rand(distr)
    end
    return A
end


text:
using OrdinaryDiffEq, Test, DiffEqDevTools, DiffEqOperators, Random, LinearAlgebra
using OrdinaryDiffEq: alg_order

@testset "Caching Out-of-place" begin
  println("Caching Out-of-place")
  μ = 1.01
  linnonlin_f2 = (u,p,t) -> μ * u
  linnonlin_f1 = DiffEqScalar(μ)
  linnonlin_fun = SplitFunction(linnonlin_f1, linnonlin_f2; analytic=(u0,p,t)->u0.*exp.(2μ*t))
  prob = SplitODEProblem(linnonlin_fun,1/2,(0.0,1.0))

  Random.seed!(100)
  dts = 1 ./2 .^(7:-1:4) #14->7 good plot
  for Alg in [LawsonEuler,NorsettEuler,ETDRK2,ETDRK3,ETDRK4,HochOst4,ETD2,KenCarp3]
    sim  = test_convergence(dts,prob,Alg())
    @test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.2
  end
  # Dense test
  sim  = test_convergence(dts,prob,ETDRK4(),dense_errors=true)
  @test abs(sim.𝒪est[:l2]-4) < 0.2
end

@testset "Caching Inplace" begin
  println("Caching Inplace")
  μ = 1.01
  u0 = rand(2)
  A = [2.0 -1.0; -1.0 2.0]
  linnonlin_f1 = DiffEqArrayOperator(A)
  linnonlin_f2 = (du,u,p,t) -> du .= μ .* u
  linnonlin_fun_iip = SplitFunction(linnonlin_f1,linnonlin_f2;analytic=(u0,p,t)->exp((A+μ*I)*t)*u0)
  prob = SplitODEProblem(linnonlin_fun_iip,u0,(0.0,1.0))

  dts = 1 ./2 .^(8:-1:4) #14->7 good plot
  for Alg in [LawsonEuler(),NorsettEuler(),ETDRK2(),ETDRK3(),ETDRK4(),HochOst4(),ETD2(),KenCarp3(linsolve=LinSolveGMRES(tol=1e-6))]
    sim  = test_convergence(dts,prob,Alg)
    @test sim.𝒪est[:l2] ≈ alg_order(Alg) atol=0.1
  end
  sim  = test_convergence(dts,prob,ETDRK4(),dense_errors=true)
  @test sim.𝒪est[:l2] ≈  4 atol=0.1
  @test sim.𝒪est[:L2] ≈ 4 atol=0.1
end

@testset "EPIRK Out-of-place" begin
  println("EPIRK Out-of-place")
  # Setup nonlinear problem
  A = [-2.0 1.0; 1.0 -2.0]
  f = (u,p,t) -> A*u - u.^3
  jac = (u,p,t) -> A - [3u[1]^2 0.0; 0.0 3u[2]^2]
  fun = ODEFunction(f; jac=jac)
  Random.seed!(0); u0 = rand(2); tspan = (0.0, 1.0)
  prob = ODEProblem(fun, u0, tspan)
  # Setup approximate solution
  test_setup = Dict(:alg=>Vern9(), :reltol=>1e-16, :abstol=>1e-16)
  # Convergence simulation
  dts = 1 ./2 .^(7:-1:4)
  Algs = [Exp4, EPIRK4s3A, EPIRK4s3B, EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2]
  for Alg in Algs
    sim = analyticless_test_convergence(dts, prob, Alg(adaptive_krylov=false), test_setup)
    if Alg == EPIRK5s3
      @test_broken sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
    else
      @test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
    end
  end
end

@testset "EPIRK Inplace" begin
  println("EPIRK Inplace")
  # Setup nonlinear problem
  A = [-2.0 1.0; 1.0 -2.0]
  f = (du,u,p,t) -> (mul!(du, A, u); du .-= u.^3)
  jac_update = (J,u,p,t) -> (copyto!(J,A); J[1,1] -= 3u[1]^2; J[2,2] -= 3u[2]^2)
  jac_prototype = DiffEqArrayOperator(zeros(2,2); update_func=jac_update)
  fun = ODEFunction(f; jac_prototype=jac_prototype)
  Random.seed!(0); u0 = rand(2); tspan = (0.0, 1.0)
  prob = ODEProblem(fun, u0, tspan)
  # Setup approximate solution
  test_setup = Dict(:alg=>Vern9(), :reltol=>1e-16, :abstol=>1e-16)
  # Convergence simulation
  dts = 1 ./2 .^(7:-1:4)
  Algs = [Exp4, EPIRK4s3A, EPIRK4s3B, EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2]
  for Alg in Algs
    sim = analyticless_test_convergence(dts, prob, Alg(adaptive_krylov=false), test_setup)
    if Alg == EPIRK5s3
      @test_broken sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
    else
      @test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
    end
  end
end

@testset "ExpRK Autodiff" begin
  println("ExpRK Autodiff")
  # Setup nonlinear problem (without explicit jacobian)
  A = [-2.0 1.0; 1.0 -2.0]
  f = (u,p,t) -> A*u - u.^3
  f_ip = (du,u,p,t) -> (mul!(du, A, u); du .-= u.^3)
  fun = ODEFunction(f)
  fun_ip = ODEFunction(f_ip)
  Random.seed!(0); u0 = rand(2); tspan = (0.0, 1.0)
  prob = ODEProblem(fun, u0, tspan)
  prob_ip = ODEProblem(fun_ip, u0, tspan)
  # Setup approximate solution
  test_setup = Dict(:alg=>Vern9(), :reltol=>1e-16, :abstol=>1e-16)
  # Convergence simulation
  dts = 1 ./2 .^(7:-1:4)
  sim = analyticless_test_convergence(dts, prob, HochOst4(krylov=true), test_setup)
  @test sim.𝒪est[:l2] ≈ 4 atol=0.1
  sim = analyticless_test_convergence(dts, prob_ip, HochOst4(krylov=true), test_setup)
  @test sim.𝒪est[:l2] ≈ 4 atol=0.1
  sim = analyticless_test_convergence(dts, prob, EPIRK5P1(adaptive_krylov=false), test_setup)
  @test sim.𝒪est[:l2] ≈ 5 atol=0.1
  sim = analyticless_test_convergence(dts, prob_ip, EPIRK5P1(adaptive_krylov=false), test_setup)
  @test sim.𝒪est[:l2] ≈ 5 atol=0.1
end

@testset "Adaptive Exprb Out-of-place" begin
  println("Adaptive Exprb Out-of-place")
  # Setup nonlinear problem
  A = [-2.0 1.0; 1.0 -2.0]
  f = (u,p,t) -> A*u - u.^3
  jac = (u,p,t) -> A - [3u[1]^2 0.0; 0.0 3u[2]^2]
  fun = ODEFunction(f; jac=jac)
  Random.seed!(0); u0 = rand(2); tspan = (0.0, 1.0)
  prob = ODEProblem(fun, u0, tspan)
  # Setup approximate solution
  test_setup = Dict(:alg=>Vern9(), :reltol=>1e-16, :abstol=>1e-16)
  # Convergence simulation
  dts = 1 ./2 .^(7:-1:4)
  Algs = [Exprb32, Exprb43]
  for Alg in Algs
    sim = analyticless_test_convergence(dts, prob, Alg(), test_setup)
    @test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
  end
end

@testset "Adaptive Exprb Inplace" begin
  println("Adaptive Exprb Inplace")
  # Setup nonlinear problem
  A = [-2.0 1.0; 1.0 -2.0]
  f = (du,u,p,t) -> (mul!(du, A, u); du .-= u.^3)
  jac_update = (J,u,p,t) -> (copyto!(J,A); J[1,1] -= 3u[1]^2; J[2,2] -= 3u[2]^2)
  jac_prototype = DiffEqArrayOperator(zeros(2,2); update_func=jac_update)
  fun = ODEFunction(f; jac_prototype=jac_prototype)
  Random.seed!(0); u0 = rand(2); tspan = (0.0, 1.0)
  prob = ODEProblem(fun, u0, tspan)
  # Setup approximate solution
  test_setup = Dict(:alg=>Vern9(), :reltol=>1e-16, :abstol=>1e-16)
  # Convergence simulation
  dts = 1 ./2 .^(7:-1:4)
  Algs = [Exprb32, Exprb43]
  for Alg in Algs
    sim = analyticless_test_convergence(dts, prob, Alg(), test_setup)
    @test sim.𝒪est[:l2] ≈ alg_order(Alg()) atol=0.1
  end
end


text:
# The purpose of this script is to compute the empirical mean squared error of
# various algorithms for line spetral estimation for separating three closeby
# sinusoids with frequencies [f_1, f_1 + delta, f_1 + 2*delta] for varying
# delta.
#
# Parameters:
# N: number of measurements per snapshot
# L: number of snapshots
# Niters: number of Monte Carlo iterations
# m: number of sparse measurements which can be used in some of the algorithms;
#    set m = N for full measurements
# SNR_db: signal-to-noise-ratio [dB] of the AWGN model
# f_true: location of the frequencies of the sinusoids

include("MUSIC.jl")
include("rootMUSIC.jl")
include("ESPRIT.jl")
include("TLSesprit.jl")
include("MPM.jl")
include("ANM.jl")
include("EMAC.jl")
include("periodogram.jl")
include("findMinWrapDist.jl")
include("findMinWrapDistUnique.jl")
include("WrapDist.jl")

using Plots, Random

N = 10
m = N
Niters = 3
L = 1
SNR_db = 30

SNR = 10^(SNR_db/10)

K = 3
f1 = 0.3
f_range = range(0.0,0.1,length = 5)
MSE = zeros(length(f_range),3)

A(f,N) = exp.(2*pi*1im*kron(f,(0:N-1)))

for ii = 1:length(f_range)
    f_true = [f1 f1+f_range[ii] f1 + 2*f_range[ii]]
    for jj = 1:Niters

        noise = (randn(N,L)+1im*randn(N,L))/ (sqrt(2)*10^(SNR_db/10))
        omega = sort(randperm(m)[1:N])
        S = exp.(2*pi*rand(K,L)*1im)

        Y1 = A(f_true,N)*S + noise
        Y2 = A(f_true,m)[omega,:]*S + noise
        #t0 = time_ns()
        #f_MUSIC = rootMUSIC(Y,K,SNR)
        #t1 = (time_ns()-t0)/1e9
        #MinWrapDist = findMinWrapDistUnique(f_true,f_MUSIC)
        #MSE[ii,1] += norm(MinWrapDist)^2

        t0 = time_ns()
        f_ESPRIT = rootMUSIC(Y1,K)
        t2 = (time_ns()-t0)/1e9
        MinWrapDist = findMinWrapDistUnique(f_true,f_ESPRIT)
        MSE[ii,1] += norm(MinWrapDist)^2

        #t0 = time_ns()
        #f_MPM = MPM(Y,K)
        #t3 = (time_ns()-t0)/1e9
        #MinWrapDist = findMinWrapDistUnique(f_true,f_MPM)
        #MSE[ii,3] += norm(MinWrapDist)^2

        f_EMAC = EMAC(Y2,m,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_EMAC)
        MSE[ii,2] += norm(MinWrapDist)^2

        f_ANM = ANM(Y2,m,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_ANM)
        MSE[ii,3] += norm(MinWrapDist)^2

    end
end

MSE /= Niters

plot(f_range, MSE[:,1],yaxis=:log,tickfontsize = 6,
label="root MUSIC",
#xlims = (0,1),
#xticks = round.(f_range.-f1; digits=2),
xlabel = "delta f",
ylabel = "MSE"
)
plot!(f_range,MSE[:,2],label="EMAC")
plot!(f_range,MSE[:,3],label="ANM")
#savefig("/home/thomas/atom-amd64/usr/share/atom/example/mse_delta_f.png")


text:
mutable struct OutlinedTriangle <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    color::Rendering.Palette

    polygon::Geometry.Polygon

    aabb::Geometry.AABB{Float64}

    detection::Nodes.Detection
    inside::Bool

    function OutlinedTriangle(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.color = Rendering.White()

        o.polygon = Geometry.Polygon{Float64}()

        o.aabb = Geometry.AABB{Float64}()
        o.detection = Nodes.Detection(Rendering.Lime(), Rendering.Red())
        o.inside = false

        Geometry.add_vertex!(o.polygon, -0.5, 0.5)
        Geometry.add_vertex!(o.polygon, 0.5, 0.5)
        Geometry.add_vertex!(o.polygon, 0.0, -0.5)

        Geometry.build!(o.polygon)

        Nodes.set_dirty!(o, true)

        o
    end
end

# --------------------------------------------------------
# Timing: either called by NodeManager or by a parent node
# --------------------------------------------------------
function Nodes.update(node::OutlinedTriangle, dt::Float64)
    Nodes.update!(node.detection, node)
end

function Nodes.draw(node::OutlinedTriangle, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        Rendering.transform!(context, node.polygon)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)
    Rendering.render_outlined_polygon(context, node.polygon, Rendering.CLOSED);

    aabb_color = Nodes.highlight_color(node.detection, node.inside)
    Rendering.set_draw_color(context, aabb_color)

    Geometry.set!(node.aabb, Nodes.get_bucket(node))
    Rendering.render_aabb_rectangle(context, node.aabb)
end

# --------------------------------------------------------
# Events
# --------------------------------------------------------
function Nodes.io_event(node::OutlinedTriangle, event::Events.MouseMotionEvent)
    # println("io_event ", event, ", node: ", node)
    Nodes.set_device_point!(node.detection, Float64(event.x), Float64(event.y))
    node.inside = Nodes.check!(node.detection, node, node.base.world)
end

function set!(node::OutlinedTriangle, v1::Point{Float64}, v2::Point{Float64}, v3::Point{Float64})
    Geometry.set!(node.polygon.mesh.vertices[1], v1)
    Geometry.set!(node.polygon.mesh.vertices[2], v2)
    Geometry.set!(node.polygon.mesh.vertices[3], v3)

    Nodes.set_dirty!(node, true)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function get_nonlinearcoeff(res, mode, source)
    getω(source)/c * res.material.nonlinearindex * mode.corefraction(source) / mode.effectivearea(source)
end

function get_ωmesh(tmesh)
    dt_mesh = tmesh[2]-tmesh[1]
    num_pts = length(tmesh)
    ω = collect(2π*(-num_pts/2:num_pts/2-1)/(num_pts*dt_mesh) .|> THz2π)
    ω
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 429806e2-ff34-11ea-184e-c772958b8fbd
using Pkg, DrWatson

# ╔═╡ 42984312-ff34-11ea-0ad3-2f757db72417
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 214f9dc2-ff34-11ea-07d3-793ec3053138
md"## Clip-06-18-19s.jl"

# ╔═╡ 4298b220-ff34-11ea-2b17-e7f249b9e5d4
md"##### Not yet done."

# ╔═╡ 429f6e6c-ff34-11ea-21a2-b9be4f78fbae
md"##### Will for now be attempted using StructuralCausalModels.jl."

# ╔═╡ 429fe194-ff34-11ea-0a4e-b3ce9a78ada5
md"## End of clip-06-18-19s.jl"

# ╔═╡ Cell order:
# ╟─214f9dc2-ff34-11ea-07d3-793ec3053138
# ╠═429806e2-ff34-11ea-184e-c772958b8fbd
# ╠═42984312-ff34-11ea-0ad3-2f757db72417
# ╠═4298b220-ff34-11ea-2b17-e7f249b9e5d4
# ╠═429f6e6c-ff34-11ea-21a2-b9be4f78fbae
# ╟─429fe194-ff34-11ea-0a4e-b3ce9a78ada5


text:
using DiffEqDevTools
using Base.Test

# write your own tests here
@time @testset "Benchmark Tests" begin include("benchmark_tests.jl") end
@time @testset "ODE AppxTrue Tests" begin include("ode_appxtrue_tests.jl") end
!is_windows() && @time @testset "ODE Tableau Convergence Tests" begin include("ode_tableau_convergence_tests.jl") end ## Windows 32-bit fails on Butcher62 convergence test

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
#Eye diagram generation
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
#=TODO: Restructure DataEye to be a vector of (x,y) vector pairs?
        - Less likely to be in an invalid state.
        - But less compatible with most plotting tool API (which tend to accept
          vectors of x & y vectors).
        - Possibly have an accessor function to "export" contents of DataEye
          into an (x, y) tuple storing vectors of Float64[] for direct
          compatibility with said plotting tools??
=#
mutable struct DataEye
	vx::Vector{Vector{Float64}}
	vy::Vector{Vector{Float64}}
end
DataEye() = DataEye([], [])

#==
===============================================================================#
#=TODO: - kwargs tbit, teye?:
          Preferable to use named arguments, but can it be done cleanly without
          being confusing to user?
Advanced features:
        - Define algorithm to detect first crossing and center eye
          (select `tstart`) accordingly.
=#

"""
    buildeye(x::Vector, y::Vector, tbit::Number, teye::Number; tstart::Number=0)

Builds an eye diagram by folding `x` values of provided `(x,y)` into multiple windows of `teye` that start (are "triggered") every `tbit`.

Trace data is stored in a `DataEye` object as vectors of `(x,y)` sub-vectors.

Inputs:
  - tbit: Bit period (s).  Defines "trigger point" of each bit start.
  - teye: Eye period (s).  Defines how much of eye data (along x-axis) algorithm is to collect after "trigger point".  Typically, this is set to `2.0*tbit`.
  - tstart: Time of first "trigger point" (s).

Example plotting with Plots.jl:

    #Assumption: (x, y) data generated here.
    tbit = 1e-9 #Assume data bit period is 1ns.

    #Build eye & use tstart to center data.
    eye = buildeye(x, y, tbit, 2.0*tbit, tstart=0.2*tbit)

    plot(eye.vx, eye.vy)
"""
function buildeye(x::Vector, y::Vector, tbit::Number, teye::Number; tstart::Number=0)
	eye = DataEye()

	i = 1
	#skip initial data:
	while i <= length(x) && x[i] < tstart
		i+=1
	end

	wndnum = 0
	inext = i
	while true
		if inext > length(x); break; end #Nothing else to add
		wndstart = tstart+wndnum*tbit
		nexteye = wndstart+tbit
		wndend = wndstart+teye
		istart = inext
		i = istart
		while i <= length(x) && x[i] < nexteye
			i+=1
		end
		inext = i
		while i <= length(x) && x[i] <= wndend
			i+=1
		end
		if i > length(x)
			i = length(x)
		end
		if x[i] > wndend
			i -= 1
		end
		if i > istart
			push!(eye.vx, x[istart:i].-wndstart)
			push!(eye.vy, y[istart:i])
		end
		wndnum += 1
	end
	return eye
end

#Last line


text:
using MicrobiomeAnalysis
# using SummarizedExperiments
using Test

@testset "MicrobiomeAnalysis.jl" begin

    se = SummarizedExperiments.exampleobject(40, 20)

    coldata(se).time = 1:20

    p = abundance_plot(se, "foo")

    @test @isdefined p

end


text:
function ChainRulesCore.rrule(::typeof(getproperty), x::ComponentArray, s::Union{Symbol, Val})
    function getproperty_adjoint(Δ)
        zero_x = zero(x)
        setproperty!(zero_x, s, Δ)
        return (ChainRulesCore.NoTangent(), zero_x, ChainRulesCore.NoTangent())
    end

    return getproperty(x, s), getproperty_adjoint
end

ChainRulesCore.rrule(::typeof(getdata), x::ComponentArray) = getdata(x), Δ->(ChainRulesCore.NoTangent(), ComponentArray(Δ, getaxes(x)))

ChainRulesCore.rrule(::Type{ComponentArray}, data, axes) = ComponentArray(data, axes), Δ->(ChainRulesCore.NoTangent(), getdata(Δ), ChainRulesCore.NoTangent())

max_forks_repo_name: jonniedie/ComponentArrays


text:
ppi_init() = register(DataDep(
    "PPI",
    """
    The protein-protein interaction networks from the `"Predicting
    Multicellular Function through Multi-layer Tissue Networks"
    <https://arxiv.org/abs/1707.04638>`_ paper, containing positional gene
    sets, motif gene sets and immunological signatures as features (50 in
    total) and gene ontology sets as labels (121 in total).
    """,
    "https://data.dgl.ai/dataset/ppi.zip",
    "1f5b2b09ac0f897fa6aa1338c64ab75a5473674cbba89380120bede8cddb2a6a";
    post_fetch_method=preprocess_ppi,
))

function preprocess_ppi(local_path)
    unzip_zip(local_path)

    for phase in ["train", "test", "valid"]
        graph_file = @datadep_str "PPI/$(phase)_graph.json"
        id_file = @datadep_str "PPI/$(phase)_graph_id.npy"
        X_file = @datadep_str "PPI/$(phase)_feats.npy"
        y_file = @datadep_str "PPI/$(phase)_labels.npy"

        py"""
        import numpy as np
        ids = np.load($id_file)
        X = np.load($X_file)
        y = np.load($y_file)
        """

        X = Matrix{Float32}(py"X")
        y = SparseMatrixCSC{Int32,Int64}(Array(py"y"))
        ids = Array(py"ids")
        graph = read_ppi_graph(graph_file)

        jld2file = replace(local_path, "ppi.zip"=>"ppi.$(phase).jld2")
        @save jld2file graph X y ids
    end
end

function read_ppi_graph(filename::String)
    d = JSON.Parser.parsefile(filename)
    g = SimpleDiGraph{Int32}(length(d["nodes"]))

    for pair in d["links"]
        add_edge!(g, pair["source"], pair["target"])
    end
    g
end

struct PPI <: Dataset
end

function traindata(::PPI)
    file = datadep"PPI/ppi.train.jld2"
    @load file graph X y ids
    graph, X, y, ids
end

function validdata(::PPI)
    file = datadep"PPI/ppi.valid.jld2"
    @load file graph X y ids
    graph, X, y, ids
end

function testdata(::PPI)
    file = datadep"PPI/ppi.test.jld2"
    @load file graph X y ids
    graph, X, y, ids
end


text:
# Latex tables
# Add: write text table +++

"""
    $(SIGNATURES)

Defines Latex cell colors. Stores name and intensity of a color.

# Example

```
CellColor("blue", 75)
```
"""
mutable struct CellColor
    name :: AbstractString
    intensity :: Integer
end

"""
	$(SIGNATURES)

Returns a string of the form "\\cellcolor{blue!75}".
"""
function color_string(c :: CellColor)
    iStr = c.intensity;
    if iStr <= 0
        return ""
    else
        return "\\cellcolor{" * c.name * "!" * "$iStr" * "}"
    end
end

"""
    $(SIGNATURES)

Single or multicolumn cell. Stores text, width, alignment, color.
"""
mutable struct Cell{T1 <: AbstractString, T2 <: Integer}
    text :: T1
    width :: T2
    align :: Char
    color :: CellColor
end

"""
    $(SIGNATURES)

Single `Cell` with default color and alignment.
"""
function Cell(txt :: T1) where T1 <: AbstractString
    return Cell(txt, 1, 'l', CellColor("blue", 0))
end

function Cell(txt :: T1, align :: Char) where T1 <: AbstractString
    return Cell(txt, 1, align, CellColor("blue", 0))
end

"""
	$(SIGNATURES)

Specify a cell that spans multiple columns.

# Example
```
Cell("Heading", 3)
```
"""
function Cell(txt :: T1, width :: T2) where
    {T1 <: AbstractString, T2 <: Integer}

    return Cell(txt, width, 'c', CellColor("blue", 0))
end

function Cell(txt :: T1, w :: T2, align :: Char) where
    {T1 <: AbstractString, T2 <: Integer}

    return Cell(txt, w, align, CellColor("blue", 0))
end

"""
	$(SIGNATURES)

Returns a Latex string that makes a cell. Either just the (colored) contents or a `multicolumn` command.
"""
function cell_string(c :: Cell)
    contentStr = color_string(c.color) * c.text;
    if c.width > 1
        return "\\multicolumn{$(c.width)}{$(c.align)}{$contentStr}"
    else
        return contentStr
    end
end

"""
    Table

Holds a Latex table, constructed as a vector of rows.

Flow
1. Constructor
2. Make and add rows (some are just provided by user as strings)
3. Write table

# Example
```
tb = Table(5, "lSSSS");
add_row!(tb, " & \\multicolumn{3}{c}{Heading 1} & ");
write_table(tb, "test.tex");
```
"""
mutable struct Table{T1 <: AbstractString, T2 <: Integer}
    # nRows :: T2
    nCols :: T2
    alignV :: T1
    # headerV :: Vector{T1}
    # Body is a vector of strings (rows)
    bodyV  :: Vector{T1}
end

function Table(alignV :: T1) where T1 <: AbstractString
    return Table(length(alignV),  alignV)
end

function Table(nc :: T1, alignV :: T2) where
    {T1 <: Integer, T2 <: AbstractString}

    return Table(nc, alignV, Vector{String}(undef, 0))
end

"""
    $(SIGNATURES)

Number of table rows.
"""
function nrows(tb :: Table)
    return length(tb.bodyV)
end

"""
	$(SIGNATURES)

Add a row to a table.

# Example
```
tb = Table(3, "lll");
add_row!(tb, "Column1 & Column2 & Column3");
add_row!(tb, "\\cline{2-3}")
```
"""
function add_row!(tb :: Table, rowStr :: T1) where T1 <: AbstractString
    push!(tb.bodyV, rowStr)
    return nothing
end

"""
	$(SIGNATURES)

Make a row from a vector of cells.

# Example
```
tb = Table(5, "lllll");
# One single cell + one multicolumn cell that spans 4 columns.
cellV = [LatexLH.Cell("cell1"),  LatexLH.Cell("cell2", 4, 'c')]
rowStr = LatexLH.make_row(tb, cellV);
```
"""
function make_row(tb :: Table, cellV :: Vector{T1}) where T1
    nCols = 0;
    rowStr = "";
    for i1 in 1 : length(cellV)
        nCols += cellV[i1].width;
        rowStr = rowStr * cell_string(cellV[i1]);
        if i1 < length(cellV)
            rowStr = rowStr * " & ";
        end
    end
    @assert nCols == tb.nCols  "Number of columns does not match table"
    return rowStr
end

function header(tb :: Table)
    return ["\\begin{tabular}{$(tb.alignV)}",
        "\\toprule"]
end

function footer(tb :: Table)
    return ["\\bottomrule",  "\\end{tabular}"]
end

"""
    $(SIGNATURES)

Write table to file.
"""
function write_table(tb :: Table, filePath :: T1) where T1 <: AbstractString
    open(filePath, "w") do io
        for lineStr in header(tb)
            write_line(io, lineStr)
        end
        for lineStr in tb.bodyV
            write_line(io, lineStr)
        end
        for lineStr in footer(tb)
            write_line(io, lineStr)
        end
    end

    # pathV = splitpath(filePath);
    d, fn = splitdir(filePath);
    println("Saved table  $fn  to dir  $d")
    return nothing
end

function write_line(io :: IO, lineStr)
    write(io, lineStr);
    if lineStr[1] == '\\'
        write(io, " \n")
    else
        # Not a command; neeed newline at end
        write(io, " \\\\ \n")
    end
end

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


text:
@inline get_data(A::AbstractArray) = A
@inline get_data(A::SizedArray) = A.data

"""
	ConVal{C,V,M,W}

Holds information about a constraint
"""
struct ConVal{C,V,M,W}
    con::C
    inds::Vector{Int}
    vals::Vector{V}
	vals2::Vector{V}
    jac::Matrix{M}
    ∇x::Matrix{W}
    ∇u::Matrix{W}
    c_max::Vector{Float64}
	is_const::Vector{Vector{Bool}}  # are the Jacobians constant
	iserr::Bool  # are the Jacobians on the error state
	function ConVal(n::Int, m::Int, con::AbstractConstraint, inds::Union{Vector{Int},UnitRange},
			jac, vals, iserr::Bool=false)
		if !iserr && size(gen_jacobian(con)) != size(jac[1])
			throw(DimensionMismatch("size of jac[i] $(size(jac[1])) does not match the expected size of $(size(gen_jacobian(con)))"))
		end
		vals2 = deepcopy(vals)
        p = length(con)
        P = length(vals)
        ix = 1:n
        iu = n .+ (1:m)
		views = [gen_views(∇c, con, n, m) for ∇c in jac]
		∇x = [v[1] for v in views]
		∇u = [v[2] for v in views]
        c_max = zeros(P)
		is_const = [zeros(Bool,P), zeros(Bool,P)]
        new{typeof(con), eltype(vals), eltype(jac), eltype(∇x)}(con,
			collect(inds), vals, vals2, jac, ∇x, ∇u, c_max, is_const, iserr)
    end
end

function ConVal(n::Int, m::Int, cval::ConVal)
	# create a ConVal for the "raw" Jacobians, if needed
	# 	otherwise return the same ConVal
	if cval.iserr
		p = length(cval.con)
		ws = widths(cval.con, n, m)
		jac = [SizedMatrix{p,w}(zeros(p,w)) for k in cval.inds, w in ws]
		ConVal(n, m, cval.con, cval.inds, jac, cval.vals, false)
	else
		return cval
	end
end

function ConVal(n::Int, m::Int, con::AbstractConstraint, inds::UnitRange{Int}, iserr::Bool=false)
	C,c = gen_convals(n,m,con,inds)
	ConVal(n, m, con, inds, C, c)
end

function _index(cval::ConVal, k::Int)
	if k ∈ cval.inds
		return k - cval.inds[1] + 1
	else
		return 0
	end
end

function evaluate!(cval::ConVal, Z::AbstractTrajectory)
	evaluate!(cval.vals, cval.con, Z, cval.inds)
end

function jacobian!(cval::ConVal, Z::AbstractTrajectory, init::Bool=false)
	if cval.iserr
		throw(ErrorException("Can't evaluate Jacobians directly on the error state Jacobians"))
	else
		jacobian!(cval.jac, cval.con, Z, cval.inds)
		# is_const = cval.is_const
	    # for (i,k) in enumerate(cval.inds)
		# 	if init || !is_const[i]
	    #     	is_const[i] = jacobian!(cval.jac[i], cval.con, Z[k])
		# 	end
	    # end
	end
end

function ∇jacobian!(G, cval::ConVal, Z::AbstractTrajectory, λ, init::Bool=false)
	∇jacobian!(G, cval.con, Z, λ, cval.inds, cval.is_const[2], init)
end

@inline violation(::Equality, v) = norm(v,Inf)
@inline violation(::Inequality, v) = maximum(v)
@inline violation(::Equality, v::Real) = abs(v)
@inline violation(::Inequality, v::Real) = v > 0 ? v : 0.0

function max_violation(cval::ConVal)
	max_violation!(cval)
    return maximum(cval.c_max)
end

function max_violation!(cval::ConVal)
	s = sense(cval.con)
    for i in eachindex(cval.inds)
        cval.c_max[i] = violation(s, cval.vals[i])
    end
end

@inline norm_violation(::Equality, v, p=2) = norm(v,p)

@inline function norm_violation(::Inequality, v, p=2)
	# TODO: try this with LazyArrays?
	if p == 1
		a = zero(eltype(v))
		for x in v
			a += max(x,0)
		end
		return a
	elseif p == 2
		a = zero(eltype(v))
		for x in v
			a += max(x, 0)^2
		end
		return sqrt(a)
	elseif p == Inf
		return maximum(v)
	else
		throw(ArgumentError("$p is not a valid norm value. Must be 1,2 or Inf"))
	end
end

function norm_violation(cval::ConVal, p=2)
	norm_violation!(cval, p)
	return norm(cval.c_max, p)
end

function norm_violation!(cval::ConVal, p=2)
	s = sense(cval.con)
	for i in eachindex(cval.inds)
		cval.c_max[i] = norm_violation(s, cval.vals[i], p)
	end
end

function norm_dgrad!(cval::ConVal, Z::AbstractTrajectory, p=1)
	for (i,k) in enumerate(cval.inds)
		zs = RobotDynamics.get_z(cval.con, Z, k)
		mul!(cval.vals2[i], cval.jac[i,1], zs[1])
		if length(zs) > 1
			mul!(cval.vals2[i], cval.jac[i,2], zs[2], 1.0, 1.0)
		end
		cval.c_max[i] = norm_dgrad(cval.vals[i], cval.vals2[i], p)
	end
	return nothing
end
"""
	dgrad(x, dx, p=1)
Directional derivative of `norm(x, p)` in the direction `dx`
"""
function norm_dgrad(x, dx, p=1)
	g = zero(eltype(x))
	if p == 1
		@assert length(x) == length(dx)
		g = zero(eltype(x))
		for i in eachindex(x)
			if x[i] < 0
				g += -dx[i]
			elseif x[i] > 0
				g += dx[i]
			else
				g += abs(dx[i])
			end
		end
	else
		throw("Directional derivative of $p-norm isn't implemented yet")
	end
	return g
end

function norm_residual!(res, cval::ConVal, λ::Vector{<:AbstractVector}, p=2)
	for (i,k) in enumerate(cval.inds)
		mul!(res[i], cval.jac[i,1], λ[i])
		if size(cval.jac,2) > 1
			mul!(res[i], cval.jac[i,2], λ[i], 1.0, 1.0)
		end
		cval.c_max[i] = norm(res[i], p)
	end
	return nothing
end

function error_expansion!(errval::ConVal, conval::ConVal, model::AbstractModel, G)
	if errval.jac !== conval.jac
		for (i,k) in enumerate(conval.inds)
			mul!(errval.∇x[i], conval.∇x[i], get_data(G[k]))
			errval.∇u[i] .= conval.∇u[i]
		end
	end
end

function error_expansion!(con::AbstractConstraint, err, jac, G)
	mul!(err, jac, G)
end

function gen_convals(n̄::Int, m::Int, con::AbstractConstraint, inds)
    # n is the state diff size
    p = length(con)
	ws = widths(con, n̄,m)
    C = [SizedMatrix{p,w}(zeros(p,w)) for k in inds, w in ws]
    c = [@MVector zeros(p) for k in inds]
    return C, c
end

function gen_convals(D::AbstractMatrix, d::AbstractVector, cinds, zinds, con::AbstractConstraint, inds)
    P = length(inds)
    p = length(con)
	n,m = get_dims(con, length(zinds[1]))
    ws = widths(con, n, m)

    C = [begin
		view(D, cinds[i], zinds[k+(j-1)][1:ws[j]])
	end for (i,k) in enumerate(inds), j = 1:length(ws)]
    c = [view(d, cinds[i]) for i in 1:P]
    return C,c
end

function gen_convals(blocks::Vector, cinds, con::AbstractConstraint, inds)
	# assumes all cinds are contiguous indices (i.e. can be represented as a UnitRange)
    C1 = map(enumerate(inds)) do (i,k)
        nm = size(blocks[k].Y,2)
		if con isa StateConstraint
			iz = 1:width(con)
		elseif con isa ControlConstraint
			m = control_dim(con)
			n = nm - m
			iz = n .+ (1:m)
		else
			iz = 1:nm
		end
		ic = cinds[i][1]:cinds[i][end]
		n1 = size(blocks[k].D2, 1)
        view(blocks[k].Y, n1 .+ (ic), iz)
    end
	C2 = map(enumerate(inds)) do (i,k)
		if con isa StageConstraint
			w = size(blocks[k].Y,2)
			view(blocks[k].Y,1:0,1:w)
		else
			w = size(blocks[k+1].Y,2)
			n = state_dim(con)
			view(blocks[k+1].Y,1:n,1:w)
		end
	end
	C = [C1 C2]
    c = map(enumerate(inds)) do (i,k)
		ic = cinds[i][1]:cinds[i][end]
        view(blocks[k].y, ic)
    end
    return C,c
end


text:
using ComparisonDataset
using Base.Test

using ComparisonDataset

@testset "Rect Tests" begin include("rect_tests.jl") end
@testset "Patch Tests" begin include("patch_tests.jl") end
@testset "Data Generation Tests" begin include("data_generation_tests.jl") end
@testset "Sample Tests" begin include("sample_tests.jl") end


text:
"""
struct Job

Struct representing a single job in the GCSPET. Consists of a transshipment to be executed on a particular location on the non-crossing axis of the gantry cranes. A job takes `t_processing` time units and can only be started as soon as `t_arrival`. Each job has a `jobtype` and `movetype`, which determine precedences between jobs and contribution to the objective function.

# Functionality (not exported)
- id(j)
- loc(j)
- t_arrival(j)
- t_processing(j)
- jobtype(j)
- movetype(j)
- istruck(j)
- istrain(j)
- isunload(j)
- isload(j)
- t_waiting(j)
- t_truckwaiting(j)
- precedes(j1, j2)
"""
struct Job
id::Int
l::Int
a::Int
p::Int
jt::Bool  # 0 = truck
mt::Bool  # 0 = unload
end

id(j::Job) = j.id
loc(j::Job) = j.l
t_arrival(j::Job) = j.a
t_processing(j::Job) = j.p
jobtype(j::Job) = j.jt
movetype(j::Job) = j.mt

istruck(j::Job) = !j.jt
istrain(j::Job) = j.jt

isunload(j::Job) = !j.mt
isload(j::Job) = j.mt

"""
    t_waiting(j::Job, t_compl)

Computes the time Job `j` spent waiting before being processed if it is completed at time `t_compl`.
"""
t_waiting(j::Job, t_compl) = t_compl - t_processing(j) - t_arrival(j)

"""
    t_truckwaiting(j::Job, t_compl)

Compute the time the (potential) truck related to Job `j` spent waiting in the terminal, if `j` is completed at time `t_compl`. This is equal to `jobtype(j) * (t_compl - t_arrival(j)`, or `jobtype(j)*(t_waiting(j, t_compl) + t_processing(j))`.
"""
t_truckwaiting(j::Job, t_compl) = istruck(j) * (t_compl - t_arrival(j))

function Base.show(io::IO, ::MIME"text/plain", j::Job)
    if istruck(j)
        a = 'V'
    else
        a = 'T'
    end
    if isload(j)
        b = 'L'
    else
        b = 'U'
    end

    pstr = join([loc(j),t_arrival(j), t_processing(j)],",")

    jstring = "$(id(j)): $pstr [$(a*b)]"
    println(io, jstring)
end

function Base.show(io::IO, ::MIME"text/plain", v::Vector{Job})
    println(io, "$(length(v))-element $(typeof(v)):")
    for j in v
        println(io, j)
    end
end

function Base.isless(j1::Job, j2::Job)
    istruck(j1) && istrain(j2) && return true
    istruck(j2) && istrain(j1) && return false
    isunload(j1) && isload(j2) && return true
    isunload(j2) && isload(j1) && return false
end

"""
    precedes(j1::Job, j2::Job)

Check if `j1` must precede `j2`. Use `≺` for short (infix) notation.
"""
function precedes(j1::Job, j2::Job)
    loc(j1) != loc(j2) && return false
    istruck(j1) && istrain(j2) && return true
    istruck(j2) && istrain(j1) && return false
    isunload(j1) && isload(j2) && return true
    isunload(j2) && isload(j1) && return false
end

"""
    ≺(j1, j2)

Shorthand for `precedes(j1,j2)`
"""
≺(j1, j2) = precedes(j1, j2)


text:
#####################################
#InnerProduct/OuterProduct###########
#####################################

immutable OuterProduct{K<:Ket,B<:Bra} <: AbstractOperator
	ket::State{K}
	bra::State{B}
end

immutable InnerProduct <: AbstractScalar
	bra::State
	ket::State
	InnerProduct{B<:Bra,K<:Ket}(b::State{B}, k::State{K}) = new(b,k)
end

promote_rule{N<:Number, I<:InnerProduct}(::Type{I}, ::Type{N}) = ScalarExpr
promote_rule{N<:Number, I<:InnerProduct, T}(::Type{Array{I,T}}, ::Type{Array{N,T}}) = Array{ScalarExpr,T}
promote_rule{N<:Number, I<:InnerProduct, T}(::Type{SparseMatrixCSC{I,T}}, ::Type{SparseMatrixCSC{N,T}}) = SparseMatrixCSC{ScalarExpr,T}

#####################################
#Misc Functions######################
#####################################

=={O<:OuterProduct}(a::O, b::O) = a.ket==b.ket && a.bra==b.bra
isequal{O<:OuterProduct}(a::O, b::O) = isequal(a.ket,b.ket) && isequal(a.bra,b.bra)

=={I<:InnerProduct}(a::I, b::I) = a.ket==b.ket && a.bra==b.bra
isequal{I<:InnerProduct}(a::I, b::I) = isequal(a.ket,b.ket) && isequal(a.bra,b.bra)

bsym(o::OuterProduct) = [bsym(o.ket), bsym(o.bra)]
label(o::OuterProduct) = [label(o.ket), label(o.bra)]
bsym(i::InnerProduct) = [bsym(i.bra), bsym(i.ket)]
label(i::InnerProduct) = [label(i.bra), label(o.ket)]
conj(i::InnerProduct) = InnerProduct(i.ket', i.bra')
ctranspose(o::OuterProduct) = OuterProduct(o.bra', o.ket')
ctranspose(i::InnerProduct) = conj(i)

isdual(a::OuterProduct, b::OuterProduct) = isdual(a.ket, b.bra) && isdual(b.ket, a.bra)
isdual(a::InnerProduct, b::InnerProduct) = isdual(a.ket, b.bra) && isdual(b.ket, a.bra)

#####################################
#Show Functions######################
#####################################

show(io::IO, o::OuterProduct) = print(io, "$(repr(o.ket))$(repr(o.bra))");
show(io::IO, i::InnerProduct) = print(io, "$(repr(i.bra))$(repr(i.ket)[2:end])");

#####################################
#Arithmetic Operations###############
#####################################

*{K<:Ket}(o::OuterProduct, s::State{K}) = dvec([inner(o.bra,s)], basis(o.ket))
*{B<:Bra}(s::State{B}, o::OuterProduct) = dvec([inner(s,o.ket)], basis(o.bra))

kron{B<:Bra}(o::OuterProduct, s::State{B}) = OuterProduct(o.ket, tensor(o.bra,s))
kron{B<:Bra}(s::State{B}, o::OuterProduct) = OuterProduct(o.ket, tensor(s,o.bra))
kron{K<:Ket}(o::OuterProduct, s::State{K}) = OuterProduct(tensor(o.ket,s), o.bra)
kron{K<:Ket}(s::State{K}, o::OuterProduct) = OuterProduct(tensor(s,o.ket), o.bra)
kron(c::DiracCoeff, o::OuterProduct) = dmat([c], basis(o.ket), basis(o.bra))
kron(o::OuterProduct, c::DiracCoeff) = kron(c,o)

kron(a::OuterProduct, b::OuterProduct) = OuterProduct(tensor(a.ket, b.ket), tensor(a.bra, b.bra))

*(a::OuterProduct, b::OuterProduct) = dmat([a.bra*b.ket], basis(a.ket), basis(b.bra))

*(c::DiracCoeff, o::OuterProduct) = kron(c,o)
*(o::OuterProduct, c::DiracCoeff) = kron(c,o)
-(o::OuterProduct) = kron(-1,o)

function +{K<:Ket, B<:Bra}(a::OuterProduct{K,B}, b::OuterProduct{K,B})
	if a==b
		return dmat([2.0], basis(a.ket), basis(b.bra))
	else
		rowb = basis(a.ket, b.ket)
		colb = basis(a.bra, b.bra)
		res = dmat(zeros(2,2)[1:length(rowb), 1:length(colb)], rowb, colb)
		res[getpos(res, a)...] = 1.0
		res[getpos(res, b)...] = 1.0
		return res
	end
end

-(a::OuterProduct, b::OuterProduct) = a+(-b)


text:
"""
    Parameters(boundary, U_field)

Struct that represents problem parameters that are determined by the geometry. A
different Parameters object should be created (or a single one updated) for every
geometry (boundary) or whenever the velocity field (U_field) changes.

Calls CalcRhoCoefs and CalcVSCoefs to update the parameters.

# ARGUMENTS
* `boundary`    : Boundary object representing geometry of boundary.
* `U_field`     : Function of the form U_field(X). Gives the velocity at every
                  point in the domain for point X.
* `sigma`       : Gaussian width. DEFAULT = 0.2.
* `verbose`     : Boolean for printing status of calculation.

# PROPERTIES
* `boundary`    : Boundary object representing geometry of boundary.
* `U_field`     : Function of the form U_field(X). Gives the velocity at every
                  point in the domain for point X.
* `alphas`      : Array containing RBF coefficients for vortex sheet strength.
* `etas`        : Array containing RBF coefficients for geometric eigen decomposition
                  function.
"""
mutable struct Parameters

    boundary::Boundary
    U_field::Union{Nothing, Function}
    alphas::Array{T} where {T<:Real}
    etas::Array{T} where {T<:Real}
    verbose::Bool

    function Parameters(boundary, U_field; sigma = 0.2, verbose = false)

        # === Calculate etas ===
        if verbose
            println("   CALCULATING ETAS IN PARAMETERS...")
        end
        etas = CalcRhoCoefs(boundary, sigma=sigma)

        # === Calculate alphas ===
        if verbose
            println("   CALCULATING ALPHAS IN PARAMETERS...")
        end
        alphas = CalcVSCoefs(boundary, etas, U_slip = U_field, sigma=sigma)

        new(boundary, U_field, alphas, etas)

    end

end


text:
"""
TolManager: convergence based only on tolerance

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""

struct TolManager <: IterationManager
    tol::Float64
    verbose::Bool
    print_skip::Int
end

TolManager(tol::Float64) = TolManager(tol, true, 50)

function Base.show(io::IO, ::MIME"text/plain", tm::TolManager)
    @printf(io,
            "TolManager: %2.4e tolerance level (%s)",
            tm.tol, tm.verbose)
end

finished(mgr::TolManager, istate::IterationState) =
    norm(istate) <= mgr.tol


text:
using Test
using RateLimiter

""" Execute a function multiple times with a limiter and return the number of seconds required to complete. """
function time_exec(limiter, cost, n)

    f() = 3

    return @elapsed for _ in 1:n
        v = @rate_limit limiter cost f()
        @assert v == 3
    end
end

#######################################################################

dt = time_exec(NoLimitRateLimiter(), 1, 1000)
println("NoLimit: dt(1)=$dt")
@assert dt < 0.001

dt = time_exec(NoLimitRateLimiter(), 2, 1000)
println("NoLimit: dt(2)=$dt")
@assert dt < 0.001

#######################################################################

tokens_per_second = 100
max_tokens = 100
initial_tokens = 0

dt = time_exec(TokenBucketRateLimiter(tokens_per_second, max_tokens, initial_tokens), 1, 1000)
println("TokenBucket: dt(1)=$dt")
@assert dt > 9.8

dt = time_exec(TokenBucketRateLimiter(tokens_per_second, max_tokens, initial_tokens), 2, 1000)
println("TokenBucket: dt(2)=$dt")
@assert dt > 19.8