text:
"""
$(SIGNATURES)

Checks for a `config.md` file in `PATHS[:src]` and uses it to set the global
variables referenced in `GLOBAL_VARS` it also sets the global latex commands
via `GLOBAL_LXDEFS`. If the configuration file is not found a warning is shown.
The keyword `init` is used internally to distinguish between the first call
where only structural variables are considered (e.g. controlling folder
structure).
"""
function process_config(; init::Bool=false)::Nothing
    if init
        # initially the paths variable aren't set, try to find a config.md
        # and read definitions in it; in particular folder_structure.
        config_path_v1 = joinpath(FOLDER_PATH[], "src", "config.md")
        config_path_v2 = joinpath(FOLDER_PATH[], "config.md")
        if isfile(config_path_v2)
            convert_md(read(config_path_v2, String); isconfig=true)
        elseif isfile(config_path_v1)
            convert_md(read(config_path_v1, String); isconfig=true)
        else
            @warn "I didn't find a config file. Ignoring."
        end
    else
        key = ifelse(FD_ENV[:STRUCTURE] < v"0.2", :src, :folder)
        dir = path(key)
        config_path = joinpath(dir, "config.md")
        if isfile(config_path)
            convert_md(read(config_path, String); isconfig=true)
        else
            @warn "I didn't find a config file. Ignoring."
        end
    end
    return nothing
end

"""
$(SIGNATURES)

Take a path to an input markdown file (via `root` and `file`), then construct
the appropriate HTML page (inserting `head`, `pg_foot` and `foot`) and finally
write it at the appropriate place.
"""
function write_page(root::String, file::String, head::String,
                    pg_foot::String, foot::String, out_path::String;
                    prerender::Bool=false, isoptim::Bool=false,
                    on_write::Function=(_,_)->nothing)::Nothing
    # 1. read the markdown into string, convert it and extract definitions
    # 2. eval the definitions and update the variable dictionary, also retrieve
    # document variables (time of creation, time of last modif) and add those
    # to the dictionary.
    fpath = joinpath(root, file)
    # The curpath is the relative path starting after /src/ so for instance:
    # f1/blah/page1.md or index.md etc... this is useful in the code evaluation
    # and management of paths
    set_cur_rpath(fpath)
    # conversion
    content = convert_md(read(fpath, String))

    # Check if should add item
    #   should we generate ? otherwise no
    #   are we in the full pass ? otherwise no
    #   is there a `rss` or `rss_description` ? otherwise no
    cond_add = GLOBAL_VARS["generate_rss"].first && # should we generate?
                FD_ENV[:FULL_PASS] &&               # are we in the full pass?
                !all(e -> isempty(locvar(e)), ("rss", "rss_description"))
    # otherwise yes
    cond_add && add_rss_item()

    # adding document variables to the dictionary
    # note that some won't change and so it's not necessary to do this every
    # time but it takes negligible time to do this so ¯\_(ツ)_/¯
    # (and it's less annoying than keeping tabs on which file has
    # already been treated etc).
    s = stat(fpath)
    set_var!(LOCAL_VARS, "fd_ctime", fd_date(unix2datetime(s.ctime)))
    set_var!(LOCAL_VARS, "fd_mtime", fd_date(unix2datetime(s.mtime)))

    # 3. process blocks in the html infra elements based on `LOCAL_VARS`
    # (e.g.: add the date in the footer)
    content = convert_html(str(content))
    head, pg_foot, foot = (e -> convert_html(e)).([head, pg_foot, foot])

    # 4. construct the page proper & prerender if needed
    pg = build_page(head, content, pg_foot, foot)
    if prerender
        # KATEX
        pg = js_prerender_katex(pg)
        # HIGHLIGHT
        if FD_CAN_HIGHLIGHT
            pg = js_prerender_highlight(pg)
            # remove script
            pg = replace(pg, r"


text:
using Test, NeuralFieldEq

@testset "Convergence space" begin
    # Test the convergence using an analytical bump solution
    # Laing, Troy, Gutkin, and Ermentrou - Multiple bumps in a neuronal model of working memory
    tol = 0.1
    function W(x)
        if x != 0
            return (((1-exp(-1.8*abs(x)))*(3.5/1.8) + (-1+exp(-1.52*abs(x)))*(3/1.52))*x)/abs(x)
        else
            return 0.0
        end
    end
    a = 2.28978
    exact_solution(x) = W(x) - W(x-a)
    I(x,t) = 0.0
    K(x) = 3.5*exp(-1.8*abs(x))-3*exp(-1.52*abs(x))
    S(V) = convert(Float64,V>0.0) # Heavyside function H(V)
    α  = 1.0
    v  = 900000.0
    V0 = exact_solution
    L  = 20
    T  = 20.0
    n  = 200
    tj = [0.2,1.0,5.0,10.0,20.0]

    in02   = Input1D(α,v,V0,L,100,T,n,I,K,S);
    prob02 = probNFE(in02)
    in01   = Input1D(α,v,V0,L,200,T,n,I,K,S);
    prob01 = probNFE(in01)
    in005  = Input1D(α,v,V0,L,400,T,n,I,K,S);
    prob005= probNFE(in005)

    V02  = solveNFE(prob02,tj)  # Solution obtained with Δx = 0.2
    V01  = solveNFE(prob01,tj)  # Solution obtained with Δx = 0.1
    V005 = solveNFE(prob005,tj) # Solution obtained with Δx = 0.05

    @test isapprox(0.0,V005(201,1),atol=0.1)   # Check if calable function is correct
    @test isapprox(0.0,V005(0.0,0.2),atol=0.1) # Check if calable function is correct

    # e_x is the absolute difference between numerical and exact solutions
    e_x = zeros(3,length(tj))
    for i = 1:length(tj)
        e_x[1,i] = abs(V02(i)[51]-exact_solution(0))
        e_x[2,i] = abs(V01(i)[101]-exact_solution(0))
        e_x[3,i] = abs(V005(i)[201]-exact_solution(0))

        @test isapprox(2.0,e_x[1,i]./e_x[2,i],atol = tol) # e_0.2/e_0.1  ≈ 2
        @test isapprox(2.0,e_x[2,i]./e_x[3,i],atol = tol) # e_0.1/e_0.05 ≈ 2
    end
end

@testset "Convergence time" begin
    # Test the convergence in time using an analytical solution obtained in
    # Pedro Lima and Evelyn Buckwar - Numerical Solution of the Neural Field Equation
    # in the two-dimensional case. Example 1, page B975
    # The exact solution is V(x,t) = exp(-t)
    tol = 0.1

    function erf(x)
        # Method to compute the error function

        signal = sign(x) # save the sign of x
        x = abs(x)

        # constants
        a1 =  0.254829592
        a2 = -0.284496736
        a3 =  1.421413741
        a4 = -1.453152027
        a5 =  1.061405429
        p  =  0.3275911

        t = 1.0/(1.0 + p*x)
        y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x)
        return signal*y # erf(-x) = -erf(x)
    end

    function I(x,y,t)
        λ = 8.0
        b = pi/λ*(erf(sqrt(λ)))^2 #-1 to 1
        return -tanh(exp(-t))*b
    end
    function K(x,y)
        λ = 8.0
        return exp(-λ*(x^2+y^2))
    end
    S(V) = tanh(V)

    α  = 1.0
    v  = 900000.0
    V0 = 1.0
    L  = 2
    N  = 200
    T  = 0.1
    tj = [0.02,0.04,0.06,0.08,0.1]

    in002   = Input2D(α,v,V0,L,N,T,5,I,K,S);
    in001   = Input2D(α,v,V0,L,N,T,10,I,K,S);
    prob002 = probNFE(in002)
    prob001 = probNFE(in001)

    V002 = solveNFE(prob002,tj) # Solution obtained with Δt = 0.02
    V001 = solveNFE(prob001,tj) # Solution obtained with Δt = 0.01

    # e_t is the absolute difference between numerical and exact solutions
    e_t = zeros(2,length(tj))
    for i = 1:length(tj)
        e_t[1,i] = abs(V002(i)[101]-exp(-tj[i]))
        e_t[2,i] = abs(V001(i)[101]-exp(-tj[i]))

        @test isapprox(2.0,e_t[1,i]./e_t[2,i],atol = tol) # e_0.02/e_0.01 ≈ 2
    end
end

@testset "Delayed NFE 2D" begin
    # Neural Field called breather
    # Example inspired in A. Hutt and N. Rougier - Numerical simulation scheme of one-
    # and two-dimensional neural fields involving space-dependent delays
    function K(x,y)
        A = 20.0/(10.0*pi)
        B = 14.0/(18.0*pi)

        return A*exp(-sqrt(x^2+y^2)) - B*exp(-sqrt(x^2+y^2)/3.0)
    end
    S(V)=convert(Float64,V>0.005) # Heavyside function H(V-Vthresh)
    I(x,y,t)=(5.0/(32.0*pi))*exp(-(x^2+y^2)/32.0)

    α  = 1.0
    v  = 20.0
    V0 = 0.0
    L  = 20
    N  = 256
    T  = 40.0
    n  = 400
    tj = collect(0:0.2:T);
    in2d = Input2D(α,v,V0,L,N,T,n,I,K,S);

    prob = probNFE(in2d)
    V    = solveNFE(prob,tj)

    # Test calable structures deterministic 2D
    @test maximum(V(40.0)) == maximum(V(201))
    @test V(0.0,0.0,40.0) == V(129,129,201)
end

@testset "Delayed SNFE 2D" begin
    # Neural Field called breather
    # Example inspired in A. Hutt and N. Rougier - Numerical simulation scheme of one-
    # and two-dimensional neural fields involving space-dependent delays
    function K(x,y)
        A = 20.0/(10.0*pi)
        B = 14.0/(18.0*pi)

        return A*exp(-sqrt(x^2+y^2)) - B*exp(-sqrt(x^2+y^2)/3.0)
    end
    S(V)=convert(Float64,V>0.005) # Heavyside function H(V-Vthresh)
    I(x,y,t)=(5.0/(32.0*pi))*exp(-(x^2+y^2)/32.0)

    α  = 1.0
    v  = 20.0
    V0 = 0.0
    L  = 20
    N  = 256
    T  = 40.0
    n  = 400
    tj = collect(0:0.2:T);
    in2d = Input2D(α,v,V0,L,N,T,n,I,K,S);

    prob = probNFE(in2d)
    Vsto = solveNFE(prob,tj,0.0005,10) # 10 trajectories, ϵ=0.0005, ξ=0.1 (default value)

    # Test calable structures stochastic 2D
    # Mean sample solution
    @test maximum(Vsto(40.0)) == maximum(Vsto(201))
    @test Vsto(0.0,0.0,40.0) == Vsto(129,129,201)

    # Trajectory 4
    @test maximum(Vsto(40.0,4)) == maximum(Vsto(201,4))
    @test Vsto(0.0,0.0,40.0,4) == Vsto(129,129,201,4)
end

@testset "Delayed NFE 1D" begin
    # Example in Kulikov, Kulikova and Lima - Numerical simulation of Neural Fields
    # with Finite Transmission Speed and Random Disturbance
    function I(x,t)
        σ  = 3.0   # External input parameter
        if t <= 1.0
            return -2.89967 + 8.0*exp(-x^2/(2.0*σ^2)) - 0.5
        else
            return -2.89967
        end
    end
    K(x) = 2*exp(-0.08*sqrt(x^2))*(0.08*sin(pi*sqrt(x^2)/10)+cos(pi*sqrt(x^2)/10))
    S(V) = convert(Float64,V>0.0) # Heavyside function H(V)

    α  = 1.0
    v  = 25.0
    V0 = 0.0
    L  = 200
    N  = 1000
    T  = 10.0
    n  = 500
    in1D = Input1D(α,v,V0,L,N,T,n,I,K,S);

    prob = probNFE(in1D)
    V    = solveNFE(prob,[1.0,5.0,10.0])

    # Test calable structures deterministic 1D
    @test maximum(V(10.0)) == maximum(V(3))
    @test V(0.0,10.0) == V(501,3)
end

@testset "Delayed SNFE 1D" begin
    # Example in Kulikov, Kulikova and Lima - Numerical simulation of Neural Fields
    # with Finite Transmission Speed and Random Disturbance
    function I(x,t)
        σ  = 3.0   # External input parameter
        if t <= 1.0
            return -2.89967 + 8.0*exp(-x^2/(2.0*σ^2)) - 0.5
        else
            return -2.89967
        end
    end
    K(x) = 2*exp(-0.08*sqrt(x^2))*(0.08*sin(pi*sqrt(x^2)/10)+cos(pi*sqrt(x^2)/10))
    S(V) = convert(Float64,V>0.0) # Heavyside function H(V)

    α  = 1.0
    v  = 15.0
    V0 = 0.0
    L  = 200
    N  = 1000
    T  = 10.0
    n  = 500
    in1D = Input1D(α,v,V0,L,N,T,n,I,K,S);

    prob = probNFE(in1D)
    Vsto = solveNFE(prob,[1.0,5.0,10.0],0.01,10) # 10 trajectories, ϵ=0.01, ξ=0.1 (default value)

    # Test calable structures stochastic 1D
    # Mean sample solution
    @test maximum(Vsto(10.0)) == maximum(Vsto(3))
    @test Vsto(0.0,10.0) == Vsto(0.0,10.0)

    # Trajectory 4
    @test maximum(Vsto(10.0,4)) == maximum(Vsto(3,4))
    @test Vsto(0.0,10.0,4) == Vsto(501,3,4)
end

ext: jl


text:
"""
    check_duplicate_reaction(rxn::Reaction, rxns::Dict{String, Reaction}; only_metabolites=true)

Check if `rxn` already exists in `rxns` but has another `id`.
If `only_metabolites` is `true` then only the metabolite `id`s are checked.
Otherwise, compares metabolite `id`s and the absolute value of their stoichiometric coefficients to those of `rxn`.
If `rxn` has the same reaction equation as another reaction in `rxns`, the return the `id`.
Otherwise return `nothing`.

See also: [`is_mass_balanced`](@ref)
"""
function check_duplicate_reaction(
    crxn::Reaction,
    rxns::OrderedDict{String,Reaction};
    only_metabolites = true,
)
    for (k, rxn) in rxns
        if rxn.id != crxn.id # skip if same ID
            if only_metabolites # only check if metabolites are the same
                if issetequal(keys(crxn.metabolites), keys(rxn.metabolites))
                    return k
                end
            else # also check the stoichiometric coefficients
                reaction_checker = true
                for (kk, vv) in rxn.metabolites # get reaction stoich
                    if abs(get(crxn.metabolites, kk, 0)) != abs(vv) # if at least one stoich doesn't match
                        reaction_checker = false
                        break
                    end
                end
                if reaction_checker &&
                   issetequal(keys(crxn.metabolites), keys(rxn.metabolites))
                    return k
                end
            end
        end
    end
    return nothing
end

"""
    is_boundary(rxn::Reaction)

Return true if reaction is a boundary reaction, otherwise return false.
Checks if boundary by inspecting number of metabolites in reaction equation.
Boundary reactions have only one metabolite, e.g. an exchange reaction, or a sink/demand reaction.
"""
function is_boundary(rxn::Reaction)::Bool
    length(keys(rxn.metabolites)) == 1 ? true : false
end

"""
    is_mass_balanced(rxn::Reaction, model::StandardModel)

Checks if `rxn` is atom balanced. Returns a boolean for whether the reaction is balanced,
and the associated balance of atoms for convenience (useful if not balanced).

See also: [`get_atoms`](@ref), [`check_duplicate_reaction`](@ref)
"""
function is_mass_balanced(rxn::Reaction, model::StandardModel)
    atom_balances = Dict{String,Float64}() # float here because stoichiometry is not Int
    for (met, stoich) in rxn.metabolites
        atoms = get_atoms(model.metabolites[met])
        isempty(atoms) && continue # ignore blanks
        for (k, v) in atoms
            atom_balances[k] = get(atom_balances, k, 0) + v * stoich
        end
    end

    return all(sum(values(atom_balances)) == 0), atom_balances
end


text:
include(joinpath(pwd(), "src/models/car.jl"))
include(joinpath(pwd(), "src/constraints/obstacles.jl"))

optimize = true

# Horizon
T = 51

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds

# ul <= u <= uu
ul, uu = control_bounds(model, T, -3.0, 3.0)

# Initial and final states
x1 = [0.0; 0.0; 0.0]
xT = [1.0; 1.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(ones(model.n)) : Diagonal(10.0 * ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T], [zeros(model.m) for t = 1:T])

# Constraints
circles = [(0.85, 0.3, 0.1),
           (0.375, 0.75, 0.1),
           (0.25, 0.2, 0.1),
           (0.75, 0.8, 0.1)]

# Constraints
function obstacles!(c, x)
    c[1] = circle_obs(x[1], x[2], circles[1][1], circles[1][2], circles[1][3])
    c[2] = circle_obs(x[1], x[2], circles[2][1], circles[2][2], circles[2][3])
    c[3] = circle_obs(x[1], x[2], circles[3][1], circles[3][2], circles[3][3])
    c[4] = circle_obs(x[1], x[2], circles[4][1], circles[4][2], circles[4][3])
    nothing
end

n_stage = 4
n_con = n_stage * T
con_obstacles = ObstacleConstraints(n_con, (1:n_con), n_stage)

# Problem
prob = trajectory_optimization_problem(model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
           con = con_obstacles
           )

# Trajectory initialization
X0 = linear_interp(x1, xT, T) # linear interpolation on state
U0 = random_controls(model, T, 0.001) # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

Z̄ = solve(prob, copy(Z0))

include(joinpath(pwd(), "src/direct_policy_optimization/dpo.jl"))

# Additive noise model
model = Car(3, 2, 3)

function fd(model::Car, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, h) - w
end

# Nominal solution
X̄, Ū = unpack(Z̄, prob)
prob_nom = prob.prob

# DPO
N = 2 * model.n
D = 2 * model.d

β = 1.0
δ = 1.0e-3

# initial samples
x1_sample = resample(x1, Diagonal([1.0; 1.0; 0.1]), 1.0e-1)

# mean problem
prob_mean = trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				h = h,
				dynamics = false
				)

# sample problems
prob_sample = [trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				h = h,
				xl = state_bounds(model, T, x1 = x1_sample[i])[1],
				xu = state_bounds(model, T, x1 = x1_sample[i])[2],
				ul = ul,
				uu = uu,
				dynamics = false,
				con = con_obstacles
				) for i = 1:N]

# sample objective
Q = [(t < T ? Diagonal([10.0; 10.0; 1.0])
	: Diagonal(100.0 * ones(model.n))) for t = 1:T]
R = [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(model.n, model.m)
dist = disturbances([Diagonal(δ * [1.0; 1.0; 0.1]) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

# TVLQR policy
K = tvlqr(model, X̄, Ū, Q, R, h)

z0_dpo = zeros(prob_dpo.num_var)
z0_dpo[prob_dpo.prob.idx.nom] = pack(X̄, Ū, prob_nom)
z0_dpo[prob_dpo.prob.idx.mean] = pack(X̄, Ū, prob_nom)
for i = 1:N
	z0_dpo[prob_dpo.prob.idx.sample[i]] = pack(X̄, Ū, prob_nom)
end
for t = 1:T-1
	z0_dpo[prob_dpo.prob.idx.policy[prob_dpo.prob.idx.θ[t]]] = vec(copy(K[t]))
end

# Solve
# Z = solve(prob_dpo, copy(z0_dpo),
# 	tol = 1.0e-3, c_tol = 1.0e-3,
# 	max_iter = 1000)

function sample_dynamics(model, xt, ut, μ, ν, w, h, t, β)
	N = 2 * model.n
	D = 2 * model.d
	M = N + D

	w0 = zeros(model.d)

	# propagate samples
	s = zeros(model.n * M)

	for j = 1:M
		if j <= N
			xi = view(xt, (j - 1) * model.n .+ (1:model.n))
			ui = view(ut, (j - 1) * model.m .+ (1:model.m))
			s[(j - 1) * model.n .+ (1:model.n)] = propagate_dynamics(model,
				xi, ui, w0, h, t)
		else
			k = j - N
			s[(j - 1) * model.n .+ (1:model.n)] = propagate_dynamics(model,
				μ, ν,
				β[t] * w[t][k], h, t)
		end
	end

	# resample
	xt⁺ = resample_vec(s, model.n, M, β[t + 1])

	return xt⁺, s
end

function sample_dynamics_jacobian(model, xt, ut, μ, ν, w, h, t, β)
	N = 2 * model.n
	D = 2 * model.d
	M = N + D

	w0 = zeros(model.d)

	dx⁺dxt = zeros(model.n * M, model.n * N)
	dx⁺dut = zeros(model.n * M, model.m * N)
	dx⁺dμ = zeros(model.n * M, model.n)
	dx⁺dν = zeros(model.n * M, model.m)

	dsdxt = zeros(model.n * M, model.n * N)
	dsdut = zeros(model.n * M, model.m * N)
	dsdμ = zeros(model.n * M, model.n)
	dsdν = zeros(model.n * M, model.m)

	xt⁺, s = sample_dynamics(model, xt, ut, μ, ν, w, h, t, β)
	r(y) = resample_vec(y, model.n, M, β[t + 1])
	dx⁺ds = real.(FiniteDiff.finite_difference_jacobian(r, s))

	# s = []
	# A = []
	# B = []

	for j = 1:M
		if j <= N
			xi = view(xt, (j - 1) * model.n .+ (1:model.n))
			ui = view(ut, (j - 1) * model.m .+ (1:model.m))
			_, _A, _B = propagate_dynamics_jacobian(model, xi, ui,
				w0, h, t)
			# push!(s, _s)
			# push!(A, _A)
			# push!(B, _B)

			dsdxt[(j - 1) * model.n .+ (1:model.n),
				(j - 1) * model.n .+ (1:model.n)] = _A
			dsdut[(j - 1) * model.n .+ (1:model.n),
				(j - 1) * model.m .+ (1:model.m)] = _B
		else
			k = j - N
			_, _A, _B = propagate_dynamics_jacobian(model, μ, ν,
				β[t] * w[t][k], h, t)

			# push!(s, _s)
			# push!(A, _A)
			# push!(B, _B)

			dsdμ[(j - 1) * model.n .+ (1:model.n), :] = _A
			dsdν[(j - 1) * model.n .+ (1:model.n), :] = _B
		end
	end

	dx⁺dxt = dx⁺ds * dsdxt
	dx⁺dut = dx⁺ds * dsdut
	dx⁺dμ = dx⁺ds * dsdμ
	dx⁺dν = dx⁺ds * dsdν

	return dx⁺dxt, dx⁺dut, dx⁺dμ, dx⁺dν
end

model = model_ft
t = 2
xt = view(z0_dpo, prob_dpo.prob.idx.xt[t])
ut = view(z0_dpo, prob_dpo.prob.idx.ut[t])
μ = view(z0_dpo, prob_dpo.prob.idx.mean[prob_dpo.prob.prob.mean.idx.x[t]])
ν = view(z0_dpo, prob_dpo.prob.idx.mean[prob_dpo.prob.prob.mean.idx.u[t]])

propagate_dynamics(model, rand(model.n), rand(model.m), rand(model.d), h, t)
propagate_dynamics_jacobian(model, rand(model.n), rand(model.m), rand(model.d), h, t)

sample_dynamics(model, xt, ut, μ, ν, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)

a1, a2, a3, a4 = sample_dynamics_jacobian(model, xt, ut, μ, ν, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)

sdx(y) = sample_dynamics(model, y, ut, μ, ν, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)[1]
sdx(xt)

sdu(y) = sample_dynamics(model, xt, y, μ, ν, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)[1]
sdu(ut)

sdμ(y) = sample_dynamics(model, xt, ut, y, ν, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)[1]
sdμ(μ)

sdν(y) = sample_dynamics(model, xt, ut, μ, y, prob_dpo.prob.dist.w, h, t,
	prob_dpo.prob.sample.β)[1]
sdν(ν)

norm(FiniteDiff.finite_difference_jacobian(sdx, xt) - a1)
norm(FiniteDiff.finite_difference_jacobian(sdu, ut) - a2)
norm(FiniteDiff.finite_difference_jacobian(sdμ, μ) - a3)
norm(FiniteDiff.finite_difference_jacobian(sdν, ν) - a4)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""
    Mattsson2014

Coefficients of the SBP operators given in
  Mattsson (2014)
  Diagonal-norm summation by parts operators for fiite difference approximations
    of third and fourth derivatives.
  Journal of Computational Physics 274, pp. 432-454.
"""
struct Mattsson2014 <: SourceOfCoefficients end

function Base.show(io::IO, ::Mattsson2014)
    print(io,
        "  Mattsson (2014) \n",
        "  Diagonal-norm summation by parts operators for fiite difference approximations\n",
        "    of third and fourth derivatives. \n",
        "  Journal of Computational Physics 274, pp. 432-454. \n")
end

function first_derivative_coefficients(source::Mattsson2014, order::Int, T=Float64, parallel=Val{:serial}())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,2}(SVector(T(-1),
                                                    T(1) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(1//2))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 1, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(-50400//35809),
                                                    T(526249//322281),
                                                    T(-75733//1933686),
                                                    T(-50767//322281),
                                                    T(-4417//71618),
                                                    T(31850//966843) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(-526249//1077057),
                                                    T(0),
                                                    T(1421209//3231171),
                                                    T(16657//239346),
                                                    T(-16934//1077057),
                                                    T(-33059//6462342) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(75733//5541372),
                                                    T(-1421209//2770686),
                                                    T(0),
                                                    T(631187//1385343),
                                                    T(400139//5541372),
                                                    T(-8789//307854) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(50767//811962),
                                                    T(-16657//180436),
                                                    T(-631187//1217943),
                                                    T(0),
                                                    T(496403//811962),
                                                    T(-308533//4871772) )),
            # d5
            DerivativeCoefficientRow{T,1,7}(SVector(T(4417//211146),
                                                    T(16934//950157),
                                                    T(-400139//5700942),
                                                    T(-496403//950157),
                                                    T(0),
                                                    T(1805647//2850471),
                                                    T(-2800//35191) )),
            # d6
            DerivativeCoefficientRow{T,1,8}(SVector(T(-31850//2713743),
                                                    T(33059//5427486),
                                                    T(8789//301527),
                                                    T(308533//5427486),
                                                    T(-1805647//2713743),
                                                    T(0),
                                                    T(22400//33503),
                                                    T(-2800//33503) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(2//3), T(-1//12))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(35809//100800),
                                T(13297//11200),
                                T(5701//5600),
                                T(45109//50400),
                                T(35191//33600),
                                T(33503//33600) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 1, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,8}(SVector(T(-508032//318365),
                                                    T(113221347//55519750),
                                                    T(-3338172//27759875),
                                                    T(-1002751721//2731571700),
                                                    T(-46815789//455261950),
                                                    T(1228638//9105239),
                                                    T(7120579//166559250),
                                                    T(-10874619//350201500) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(-4642075227//10228748530),
                                                    T(0),
                                                    T(69095487//265681780),
                                                    T(32352081//146124979),
                                                    T(28592150//438374937),
                                                    T(-137946309//1461249790),
                                                    T(-720387//44961532),
                                                    T(538846039//30686245590) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(45621684//696580885),
                                                    T(-23031829//36186020),
                                                    T(0),
                                                    T(28368209//39804622),
                                                    T(-9693137//39804622),
                                                    T(3868089//30618940),
                                                    T(-2468403//199023110),
                                                    T(-1686470//139316177) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(1002751721//11591207628),
                                                    T(-32352081//137990567),
                                                    T(-85104627//275981134),
                                                    T(0),
                                                    T(17499453//42458636),
                                                    T(13059537//275981134),
                                                    T(4924918//413971701),
                                                    T(-29088585//1931867938) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(46815789//1188140954),
                                                    T(-28592150//254601633),
                                                    T(29079411//169734422),
                                                    T(-17499453//26112988),
                                                    T(0),
                                                    T(112822635//169734422),
                                                    T(-9119079//84867211),
                                                    T(103152839//7128845724) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(-2047730//52061009),
                                                    T(45982103//371864350),
                                                    T(-3868089//57209900),
                                                    T(-4353179//74372870),
                                                    T(-7521509//14874574),
                                                    T(0),
                                                    T(474569879//743728700),
                                                    T(-138109864//1301525225),
                                                    T(4032//260045) )),
            # d7
            DerivativeCoefficientRow{T,1,10}(SVector(T(-291943739//21305233950),
                                                    T(720387//31216460),
                                                    T(7405209//1014534950),
                                                    T(-4924918//304360485),
                                                    T(9119079//101453495),
                                                    T(-1423709637//2029069900),
                                                    T(0),
                                                    T(5309800707//7101744650),
                                                    T(-108864//709465),
                                                    T(12096//709465) )),
            # d8
            DerivativeCoefficientRow{T,1,11}(SVector(T(10874619//1121684300),
                                                    T(-538846039//21872843850),
                                                    T(1011882//145818959),
                                                    T(5817717//291637918),
                                                    T(-103152839//8749137540),
                                                    T(414329592//3645473975),
                                                    T(-5309800707//7290947950),
                                                    T(0),
                                                    T(762048//1019713),
                                                    T(-762048//5098565),
                                                    T(84672//5098565) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(318365//1016064),
                                T(145979//103680),
                                T(139177//241920),
                                T(964969//725760),
                                T(593477//725760),
                                T(52009//48384),
                                T(141893//145152),
                                T(1019713//1016064) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 1, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

function second_derivative_coefficients(source::Mattsson2014, order::Int, T=Float64, parallel=Val{:serial}())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(1))
        central_coef = T(-2)
        lower_coef = upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
                                                    T(2),
                                                    T(-1//2) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 2, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(8027765//3867372),
                                                    T(-1690702//322281),
                                                    T(8240267//1933686),
                                                    T(-1030334//966843),
                                                    T(-6817//143236),
                                                    T(21380//966843) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(1030898//1077057),
                                                    T(-23780195//12924684),
                                                    T(2540018//3231171),
                                                    T(26257//239346),
                                                    T(-12268//3231171),
                                                    T(-119459//12924684) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(75467//5541372),
                                                    T(1270009//1385343),
                                                    T(-558115//307854),
                                                    T(1111174//1385343),
                                                    T(551339//5541372),
                                                    T(-8789//461781) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(-61567//1217943),
                                                    T(26257//180436),
                                                    T(1111174//1217943),
                                                    T(-5126635//2435886),
                                                    T(464003//405981),
                                                    T(-222133//4871772) )),
            # d5
            DerivativeCoefficientRow{T,1,7}(SVector(T(-6817//422292),
                                                    T(-12268//2850471),
                                                    T(551339//5700942),
                                                    T(928006//950157),
                                                    T(-25370795//11401884),
                                                    T(3568094//2850471),
                                                    T(-2800//35191) )),
            # d6
            DerivativeCoefficientRow{T,1,8}(SVector(T(21380//2713743),
                                                    T(-119459//10854972),
                                                    T(-17578//904581),
                                                    T(-222133//5427486),
                                                    T(3568094//2713743),
                                                    T(-9063745//3618324),
                                                    T(44800//33503),
                                                    T(-2800//33503) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(4//3), T(-1//12))
        central_coef = T(-5//2)
        lower_coef = upper_coef
        left_weights = SVector( T(35809//100800),
                                T(13297//11200),
                                T(5701//5600),
                                T(45109//50400),
                                T(35191//33600),
                                T(33503//33600) )
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6),
                                                    T(3),
                                                    T(-3//2),
                                                    T(1//3) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 2, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,8}(SVector(T(24055498439//8194715100),
                                                    T(-9925742373//1138154875),
                                                    T(983223468//103468625),
                                                    T(-18528585641//4097357550),
                                                    T(16391451//22207900),
                                                    T(13725804//227630975),
                                                    T(355447739//20486787750),
                                                    T(-24731721//2276309750) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(270821931//300845545),
                                                    T(-42416226217//26302496220),
                                                    T(685962357//1461249790),
                                                    T(2160993//8595587),
                                                    T(59905900//1315124811),
                                                    T(-158509509//2922499580),
                                                    T(-9667431//1461249790),
                                                    T(634102039//92058736770) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(-2479644//63325535),
                                                    T(228654119//199023110),
                                                    T(-16197861//7237204),
                                                    T(24739409//19902311),
                                                    T(-7878737//39804622),
                                                    T(1917829//18093010),
                                                    T(-7508403//398046220),
                                                    T(-41036//10716629) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(-1092927401//17386811442),
                                                    T(36736881//137990567),
                                                    T(74218227//137990567),
                                                    T(-7780367599//4967660412),
                                                    T(210256089//275981134),
                                                    T(1500627//25089194),
                                                    T(1246172//95531931),
                                                    T(-37555785//3863735876) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(-54436269//2376281908),
                                                    T(59905900//763804899),
                                                    T(-23636211//169734422),
                                                    T(210256089//169734422),
                                                    T(-7116321131//3055219596),
                                                    T(8176215//6528247),
                                                    T(-7304679//84867211),
                                                    T(84101639//10693268586) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(4575268//260305045),
                                                    T(-52836503//743728700),
                                                    T(1917829//33805850),
                                                    T(500209//6761170),
                                                    T(545081//572099),
                                                    T(-324760747//148745740),
                                                    T(461969879//371864350),
                                                    T(-11753624//118320475),
                                                    T(2688//260045) )),
            # d7
            DerivativeCoefficientRow{T,1,10}(SVector(T(355447739//63915701850),
                                                    T(-9667431//1014534950),
                                                    T(-22525209//2029069900),
                                                    T(1246172//70237035),
                                                    T(-7304679//101453495),
                                                    T(1385909637//1014534950),
                                                    T(-48284442317//18261629100),
                                                    T(5288632707//3550872325),
                                                    T(-108864//709465),
                                                    T(8064//709465) )),
            # d8
            DerivativeCoefficientRow{T,1,11}(SVector(T(-24731721//7290947950),
                                                    T(634102039//65618531550),
                                                    T(-123108//56084215),
                                                    T(-7511157//583275836),
                                                    T(84101639//13123706310),
                                                    T(-35260872//331406725),
                                                    T(5288632707//3645473975),
                                                    T(-70820489957//26247412620),
                                                    T(1524096//1019713),
                                                    T(-762048//5098565),
                                                    T(56448//5098565) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(3//2), T(-3//20), T(1//90))
        central_coef = T(-49//18)
        lower_coef = upper_coef
        left_weights = SVector( T(318365//1016064),
                                T(145979//103680),
                                T(139177//241920),
                                T(964969//725760),
                                T(593477//725760),
                                T(52009//48384),
                                T(141893//145152),
                                T(1019713//1016064) )
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,5}(SVector(T(-25/12),
                                                    T(4),
                                                    T(-3),
                                                    T(4//3),
                                                    T(-1//4) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 2, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

function third_derivative_coefficients(source::Mattsson2014, order::Int, T=Float64, parallel=Val{:serial}())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,4}(SVector(T(1//4),
                                                    T(-5//8),
                                                    T(1//2),
                                                    T(-1//8) )),
            # d2
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11//16),
                                                    T(2),
                                                    T(-31//16),
                                                    T(5//8) )),
            # d3
            DerivativeCoefficientRow{T,1,5}(SVector(T(-1//2),
                                                    T(15//16),
                                                    T(1//8),
                                                    T(-17//16),
                                                    T(1//2) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(-1), T(1//2))
        central_coef = T(0)
        lower_coef = -upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
                                                    T(2),
                                                    T(-1//2) )),
            # second derivative
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2]
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 3, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(-32200//35809),
                                                    T(188187//71618),
                                                    T(-90183//35809),
                                                    T(27988//35809),
                                                    T(-801//35809),
                                                    T(2205//71618) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(-96329//79782),
                                                    T(50400//13297),
                                                    T(-163046//39891),
                                                    T(63583//39891),
                                                    T(-1337//26594),
                                                    T(-1567//39891) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(-11939//34206),
                                                    T(5923//17103),
                                                    T(6300//5701),
                                                    T(-32543//17103),
                                                    T(29083//34206),
                                                    T(-284//5701) )),
            # d4
            DerivativeCoefficientRow{T,1,7}(SVector(T(5606//45109),
                                                    T(-89949//90218),
                                                    T(72429//45109),
                                                    T(2800//45109),
                                                    T(-77319//45109),
                                                    T(95517//90218),
                                                    T(-6300//45109) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(267//35191),
                                                    T(4011//70382),
                                                    T(-29083//35191),
                                                    T(51546//35191),
                                                    T(0),
                                                    T(-108271//70382),
                                                    T(33600//35191),
                                                    T(-4200//35191) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(-735//67006),
                                                    T(1567//33503),
                                                    T(1704//33503),
                                                    T(-31839//33503),
                                                    T(108271//67006),
                                                    T(0),
                                                    T(-54600//33503),
                                                    T(33600//33503),
                                                    T(-4200//33503) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(-13//8), T(1), T(-1//8))
        central_coef = T(0)
        lower_coef = -upper_coef
        left_weights = SVector( T(35809//100800),
                                T(13297//11200),
                                T(5701//5600),
                                T(45109//50400),
                                T(35191//33600),
                                T(33503//33600) )
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6),
                                                    T(3),
                                                    T(-3//2),
                                                    T(1//3) )),
            # second derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(2),
                                                    T(-5),
                                                    T(4),
                                                    T(-1) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2]
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 3, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,8}(SVector(T(-151704//63673),
                                                    T(15270676769//1821047800),
                                                    T(-443349971//41387450),
                                                    T(2063356637//364209560),
                                                    T(-39300617//45526195),
                                                    T(-11473393//364209560),
                                                    T(-38062741//455261950),
                                                    T(40315779//1821047800) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(-13333381409//8182998824),
                                                    T(829440//145979),
                                                    T(-8702160983//1168999832),
                                                    T(1321219979//292249958),
                                                    T(-1463113021//1168999832),
                                                    T(1240729//20874997),
                                                    T(102110955//1168999832),
                                                    T(-50022767//2045749706) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(14062931//75990642),
                                                    T(-1261072297//477655464),
                                                    T(1088640//139177),
                                                    T(-4530616889//477655464),
                                                    T(602572103//119413866),
                                                    T(-116503713//159218488),
                                                    T(-17846623//59706933),
                                                    T(343537955//3343588248) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(661223855//7727471752),
                                                    T(-214194059//275981134),
                                                    T(1209539129//1103924536),
                                                    T(645120//964969),
                                                    T(-2321979501//1103924536),
                                                    T(327603877//275981134),
                                                    T(-175223717//1103924536),
                                                    T(1353613//965933969) )),
            # d5
            DerivativeCoefficientRow{T,1,9}(SVector(T(-91064195//594070477),
                                                    T(632843581//678937688),
                                                    T(-446896583//169734422),
                                                    T(2045223021//678937688),
                                                    T(22680//593477),
                                                    T(-1804641793//678937688),
                                                    T(311038417//169734422),
                                                    T(-1932566239//4752563816),
                                                    T(21168//593477) )),
            # d6
            DerivativeCoefficientRow{T,1,10}(SVector(T(11473393//1249464216),
                                                    T(-8685103//111559305),
                                                    T(116503713//297491480),
                                                    T(-327603877//223118610),
                                                    T(1804641793//892474440),
                                                    T(0),
                                                    T(-1760949511//892474440),
                                                    T(2105883973//1561830270),
                                                    T(-72576//260045),
                                                    T(7056//260045) )),
            # d7
            DerivativeCoefficientRow{T,1,11}(SVector(T(38062741//1420348930),
                                                    T(-20422191//162325592),
                                                    T(17846623//101453495),
                                                    T(175223717//811627960),
                                                    T(-311038417//202906990),
                                                    T(1760949511//811627960),
                                                    T(0),
                                                    T(-1081094773//516490520),
                                                    T(1022112//709465),
                                                    T(-217728//709465),
                                                    T(21168//709465) )),
            # d8
            DerivativeCoefficientRow{T,1,12}(SVector(T(-40315779//5832758360),
                                                    T(50022767//1458189590),
                                                    T(-68707591//1166551672),
                                                    T(-1353613//729094795),
                                                    T(1932566239//5832758360),
                                                    T(-2105883973//1458189590),
                                                    T(1081094773//530250760),
                                                    T(0),
                                                    T(-10329984//5098565),
                                                    T(7154784//5098565),
                                                    T(-1524096//5098565),
                                                    T(148176//5098565) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(-61//30), T(169//120), T(-3//10), T(7//240))
        central_coef = T(0)
        lower_coef = -upper_coef
        left_weights = SVector( T(318365//1016064),
                                T(145979//103680),
                                T(139177//241920),
                                T(964969//725760),
                                T(593477//725760),
                                T(52009//48384),
                                T(141893//145152),
                                T(1019713//1016064) )
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12),
                                                    T(4),
                                                    T(-3),
                                                    T(4//3),
                                                    T(-1//4) )),
            # second derivative
            DerivativeCoefficientRow{T,1,5}(SVector(T(35//12),
                                                    T(-26//3),
                                                    T(19//2),
                                                    T(-14//3),
                                                    T(11//12) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2]
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 3, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

function fourth_derivative_coefficients(source::Mattsson2014, order::Int, T=Float64, parallel=Val{:serial}())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,4}(SVector(T(8//5),
                                                    T(-24//5),
                                                    T(24//5),
                                                    T(-8//5) )),
            # d2
            DerivativeCoefficientRow{T,1,4}(SVector(T(-2//5),
                                                    T(6//5),
                                                    T(-6//5),
                                                    T(2//5) )),
            # d3
            DerivativeCoefficientRow{T,1,5}(SVector(T(2//5),
                                                    T(-11//5),
                                                    T(21//5),
                                                    T(-17//5),
                                                    T(1) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(1//5),
                                                    T(2//5),
                                                    T(-17//5),
                                                    T(-29//5),
                                                    T(-4),
                                                    T(1) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(-4), T(1))
        central_coef = T(6)
        lower_coef = upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
                                                    T(2),
                                                    T(-1//2) )),
            # second derivative
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
            # third derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(-1),
                                                    T(3),
                                                    T(-3),
                                                    T(1) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2],
            -left_boundary_derivatives[3],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 4, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(-242219//644562),
                                                    T(881057//644562),
                                                    T(-183673//107427),
                                                    T(220981//322281),
                                                    T(109057//644562),
                                                    T(-29273//214854) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(578657//2154114),
                                                    T(-703457//718038),
                                                    T(1327457//1077057),
                                                    T(-544543//1077057),
                                                    T(-79457//718038),
                                                    T(204257//2154114) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(219527//307854),
                                                    T(-2754943//923562),
                                                    T(2216981//461781),
                                                    T(-559673//153927),
                                                    T(1141057//923562),
                                                    T(-120619//923562) )),
            # d4
            DerivativeCoefficientRow{T,1,7}(SVector(T(69781//811962),
                                                    T(665057//811962),
                                                    T(-584873//135327),
                                                    T(2995381//405981),
                                                    T(-4614143//811962),
                                                    T(172109//90218),
                                                    T(-8400//45109) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(8389//146178),
                                                    T(-79457//633438),
                                                    T(1141057//950157),
                                                    T(-4614143//950157),
                                                    T(557127//70382),
                                                    T(-11293343//1900314),
                                                    T(67200//35191),
                                                    T(-5600//35191) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(-29273//603054),
                                                    T(204257//1809162),
                                                    T(-120619//904581),
                                                    T(172109//100509),
                                                    T(-11293343//1809162),
                                                    T(16787381//1809162),
                                                    T(-218400//33503),
                                                    T(67200//33503),
                                                    T(-5600//33503) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(-13//2), T(2), T(-1//6))
        central_coef = T(28//3)
        lower_coef = upper_coef
        left_weights = SVector( T(35809//100800),
                                T(13297//11200),
                                T(5701//5600),
                                T(45109//50400),
                                T(35191//33600),
                                T(33503//33600) )
        right_weights = left_weights
        left_boundary_derivatives = (
            # first derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6),
                                                    T(3),
                                                    T(-3//2),
                                                    T(1//3) )),
            # second derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(2),
                                                    T(-5),
                                                    T(4),
                                                    T(-1) )),
            # third derivative
            DerivativeCoefficientRow{T,1,4}(SVector(T(-1),
                                                    T(3),
                                                    T(-3),
                                                    T(1) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2],
            -left_boundary_derivatives[3],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 4, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,8}(SVector(T(37567391168//53948541075),
                                                    T(-95834307667//35965694050),
                                                    T(252350074//65392171),
                                                    T(-58232913019//21579416430),
                                                    T(4040770588//3596569405),
                                                    T(-15248255797//35965694050),
                                                    T(4832196698//53948541075),
                                                    T(134156001//7193138810) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(29125918379//23087746682),
                                                    T(-1255810938848//242421340161),
                                                    T(1289206067431//161614226774),
                                                    T(-431078362378//80807113387),
                                                    T(494586219497//484842680322),
                                                    T(31446420748//80807113387),
                                                    T(-21701585799//161614226774),
                                                    T(334788562//242421340161) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(1308658570//3001630359),
                                                    T(-88210933529//66035867898),
                                                    T(13622370452//11005977983),
                                                    T(-27138341627//66035867898),
                                                    T(23881355534//33017933949),
                                                    T(-26412188989//22011955966),
                                                    T(21399717536//33017933949),
                                                    T(-928716467//9433695414) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(110582060185//457852701306),
                                                    T(-22954806538//76308783551),
                                                    T(-180184675067//152617567102),
                                                    T(678091654628//228926350653),
                                                    T(-378329435643//152617567102),
                                                    T(69519106966//76308783551),
                                                    T(-98928859751//457852701306),
                                                    T(4720003312//76308783551) )),
            # d5
            DerivativeCoefficientRow{T,1,9}(SVector(T(1870177580//46931567683),
                                                    T(-21945155863//281589406098),
                                                    T(45403496174//46931567683),
                                                    T(-384706366203//93863135366),
                                                    T(974238057544//140794703049),
                                                    T(-520477408939//93863135366),
                                                    T(99162460006//46931567683),
                                                    T(-99640101991//281589406098),
                                                    T(21168//593477) )),
            # d6
            DerivativeCoefficientRow{T,1,10}(SVector(T(-15248255797//123384591330),
                                                    T(31446420748//61692295665),
                                                    T(-26412188989//41128197110),
                                                    T(69519106966//61692295665),
                                                    T(-520477408939//123384591330),
                                                    T(155376599432//20564098555),
                                                    T(-772894368601//123384591330),
                                                    T(21159425698//8813185095),
                                                    T(-96768//260045),
                                                    T(7056//260045) )),
            # d7
            DerivativeCoefficientRow{T,1,11}(SVector(T(690313814//24044478315),
                                                    T(-21701585799//112207565470),
                                                    T(21399717536//56103782735),
                                                    T(-98928859751//336622696410),
                                                    T(99162460006//56103782735),
                                                    T(-772894368601//112207565470),
                                                    T(1826861184956//168311348205),
                                                    T(-915425403107//112207565470),
                                                    T(2044224//709465),
                                                    T(-290304//709465),
                                                    T(21168//709465) )),
            # d8
            DerivativeCoefficientRow{T,1,12}(SVector(T(134156001//23039395522),
                                                    T(334788562//172795466415),
                                                    T(-6501015269//115196977610),
                                                    T(4720003312//57598488805),
                                                    T(-99640101991//345590932830),
                                                    T(148115979886//57598488805),
                                                    T(-915425403107//115196977610),
                                                    T(1952118169516//172795466415),
                                                    T(-41319936//5098565),
                                                    T(14309568//5098565),
                                                    T(-2032128//5098565),
                                                    T(148176//5098565) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(-122//15), T(169//60), T(-2//5), T(7//240))
        central_coef = T(91//8)
        lower_coef = upper_coef
        left_weights = SVector( T(318365//1016064),
                                T(145979//103680),
                                T(139177//241920),
                                T(964969//725760),
                                T(593477//725760),
                                T(52009//48384),
                                T(141893//145152),
                                T(1019713//1016064) )
        right_weights = left_weights
        left_boundary_derivatives = (
        # first derivative
        DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12),
                                                T(4),
                                                T(-3),
                                                T(4//3),
                                                T(-1//4) )),
        # second derivative
        DerivativeCoefficientRow{T,1,5}(SVector(T(35//12),
                                                T(-26//3),
                                                T(19//2),
                                                T(-14//3),
                                                T(11//12) )),
        # third derivative
        DerivativeCoefficientRow{T,1,5}(SVector(T(-5//2),
                                                T(9),
                                                T(-12),
                                                T(7),
                                                T(-3//2) )),
        )
        right_boundary_derivatives = (
            -left_boundary_derivatives[1],
            left_boundary_derivatives[2],
            -left_boundary_derivatives[3],
        )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, parallel, 4, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end


text:
# This file contains an Example Lesson plan for Julia tutor.

prompt_1 =
"""Welcome to your first day of class. Today we'll be doing a classic excercise. We will be learning how to print text to our console. Along the way we'll cover some ways to handle strings in Julia.
Strings are defined by double quotes " " in julia. For example, we could define a variable, called console_output to contain the text "hello world" like the following,
console_output = "hello world"

Okay now it's your turn,
"""

request_1 = """Define a string titled 'console_output' to contain the text "hello world" """

evaluator_1 = ExpressionEvaluator(:console_output, "hello world")

hint_1 = """Please ensure your solution is cases sensitive.
IE: "HeLlO WoRlD" is not the same as "hello world". """

# These four items compose a single lesson instance.
lesson_1 = Lesson(prompt_1, request_1, evaluator_1, hint_1)

prompt_2 = """Congratulations you have instantiated your first variable in Julia! How might we display this text to console? There are two functions that let you do this: print() and println(). Print will display the String or variable as it is, but println will add a new line/carriage return to the end of it.

Now that our console_output variable is stored in memory, we can print it as follows:
print(console_output)
"""

request_2 = """Print your 'console_output' variable to the console using the `print()` function"""

evaluator_2 = CommandEvaluator("hello world")
# Notice: we did not provide any hints for this excercise. Not nice, but either they get it, or they don't.
lesson_2 = Lesson(prompt_2, request_2, evaluator_2)

# To create the lesson_plan (which is REQUIRED) we make an array of Lesson's.
lesson_plan = [lesson_1, lesson_2]


text:
######################################################################
# Additional errors used in the library.
# -----
# Licensed under MIT License

export CancellationError
struct CancellationError <: Exception
    what
end
CancellationError() = CancellationError(nothing)

function Base.showerror(io::IO, err::CancellationError)
    print(io, "Cancelled")
    if err.what !== nothing
        print(io, ": ")
        print(io, err.what)
    end
end

export TimeoutError
struct TimeoutError <: Exception
    message::String
end
TimeoutError() = TimeoutError("")

function Base.showerror(io::IO, err::TimeoutError)
    if isempty(err.message)
        print(io, "Timeout")
    else
        print(io, "Timeout: $(err.message)")
    end
end


text:
#Base.include(@__MODULE__,"../basis/BasisStructs.jl")
"""
  module JCInput
The module required for reading in and processing the selected input file.
Import this module into the script when you need to process an input file
(which will be every single calculation).
"""
module JCBasis

using JuliaChem.JCModules
using JuliaChem.JERI

using CxxWrap
using MPI
using Base.Threads
using HDF5
using Printf

Base.include(@__MODULE__, "BasisHelpers.jl")

"""
  run(args::String)
Perform the operations necessary to read in, process, and extract data from the
selected input file.

One input variable is required:
1. args = The name of the input file.

Two variables are output:
1. input_info = Information gathered from the input file.
2. basis = The basis set shells, determined from the input file.

Thus, proper use of the Input.run() function would look like this:

```
input_info, basis = Input.run(args)
```
"""
function run(molecule, model; output="none")
  comm=MPI.COMM_WORLD

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println("--------------------------------------------------------------------------------")
    println("                       ========================================                 ")
    println("                                GENERATING BASIS SET                            ")
    println("                       ========================================                 ")
    println(" ")
  end

  #== initialize variables ==#
  geometry_array::Vector{Float64} = molecule["geometry"]
  symbols::Vector{String} = molecule["symbols"]
  basis::String = model["basis"]
  charge::Int64 = molecule["molecular_charge"]

  num_atoms::Int64 = length(geometry_array)/3
  geometry_array_t::Matrix{Float64} = reshape(geometry_array,(3,num_atoms))
  geometry::Matrix{Float64} = transpose(geometry_array_t)
  geometry .*= 1.0/0.52917724924 #switch from angs to bohr

  atomic_number_mapping::Dict{String,Int64} = create_atomic_number_mapping()
  atomic_mass_mapping::Dict{String,Float64} = create_atomic_mass_mapping()
  shell_am_mapping::Dict{String,Int64} = create_shell_am_mapping()

  mol = Molecule([], StdVector{JERI.Atom}())

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println("----------------------------------------          ")
    println("          Printing basis set...                   ")
    println("----------------------------------------          ")
    println()
  end

  basis_set_shells = Vector{JCModules.Shell}([])
  shells_cxx = StdVector([ StdVector{JERI.Shell}() for i in 1:55 ])
  shells_cxx_added = [ false for i in 1:55 ]

  basis_set_nels = -charge
  basis_set_norb = 0
  pos = 1

  #== relocate center of mass of system to origin ==#
  center_of_mass = Vector{Float64}([0.0, 0.0, 0.0])

  atom_centers = Vector{Vector{Float64}}([])
  atomic_masses = Vector{Float64}([])
  for (atom_idx, symbol) in enumerate(symbols)
    push!(atom_centers, geometry[atom_idx,:])
    push!(atomic_masses, atomic_mass_mapping[symbol])

    center_of_mass .+= atomic_masses[atom_idx] .* atom_centers[atom_idx]
  end
  center_of_mass ./= sum(atomic_masses)

  for icenter in atom_centers
    icenter .-= center_of_mass
  end

  #== create basis set ==#
  h5open(joinpath(@__DIR__, "../../records/bsed.h5"),"r") do bsed
    shell_id = 1

    for (atom_idx, symbol) in enumerate(symbols)
      #== initialize variables needed for shell ==#
      atom_center = atom_centers[atom_idx]
      atomic_number = atomic_number_mapping[symbol]

      #== create atom objects ==#
      push!(mol, JCModules.Atom(atomic_number, symbol, atom_center))
      push!(mol.mol_cxx, JERI.create_atom(atomic_number, atom_center))

      basis_set_nels += atomic_number

      #== read in basis set values==#
      shells::Dict{String,Any} = read(
        bsed["$symbol/$basis"])

      #== process basis set values into shell objects ==#
      if MPI.Comm_rank(comm) == 0 && output >= 2
        println("Atom #$atom_idx ($symbol):")
        println("  Shell     Ang. Mom.     Prim.       Exp.         Coeff.")
      end

      for shell_num::Int64 in 1:length(shells)
        new_shell_dict::Dict{String,Any} = shells["$shell_num"]

        new_shell_am::Int64 = shell_am_mapping[new_shell_dict["Shell Type"]]

        new_shell_exp = new_shell_dict["Exponents"]
        nprim = length(new_shell_exp)

        new_shell_coeff = new_shell_dict["Coefficients"]
        ncoeff = length(new_shell_coeff)

        #== if L shell, divide up ==#
        if new_shell_am == -1
          #== s component ==#
          if MPI.Comm_rank(comm) == 0 && output >= 2
            for iprim in 1:nprim
              @printf("    %d        L (s)          %d     %.6f     %.6f\n",
                shell_num, iprim, new_shell_exp[iprim], new_shell_coeff[iprim])
            end
          end

          new_shell = JCModules.Shell(shell_id, atom_idx, atomic_number,
            new_shell_exp, new_shell_coeff[1:nprim],
            atom_center, 1, nprim, pos, true)
          push!(basis_set_shells, new_shell)
          #display(new_shell_coeff[:,1])
          if !shells_cxx_added[atomic_number+1]
            push!(shells_cxx[atomic_number+1], JERI.create_shell(0,
              new_shell_exp, new_shell_coeff[1:nprim], atom_center))
          end

          basis_set_norb += 1
          shell_id += 1
          pos += new_shell.nbas

          #== p component ==#
          if MPI.Comm_rank(comm) == 0 && output >= 2
            for iprim in 1:nprim
              @printf("    %d        L (p)          %d      %.6f     %.6f\n",
                shell_num, iprim, new_shell_exp[iprim],
                new_shell_coeff[nprim+iprim])
            end
          end

          new_shell = JCModules.Shell(shell_id, atom_idx, atomic_number,
            new_shell_exp, new_shell_coeff[(nprim+1):ncoeff],
            atom_center, 2, nprim, pos, true)
          push!(basis_set_shells,new_shell)
          #display(new_shell_coeff[:,2])
          if !shells_cxx_added[atomic_number+1]
            push!(shells_cxx[atomic_number+1], JERI.create_shell(1,
              new_shell_exp, new_shell_coeff[(nprim+1):ncoeff], atom_center))
          end

          basis_set_norb += 3
          shell_id += 1
          pos += new_shell.nbas
        #== otherwise accept shell as is ==#
        else
          if MPI.Comm_rank(comm) == 0 && output >= 2
            for iprim in 1:nprim
              @printf("    %d          %s            %d      %.6f     %.6f\n",
                shell_num, new_shell_dict["Shell Type"], iprim,
                new_shell_exp[iprim], new_shell_coeff[iprim])
            end
          end

          new_shell = JCModules.Shell(shell_id, atom_idx, atomic_number,
            new_shell_exp, deepcopy(new_shell_coeff),
            atom_center, new_shell_am, nprim, pos, true)
          push!(basis_set_shells,new_shell)
          if !shells_cxx_added[atomic_number+1]
            push!(shells_cxx[atomic_number+1], JERI.create_shell(new_shell_am-1,
              new_shell_exp, new_shell_coeff, atom_center))
          end

          basis_set_norb += new_shell.nbas
          shell_id += 1
          pos += new_shell.nbas
        end

        if MPI.Comm_rank(comm) == 0 && output >= 2
          println(" ")
        end
      end

      shells_cxx_added[atomic_number+1] = true
      #display(shells_cxx)

      if MPI.Comm_rank(comm) == 0 && output >= 2
        println(" ")
      end
    end
  end

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println("----------------------------------------          ")
    println("     Printing basis set metadata...               ")
    println("----------------------------------------          ")
    println()
    println("Basis set: $basis")
    println("Number of basis functions: $basis_set_norb")
    println("Number of electrons: $basis_set_nels")
  end

  #sort!(basis_set_shells, by = x->((x.nbas*x.nprim),x.am))
  #sort!(basis_set_shells, by = x->(x.atomic_number,x.atom_id))

  basis_set_cxx = JERI.BasisSet(mol.mol_cxx, shells_cxx)
  basis_set::Basis = Basis(basis_set_shells, basis_set_cxx,
    StdVector{JERI.ShellPair}(), basis,
    basis_set_norb, basis_set_nels)

  precompute_shell_pair_data(basis_set.shpdata_cxx, basis_set.basis_cxx)

  #== set up shell pair ordering ==#
  #for ish in 1:length(basis_set.shells), jsh in 1:ish
  #  push!(basis_set.shpair_ordering, ShPair(basis_set.shells[ish],
  #    basis_set.shells[jsh]))
  #end

  #sort!(basis_set.shpair_ordering, by = x->((x.nbas2*x.nprim2),x.am2,unsafe_string(x.class)))

  #for ish in 1:length(basis_set.shells), jsh in 1:ish
  #  idx = ceil(Int64, ish*(ish-1)/2) + jsh
  #  push!(basis_set.shpair_ordering,(shellpairs[idx].sh_a.shell_id,
  #    shellpairs[idx].sh_b.shell_id))
  #end

  #for shellpair in shellpairs
  #  basis_set.shpair_ordering = vcat(basis_set.shpair_ordering, [shellpair.sh_a.shell_id shellpair.sh_b.shell_id])
  #end

  #delete first row, as it is simply zeroes
  #basis_set.shpair_ordering = basis_set.shpair_ordering[setdiff(1:end, 1),:]

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println(" ")
    println("                       ========================================                 ")
    println("                                       END BASIS                                ")
    println("                       ========================================                 ")
    println("--------------------------------------------------------------------------------")
  end

  return mol, basis_set
end
export run

end


text:
using Documenter, SpatialEcology

makedocs(
    sitename = "Spatial Ecology in Julia",
    authors = "Michael Krabbe Borregaard",
    pages = [
        "Home" => "index.md",
        "Tutorial / Quick start" => "tutorial.md",
        "Manual" => Any[
            "Creating SpatialEcology objects" => "man/constructors.md",
            "Working with SpatialEcology objects" => "man/getters.md",
            "Getting and setting data" => "man/data.md",
            "Subsetting" => "man/subsetting.md",
            "Randomizations" => "man/randomization.md",
        ],
        "Example analyses" => Any[
            "Node-based analysis of species distributions" => "examples/nodebased.md"
        ],
        "Library" => "lib/public.md"
    ]
)

deploydocs(
    repo = "github.com/EcoJulia/SpatialEcology.jl.git",
    push_preview = true
)


text:
##########################################################################################
#
#    Log likelihood accumulator : type definition and derivation rules
#
##########################################################################################

# this makes the model function easier to generate compared to a Float64
#   - embeds the error throwing when log-likelihood reaches -Inf
#   - calculates the sum when logpdf() returns an Array
type OutOfSupportError <: Exception ; end

immutable LLAcc
	val::Float64
	function LLAcc(x::Real)
		isfinite(x) || throw(OutOfSupportError())
		new(x)
	end
end
+(ll::LLAcc, x::Real)           = LLAcc(ll.val + x)
+(ll::LLAcc, x::Array{Float64}) = LLAcc(ll.val + sum(x))

ReverseDiffSource.declareType(Lora.LLAcc, :LLAcc) # declares new type to Autodiff

####### derivation rules  ############
# (note : only additions are possible with LLAcc type )
ReverseDiffSource.@deriv_rule getfield(x::LLAcc, f      )      x     dx1 = ds

ReverseDiffSource.@deriv_rule +(x::LLAcc, y      )             x     dx1 += ds1
ReverseDiffSource.@deriv_rule +(x::LLAcc, y::Real)             y     dy += ds1
ReverseDiffSource.@deriv_rule +(x::LLAcc, y::AbstractArray)    y     for i in 1:length(y) ; dy[i] += ds1 ; end


text:
using Test

input = parse.(Int, readlines("01.input"))

star1(input) = sum(input)
@show star1(input)

function star2(input)
    seen = Set{Int}([0])
    freq = 0
    for i in Iterators.cycle(input)
        freq += i
        if freq ∈ seen
            return freq
        else
            push!(seen, freq)
        end
    end
end

@test [+1, -1] |> star2 == 0
@test [+3, +3, +4, -2, -4] |> star2 == 10
@test [-6, +3, +8, +5, -6] |> star2 == 5
@test [+7, +7, -2, -7, -4] |> star2 == 14

@show star2(input)

# can maybe use a BitSet to speed things up (represents a set of integers with a
# bitstring).
function star2bitset(input)
    seen = BitSet([0])
    freq = 0
    for i in Iterators.cycle(input)
        freq += i
        if freq ∈ seen
            return freq
        else
            push!(seen, freq)
        end
    end
end

@test [+1, -1] |> star2bitset == 0
@test [+3, +3, +4, -2, -4] |> star2bitset == 10
@test [-6, +3, +8, +5, -6] |> star2bitset == 5
@test [+7, +7, -2, -7, -4] |> star2bitset == 14

star2bitset(input)

using BenchmarkTools
@btime star2(input)             # 3.8 ms
@btime star2bitset(input)       # 301 μs


text:
push!(LOAD_PATH, "../src")
push!(LOAD_PATH, "../src/models")
push!(LOAD_PATH, "../src/models/c_functions")
ENV["MOCHA_USE_NATIVE_EXT"] = "true"
ENV["OMP_NUM_THREADS"] = 1
using Mocha
using SGMCMC
using MLUtilities

using MochaClassificationWrapper
using MochaClassification
using Base.Test
using GaussianMixture
using Banana
using MatrixFactorisation
using Logging
import StatsBase
Logging.configure(level=ERROR)
### test data models
@testset "Data models" begin
    include("test_banana.jl")
    include("test_gaussian_mixture.jl")
    include("test_logistic_regression.jl")
    include("test_mochaclassification.jl")
    include("test_mocharegression.jl")
    include("test_matrixfactorisation.jl")
end
### test optimizers
include("test_opt.jl")
### test samplers
@testset "Samplers" begin
    include("test_sgld.jl")
    include("test_sghmc.jl")
    include("test_sgnht.jl")
    include("test_relsghmc.jl")
    include("test_relsgnht.jl")
    include("test_msgnht.jl")
    include("test_sgldadagrad.jl")
    include("test_sgldrmsprop.jl")
    include("test_sgldadam.jl")
end


text:
using CubicEoS

mix = load(BrusilovskyEoSMixture;
    # names=("methane",),
    # names=("methane", "n-pentane"),
    # names=("nitrogen", "methane", "propane", "n-decane"),
    # names=("nitrogen","carbon_dioxide","hydrogen_sulfide","methane","ethane","propane","n-butane","n-pentane"),
    # names=("nitrogen","carbon_dioxide","hydrogen_sulfide","methane","ethane","propane","i-butane","n-butane","i-pentane","n-pentane","n-hexane","n-heptane","n-octane","n-nonane","n-decane","undecane","dodecane","tridecane","tetradecane","pentadecane","hexadecane","heptadecane","octadecane","nonadecane","icosane"),
    names=("carbon_dioxide", "n-decane"),
)

# χ = [1.0]
χ = [0.547413, 0.452587]
# χ = [0.2463, 0.2208, 0.2208, 0.3121]
# χ = [0.00325,0.01556,0.03329,0.82829,0.08552,0.01725,0.01015,0.00669]
# χ = [0.00325,0.01556,0.03329,0.82829,0.05850,0.02702,0.00543,0.00562,0.00248,0.00198,0.00174,0.00179,0.00275,0.00233,0.00191,0.00137,0.00112,0.00098,0.00087,0.00076,0.00053,0.00042,0.00039,0.00031,0.00131]

V = 1e-6
N = (4000 * V) .* χ  # 8500 phase transition c1c5, c1 4000
RT = CubicEoS.GAS_CONSTANT_SI * 250  # 170 c1, 371 c1c5

state = CubicEoS.vt_flash(mix, N, V, RT)
dump(state)

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
# Pressure plot alone a line across multiple frames.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, Glob, Printf
using Vlasiator: RE # Earth radius [m]

files = glob("bulk*.vlsv", ".")
nfile = length(files)

x1, x2 = 8.0, 29.0
point1 = [x1, 0, 0] .* RE
point2 = [x2, 0, 0] .* RE

cellids, distances, coords =
   load(files[1]) do meta
      getcellinline(meta, point1, point2)
   end

fig, ax = plt.subplots(figsize=(8, 4.8))

for (i, file) in enumerate(files)
   println("i = $i/$nfile, file = $file")
   local meta = load(file)

   p_extract = readvariable(meta, "vg_pressure", cellids) .* 1e9 |> vec # [nPa]
   rho_extract = readvariable(meta, "proton/vg_rho", cellids) |> vec
   v_extract = readvariable(meta, "proton/vg_v", cellids)
   vmag2_extract = sum(x -> x*x, v_extract, dims=1) |> vec
   pram_extract = rho_extract .* Vlasiator.mᵢ .* vmag2_extract .* 1e9 # [nPa]
   loc = range(x1, x2, length=length(rho_extract))

   ax.plot(loc, pram_extract, label="ram")
   ax.plot(loc, p_extract, label="thermal")

   ax.set_xlim(x1, x2)
   ax.set_ylim(0.0, 1.8)
   ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   grid(true)

   xlabel("x [Re]", fontsize=14)
   ylabel("Pressure [nPa]", fontsize=14)
   str_title = @sprintf "t= %4.1fs" meta.time
   title(str_title, fontsize=14)
   legend(loc="upper left", fontsize=14)
   savefig("out/"*meta.name[end-8:end-5]*".png", bbox_inches="tight")
   ax.cla()
end

close(fig)
println("plot finished!")