text:
module Bar
using ExportAll

function foo()
  1
end

function bar()
  2
end

function barBar()
  3
end

@exportAll()
end

ext: jl


text:
using Test
using StaticArrays
using RobotDynamics: DynamicsJacobian

n,m = 4,1
a = SizedMatrix{4,5}(rand(n,n+m))
F = DynamicsJacobian(a)
F[1,1] = 2
@test F[1,1] == 2
@test F[1,2] == a[1,2]
@test F.data === a

@test size(F.A) == (n,n)
@test size(F.B) == (n,m)
@test F.A == a[:,1:n]
@test F.B == a[:,n .+ (1:m)]

F.A[1,2] = 3
@test F[1,2] == 3
F.B[3,1] = 4
@test F[3,n+1] == 4

@test RobotDynamics.get_A(F) ≈ F.A
@test RobotDynamics.get_B(F) ≈ F.B
@test RobotDynamics.get_static_A(F) ≈ F.A
@test RobotDynamics.get_static_B(F) ≈ F.B
@test RobotDynamics.get_static_A(F) isa SMatrix{n,n}
@test RobotDynamics.get_static_B(F) isa SMatrix{n,m}

F_ = SMatrix(F)
@test DynamicsJacobian{n,n+m}(Tuple(F_)) isa DynamicsJacobian{n,n+m,Float64}
@test DynamicsJacobian{n,n+m,Float32}(Tuple(F_)) isa DynamicsJacobian{n,n+m,Float32}
@test Tuple(F) == Tuple(F_)
@test DynamicsJacobian(F_) isa DynamicsJacobian
@test DynamicsJacobian(F_) ≈ F

@test RobotDynamics.get_data(F_) === F_
@test RobotDynamics.get_data(F) === F.data.data
@test RobotDynamics.get_data(F.data) === F.data.data
@test RobotDynamics.get_data(F) isa Matrix
@test RobotDynamics.get_data(F.data) isa Matrix

############################################################################################
#                                    JACOBIAN-VECTOR PRODUCTS
############################################################################################
# model = Cartpole()
# λ = @SVector rand(state_dim(model))
# F = RD.DynamicsJacobian(model)
# x,u = rand(model)
# z = StaticKnotPoint(x,u,0.1)
# discrete_jacobian!(RK4, F, model, z)
# dg0 = F'λ
# jacobian!(F, model, z)
# g0 = F'λ

# RD.diffmethod(::Cartpole) = RD.ForwardAD()
# g = zeros(sum(size(model)))
# RD.discrete_jvp!(RK4, g, model, z, λ)
# @test g ≈ dg0
# @test (@allocated RD.discrete_jvp!(RK4, g, model, z, λ)) == 0

# RD.jvp!(g, model, z, λ)
# @test g ≈ g0
# @test (@allocated RD.jvp!(g, model, z, λ)) == 0

# # Finite difference
# RD.diffmethod(::Cartpole) = RD.FiniteDifference()
# cache = RD.gen_grad_cache(model)
# @test cache isa FiniteDiff.GradientCache
# g .= 0
# RD.discrete_jvp!(RK4, g, model, z, λ, cache)
# @test g ≈ dg0 atol=1e-6
# @test (@allocated RD.discrete_jvp!(RK4, g, model, z, λ, cache)) == 0

# RD.jvp!(g, model, z, λ, cache)
# @test g ≈ g0 atol=1e-6
# @test (@allocated RD.jvp!(g, model, z, λ, cache)) == 0

############################################################################################
#                               DYNAMICS HESSIANS
############################################################################################
# n,m = size(model)
# H = zeros(n+m,n+m)
# RD.∇jacobian!(H, model, z, λ)
# RD.∇discrete_jacobian!(RK4, H, model, z, λ)


text:
using Gen
using Test
import Random

"""
Compute a numerical partial derivative of `f` with respect to the `i`th
argument using finite differences.

If `broadcast` is `false`, then `args[i]` must be a scalar.

If `broadcast` is `true`, then `args[i]` may be an array.  In that case, `f` is
assumed to be the broadcast of a function whose `i`th argument is a scalar.  In
particular, `f` must still operate independently on each element of `args[i]`.
This condition cannot be checked automatically for arbitrary `f::Function`, so
the caller must guarantee it.
"""
function finite_diff(f::Function, args::Tuple, i::Int, dx::Float64;
                     broadcast=false)
    pos_args = Any[args...]
    neg_args = Any[args...]
    if broadcast
        pos_args[i] = copy(args[i]) .+ dx
        neg_args[i] = copy(args[i]) .- dx
        return (f(pos_args...) - f(neg_args...)) ./ (2. * dx)
    else
        pos_args[i] += dx
        neg_args[i] -= dx
        return (f(pos_args...) - f(neg_args...)) / (2. * dx)
    end
end

function finite_diff_vec(f::Function, args::Tuple, i::Int, j::Int, dx::Float64)
    pos_args = Any[deepcopy(args)...]
    pos_args[i][j] += dx
    neg_args = Any[deepcopy(args)...]
    neg_args[i][j] -= dx
    return (f(pos_args...) - f(neg_args...)) / (2. * dx)
end

"""
Compute a numerical partial derivative of `f` with respect to the `i`th argument, a symettric matrix, using finite differences.

Finite differences are applied symettrically to off-diagonal elements of the matrix, ensuring that the modified matrix is still symettric.

When applied to off-diagonal elements, the finite difference is scaled by a factor of 1/2.

The condition that the `i`th argument is a symettric matrix cannot be checked automatically, so the caller must guarantee it.
"""

function finite_diff_mat_sym(f::Function, args::Tuple, i::Int, j::Int, k::Int, dx::Float64)
    pos_args = Any[deepcopy(args)...]
    pos_args[i][j, k] += dx
    neg_args = Any[deepcopy(args)...]
    neg_args[i][j, k] -= dx

    if j!=k
        pos_args[i][k, j] += dx
        neg_args[i][k, j] -= dx
        return (f(pos_args...) - f(neg_args...)) / (4. * dx)
    end

    return (f(pos_args...) - f(neg_args...)) / (2. * dx)
end

const dx = 1e-6

include("autodiff.jl")
include("diff.jl")
include("selection.jl")
include("assignment.jl")
include("dynamic_dsl.jl")
include("optional_args.jl")
include("static_ir/static_ir.jl")
include("static_dsl.jl")
include("tilde_sugar.jl")
include("inference/inference.jl")
include("modeling_library/modeling_library.jl")


text:
abstract type Operator end
struct Sum <: Operator end
struct Prod <: Operator end
Base.isequal(a::Operator, b::Operator) = (typeof(a) == typeof(b))
Base.:(==)(a::Operator, b::Operator) = Base.isequal(a, b)
apply(o::Operator, diags) = error("not implemented!")

Base.show(io::IO, o::Operator) = print(io, typeof(o))
Base.show(io::IO, o::Sum) = print(io, "⨁")
Base.show(io::IO, o::Prod) = print(io, "Ⓧ")

"""
    abstract type DiagramId end

    The abstract type of all diagrams/subdiagrams/bare propagators
"""
abstract type DiagramId end

"""
    abstract type PropagatorId <: DiagramId end

    The abstract type of all bare propagators
"""
abstract type PropagatorId <: DiagramId end

# Base.Dict(x::DiagramId) = Dict{Symbol,Any}([fn => getfield(x, fn) for fn ∈ fieldnames(typeof(x))])
# Base.show(io::IO, d::DiagramId) = error("Base.show not implemented!")
# Base.isequal(a::DiagramId, b::DiagramId) = error("Base.isequal not implemented!")
Base.:(==)(a::DiagramId, b::DiagramId) = Base.isequal(a, b)

struct BareGreenId <: PropagatorId
    para::GenericPara
    type::AnalyticProperty #Instant, Dynamic, D_Instant, D_Dynamic
    extK::Vector{Float64}
    extT::Tuple{Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function BareGreenId(para::GenericPara, type::AnalyticProperty = Dynamic, order = [0, 0]; k, t)
        return new(para, type, k, Tuple(t), order)
    end
end
Base.show(io::IO, v::BareGreenId) = print(io, "$(short(v.type))#$(v.order), k$(v.extK), t$(v.extT)")

struct BareInteractionId <: PropagatorId
    para::GenericPara
    response::Response #UpUp, UpDown, ...
    type::AnalyticProperty #Instant, Dynamic, D_Instant, D_Dynamic
    permutation::Permutation
    extK::Vector{Float64}
    extT::Tuple{Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function BareInteractionId(para::GenericPara, response::Response, type::AnalyticProperty = Instant, order = [0, 0]; k, t = (0, 0), permu::Permutation = DiEx)
        return new(para, response, type, permu, k, Tuple(t), order)
    end
end
Base.show(io::IO, v::BareInteractionId) = print(io, "$(short(v.response))$(short(v.type))$(v.permutation)#$(v.order), k$(v.extK), t$(v.extT)")

struct GenericId <: DiagramId
    para::GenericPara
    extra::Any
    order::Vector{Int}
    GenericId(para::GenericPara, extra = Nothing, order = [0, 0]) = new(para, extra, order)
end
Base.show(io::IO, v::GenericId) = print(io, v.extra == Nothing ? "#$(v.order)" : "$(v.extra)#$(v.order)")

struct GreenId <: DiagramId
    para::GenericPara
    type::AnalyticProperty #Instant, Dynamic, D_Instant, D_Dynamic
    extK::Vector{Float64}
    extT::Tuple{Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function GreenId(para::GenericPara, type::AnalyticProperty = Dynamic, order = [0, 0]; k, t)
        return new(para, type, k, Tuple(t), order)
    end
end
Base.show(io::IO, v::GreenId) = print(io, "$(short(v.type))#$(v.order), k$(v.extK), t$(v.extT)")

struct SigmaId <: DiagramId
    para::GenericPara
    type::AnalyticProperty #Instant, Dynamic, D_Instant, D_Dynamic
    extK::Vector{Float64}
    extT::Tuple{Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function SigmaId(para::GenericPara, type::AnalyticProperty, order = [0, 0]; k, t = (0, 0))
        return new(para, type, k, t, order)
    end
end
Base.show(io::IO, v::SigmaId) = print(io, "$(short(v.type))#$(v.order), t$(v.extT)")

struct PolarId <: DiagramId
    para::GenericPara
    response::Response #UpUp, UpDown, ...
    extK::Vector{Float64}
    extT::Tuple{Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function PolarId(para::GenericPara, response::Response, order = [0, 0]; k, t = (0, 0))
        return new(para, response, k, t, order)
    end
end
Base.show(io::IO, v::PolarId) = print(io, "$(short(v.response))#$(v.order), k$(v.extK), t$(v.extT)")

struct Ver3Id <: DiagramId
    para::GenericPara
    response::Response #UpUp, UpDown, ...
    extK::Vector{Vector{Float64}}
    extT::Tuple{Int,Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function Ver3Id(para::GenericPara, response::Response, order = [0, 0]; k, t = (0, 0, 0))
        return new(para, response, k, Tuple(t), order)
    end
end
Base.show(io::IO, v::Ver3Id) = print(io, "$(short(v.response))#$(v.order),t$(v.extT)")

struct Ver4Id <: DiagramId
    para::GenericPara
    response::Response #UpUp, UpDown, ...
    type::AnalyticProperty #Instant, Dynamic, D_Instant, D_Dynamic
    channel::TwoBodyChannel # particle-hole, particle-hole exchange, particle-particle, irreducible
    extK::Vector{Vector{Float64}}
    extT::Tuple{Int,Int,Int,Int} #all possible extT from different interactionType
    order::Vector{Int}
    function Ver4Id(para::GenericPara, response::Response, type::AnalyticProperty = Dynamic, order = [0, 0]; k, t = (0, 0, 0, 0), chan::TwoBodyChannel = AnyChan)
        return new(para, response, type, chan, k, Tuple(t), order)
    end
end
Base.show(io::IO, v::Ver4Id) = print(io, (v.channel == AnyChan ? "" : "$(v.channel) ") * "$(short(v.response))$(short(v.type))#$(v.order),t$(v.extT)")

function vstr(r, c)
    N = length(r)
    # cstr(x) = x ? "⁺" : "⁻"
    s = ""
    for i = 1:N-1
        s *= "$(r[i])$c"
    end
    s *= "$(r[end])$c"
    return s
end

function vcstr(r, creation)
    N = length(r)
    # cstr(x) = x ? "⁺" : "⁻"
    s = ""
    for i = 1:N-1
        if creation[i]
            s *= "$(r[i])⁺"
        else
            s *= "$(r[i])⁻"
        end
    end
    if creation[end]
        s *= "$(r[end])⁺"
    else
        s *= "$(r[end])⁻"
    end
    return s
end

"""
hopping function c⁺c⁻
"""
struct BareHoppingId <: PropagatorId
    para::GenericPara
    site::Tuple{Int,Int}
    orbital::Tuple{Int,Int}
    extT::Tuple{Int,Int}
    function BareHoppingId(para::GenericPara, orbital, t, r)
        return new(para, r, orbital, t)
    end
end
Base.show(io::IO, v::BareHoppingId) = print(io, "($(vstr(v.site, "ᵣ"))|$(vstr(v.orbital, "ₒ"))|$(vcstr(v.extT, [true, false])))")

"""
time-ordered N-point Bare Green's function
"""
struct BareGreenNId <: PropagatorId
    para::GenericPara
    site::Int
    creation::Vector{Bool}
    orbital::Vector{Int}
    extT::Vector{Int}
    N::Int
    function BareGreenNId(para::GenericPara; orbital = [], t = [], creation = [], r = 0)
        @assert length(orbital) == length(t) == length(creation)
        return new(para, r, creation, orbital, t, length(orbital))
    end
end
Base.show(io::IO, v::BareGreenNId) = print(io, "($(v.site)ᵣ|$(vstr(v.orbital, "ₒ"))|$(vcstr(v.extT, v.creation)))")

"""
time-ordered N-point Composite Green's function
"""
struct GreenNId <: DiagramId
    para::GenericPara
    site::Vector{Int}
    creation::Vector{Bool}
    orbital::Vector{Int}
    extT::Vector{Int}
    N::Int
    function GreenNId(para::GenericPara; orbital = [], t = [], creation = [], r = [])
        @assert length(orbital) == length(t) == length(r) == length(creation)
        return new(para, r, creation, orbital, t, length(orbital))
    end
end
Base.show(io::IO, v::GreenNId) = print(io, "($(vstr(v.site, "ᵣ"))|$(vstr(v.orbital, "ₒ"))|$(vcstr(v.extT, v.creation)))")

"""
time-ordered N-point Composite Green's function
"""
struct ConnectedGreenNId <: DiagramId
    para::GenericPara
    site::Vector{Int}
    creation::Vector{Bool}
    orbital::Vector{Int}
    extT::Vector{Int}
    N::Int
    function ConnectedGreenNId(para::GenericPara; orbital = [], t = [], creation = [], r = [])
        @assert length(orbital) == length(t) == length(r) == length(creation)
        return new(para, r, creation, orbital, t, length(orbital))
    end
end
Base.show(io::IO, v::ConnectedGreenNId) = print(io, "($(vstr(v.site, "ᵣ"))|$(vstr(v.orbital, "ₒ"))|$(vcstr(v.extT, v.creation)))")

function Base.isequal(a::DiagramId, b::DiagramId)
    if typeof(a) != typeof(b)
        return false
    end
    for field in fieldnames(typeof(a))
        if field == :extK
            if (getproperty(a, :extK) ≈ getproperty(b, :extK)) == false
                return false
            end
        end
        if getproperty(a, field) != getproperty(b, field)
            return false
        end
    end
    return true
end


text:
using BHAtp
using Test

# write your own tests here

tests = Array{String, 1}([
  "case_01"
])

for test in tests
  println()
  println("\nRunning $(test)_test.jl.\n")
  include(test*"_test.jl")
  println()
end


text:
using Random, JuMP, JuDGE, Test, DelimitedFiles

include("solvers/setup_gurobi.jl")
#include("solvers/setup_cplex.jl")
#include("solvers/setup_coin.jl")
#include("solvers/setup_glpk.jl")

function transportation()
   mytree = narytree(5,2)

   function invest_supply_cost(node)
      if node.parent==nothing
         return Dict( zip( supply_nodes, [1.0,2.0]) )
      else
         p=node.parent
         for i in 1:length(p.children)
            if p.children[i]==node
               temp=deepcopy(invest_supply_cost(p))
               for key in keys(temp)
                  temp[key]*=(0.2*i+0.8)
               end
               return temp
            end
         end
      end
   end

   function invest_arc_cost(node)
      if node.parent==nothing
         temp=[]
         for i in supply_nodes
            for j in demand_nodes
               push!(temp,(i,j))
            end
         end
         return Dict( zip( temp, [1.0,2.0,3.0,4.0,5.0,6.0]) )
      else
         p=node.parent
         for i in 1:length(p.children)
            if p.children[i]==node
               temp=deepcopy(invest_arc_cost(p))
               for key in keys(temp)
                  temp[key]*=(0.2*i+0.8)
               end
               return temp
            end
         end
      end
   end

   function demand(node)
      if node.parent==nothing
         return d_dict
      else
         p=node.parent
         for i in 1:length(p.children)
            if p.children[i]==node
               temp=deepcopy(demand(p))
               for key in keys(temp)
                  temp[key]+=i*i
               end
               return temp
            end
         end
      end
   end

   function supply(node)
      return s_dict
   end

   data_file = "transportation.csv"
   data = readdlm(joinpath(@__DIR__,data_file),',')

   supply_nodes = data[3:end, 2]
   s = data[3:end, 1]

   demand_nodes = collect(data[2, 3:end])
   d = collect(data[1, 3:end])

   c = data[3:end, 3:end]

   # Converting arrays to dictionaries
   s_dict = Dict( zip( supply_nodes, s*2) )
   d_dict = Dict( zip( demand_nodes, d) )

   c_dict = Dict()
   for i in 1:length(supply_nodes)
     for j in 1:length(demand_nodes)
       c_dict[supply_nodes[i], demand_nodes[j]] = c[i,j]
     end
   end

   ### with judge
   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)

      @expansion(model, new_supply[supply_nodes], Bin) #invest in more supply
      @expansion(model, 0<=new_capacity[supply_nodes,demand_nodes]<=2, Int) #invest in more arc capacity

      @capitalcosts(model, sum(invest_supply_cost(node)[i] * new_supply[i] for i in supply_nodes) +
         sum(invest_arc_cost(node)[i,j] * new_capacity[i,j] for i in supply_nodes for j in demand_nodes))

      @variable(model, x[supply_nodes, demand_nodes] >= 0)

      @objective(model, Min, sum(c_dict[i,j]*x[i,j] for i in supply_nodes, j in demand_nodes))
      @constraint(model, SupplyIncrease[i in supply_nodes],
                  sum(x[i,j] for j in demand_nodes) <= supply(node)[i] + s_dict[i]*new_supply[i])

      @constraint(model, CapacityIncrease[i in supply_nodes, j in demand_nodes],
                  x[i,j] <= c_dict[i,j] + c_dict[i,j]*new_capacity[i,j])

      @constraint(model,demandCon[j in demand_nodes], sum(x[i,j] for i in supply_nodes) == demand(node)[j] )

      return model
   end

   function format_output(s::Symbol,values)
      if s==:new_capacity
         output=Dict{Tuple,Float64}()
         for i in supply_nodes
            for j in demand_nodes
               output[i,j]=values[i,j]*c_dict[i,j]
            end
         end
         return output
      elseif s==:new_supply
         output=Dict{String,Float64}()
         for i in supply_nodes
            output[i]=values[i]*s_dict[i]
         end
         return output
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver, discount_factor=0.9)
   JuDGE.solve(judy,termination=Termination(inttol=10^-7))

   println("\nObjective: "*string(objective_value(judy.master_problem))*"\n")
   JuDGE.print_expansions(judy,format=format_output)

   println("\nRe-solved Objective: " * string(resolve_subproblems(judy)))
   solution=JuDGE.solution_to_dictionary(judy)

   JuDGE.visualize_tree(mytree,solution)
   JuDGE.write_solution_to_file(judy,joinpath(@__DIR__,"transport_solution_decomp.csv"))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver, discount_factor=0.9)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(objective_value(deteq.problem)))
   JuDGE.write_solution_to_file(deteq,joinpath(@__DIR__,"transport_solution_deteq.csv"))
   return objective_value(judy.master_problem)
end

@test transportation() ≈ 1924.35 atol = 1e-2


text:
using Catalyst, DelaySSAToolkit
# 0 -->A1, => A2, ..., => AN => 0 with delay 1 for each cascade delay, N being the length of the delay reaction chain.

rn = @reaction_network begin
    5, 0 --> A1
end

delay_trigger_affect! = []
chain_len = 10
delay_trigger_affect! = function (integrator, rng)
    append!(integrator.de_chan[1], 1.)
end

u0 = zeros(chain_len)
de_chan0 = []
for _ in 1:chain_len
    push!(de_chan0, [])
end
tspan = (0.,50.)
delay_complete_affect! = []
for i in 1:chain_len-1
    push!(delay_complete_affect!, function (integrator, rng)
        integrator.u[i] -= 1 # A_prev minus 1
        integrator.u[i+1] += 1 # A plus 1
        append!(integrator.de_chan[i+1], 1.) # add to the delay channel
    end
    )
end
push!(delay_complete_affect!, function (integrator, rng)
    integrator.u[chain_len] -= 1 # A_prev minus 1
end)

delay_trigger = Dict(Pair(1, delay_trigger_affect!))
delay_complete = Dict(i=>delay_complete_affect![i] for i in 1:chain_len)
delay_interrupt = Dict()

delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)
dprob = DiscreteProblem(jumpsys, u0, tspan)

using Test
algos = [DelayRejection(),DelayDirect(), DelayMNRM(),  DelayDirectCR()]
@testset for algo in algos
    djprob = DelayJumpProblem(jumpsys, dprob, algo, delayjumpset, de_chan0, save_positions=(false,false), save_delay_channel= true)
    @info "Testing method $(algo)"
    @time sol = solve(djprob, SSAStepper(), saveat = 1.)
    for i in 1:length(sol.t)-1, j in 1:chain_len-1
        err = sum(abs2, sol.channel[i][j] .- sol.channel[i+1][j+1])
        @test err < 1e-20
    end
end


text:
include("5.5.1 u (a+b arcsec(c x))^n.jl")
include("5.5.2 Miscellaneous inverse secant.jl")

max_stars_repo_path: src/rules/5 Inverse trig functions/5.5 Inverse secant/.jl


text:
# loaders.jl
# utilities for loading in data / images

"""
`data = ir_load_brainweb_t1_256()`

load brainweb T1-weighted MRI slice of size 256x256
"""
function ir_load_brainweb_t1_256()
	file = joinpath(dirname(pathof(MIRT)), "../data/mri/brainweb_t1_256.fld")
	return Float32.(fld_read(file))' # transpose!
end


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

using Base.Test
using RBT: LiteralSyntax, ApplySyntax, ComposeSyntax, syntax

q1 = syntax("department")

@test string(q1) == "department"
@test typeof(q1) == ApplySyntax
@test q1.fn == :department
@test q1.args == []

q2 = syntax("department.name")

@test string(q2) == "department.name"
@test typeof(q2) == ComposeSyntax
@test string(q2.f) == "department"
@test string(q2.g) == "name"

q3 = syntax("""
    department
    :filter(count(employee)>100)
    :sort(count(employee):desc)
    :select(
        name,
        count(employee),
        employee)
""")

@test typeof(q3) == ApplySyntax
@test q3.fn == :select
@test string(q3.args[1]) == "department:filter(>(count(employee),100)):sort(count(employee):desc)"
@test string(q3.args[2]) == "name"
@test string(q3.args[3]) == "count(employee)"
@test string(q3.args[4]) == "employee"

q4 = syntax("""
    employee
    :group(position)
    :select(
        position,
        count(employee),
        max(employee.salary),
        mean(employee.salary))
""")

@test typeof(q4) == ApplySyntax
@test q4.fn == :select
@test string(q4.args[1]) == "employee:group(position)"
@test string(q4.args[2]) == "position"
@test string(q4.args[3]) == "count(employee)"
@test string(q4.args[4]) == "max(employee.salary)"
@test string(q4.args[5]) == "mean(employee.salary)"

q5 = syntax("""
    department:filter(count(employee)>100).employee.name
""")

@test typeof(q5) == ComposeSyntax
@test string(q5.f) == "department:filter(>(count(employee),100))"
@test string(q5.g) == "employee.name"

q6 = syntax("6 * (3 + 4)")

@test string(q6) == "*(6,+(3,4))"
@test typeof(q6) == ApplySyntax
@test typeof(q6.args[1]) == LiteralSyntax
@test q6.args[1].val == 6

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# My formulation

#=
@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1.0
          den = 1.0
          For( j=1, j <= N, Increment(j),
            begin
              num = If(j != i, num = num * (x-xi[j]), num)
              den = If(j != i, den = den * (xi[i]-xi[j]), den)
            end
          )
          sum = sum + yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end
=#

# John's better formulation, much better, need to update tests

@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1
          den = 1
          For( j=1, j <= N, Increment(j),
            begin
              If(j != i, num *= (x-xi[j]))
              If(j != i, den *= (xi[i]-xi[j]))
            end
          )
          sum += yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsSimulator API
# -----------------------------------------------------------------------------

"""
    NBSimulator <: MolecularDynamicsSimulator

A wrapper around NBodySimulator to implement the Atomistic API.

**Field descriptions**
- `Δt::Real` the time between timesteps, assumed to be in atomic units
- `steps::Integer` the number of timesteps for the simulation
- `t₀::Real` the starting time of the simulation, assumed to be in atomic units;
    defaults to 0
- `thermostat::Thermostat` the thermostat for the simulation;
    many options are defined by `NBodySimulator`, but a user could also define a custom thermostat;
    defaults to the `NullThermostat`
- `simulator::OrdinaryDiffEqAlgorithm` the algorithm to be used for the ODE;
    defaults to VelocityVerlet
- `potentials::Dict{Symbol,PotentialParameters}` dictionary of potentials;
    shouldn't be manipulated directly by the user
"""
@kwdef struct NBSimulator <: MolecularDynamicsSimulator
    Δt::Real       # in TIME_UNIT
    steps::Integer
    t₀::Real = 0.0 # in TIME_UNIT
    thermostat::Thermostat = NullThermostat()
    simulator::OrdinaryDiffEqAlgorithm = VelocityVerlet()
    potentials::Dict{Symbol,PotentialParameters} = Dict{Symbol,PotentialParameters}()
end
function NBSimulator(Δt::Real, steps::Integer; kwargs...)
    NBSimulator(; Δt = Δt, steps = steps, kwargs...)
end
function NBSimulator(Δt::Unitful.Time, steps::Integer; t₀::Unitful.Time = 0.0 * TIME_UNIT, kwargs...)
    NBSimulator(; Δt = austrip(Δt), steps = steps, t₀ = austrip(t₀), kwargs...)
end

# Extract the tuple of start_time, end_time from the simulator
time_range(simulator::NBSimulator) = (simulator.t₀, simulator.t₀ + simulator.steps * simulator.Δt)

function simulate(system::AbstractSystem, simulator::NBSimulator, potential::ArbitraryPotential)
    simulator.potentials[:custom] = InteratomicPotentialParameters(potential)
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator, potential::LennardJonesParameters)
    simulator.potentials[:lennard_jones] = potential
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator)
    nb_system = PotentialNBodySystem{ElementMassBody}(bodies(system), simulator.potentials)
    simulation = NBodySimulation(nb_system, time_range(simulator), nbs_boundary_conditions(system), simulator.thermostat, austrip(u"k"))
    NBSResult(run_simulation(simulation, simulator.simulator, dt = simulator.Δt))
end

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

# Internal struct that represents a potential in the format accepted by NBodySimulator
# Wraps the underlying InteratomicPotential with a cache of the forces for the current timestep
struct InteratomicPotentialParameters <: PotentialParameters
    potential::ArbitraryPotential
    timestep_cache::Ref{Real}
    force_cache::Ref{Vector{SVector{3,Real}}}
    InteratomicPotentialParameters(potential::ArbitraryPotential) = new(potential, Ref{Real}(), Ref{Vector{SVector{3,Real}}}())
end

function NBodySimulator.get_accelerating_function(parameters::InteratomicPotentialParameters, simulation::NBodySimulation)
    bodies = simulation.system.bodies
    boundary_conditions = simulation.boundary_conditions
    (dv, u, v, t, i) -> begin
        if !isassigned(parameters.timestep_cache) || t != parameters.timestep_cache[]
            particles = [ElementMassBody(bodies[i], SVector{3}(u[:, j]), SVector{3}(v[:, j])) for j ∈ 1:length(bodies)]
            system = DynamicSystem(FlexibleSystem(particles, boundary_conditions), t * TIME_UNIT)
            parameters.timestep_cache[] = t
            parameters.force_cache[] = force(system, parameters.potential)
        end
        dv .+= parameters.force_cache[][i] / bodies[i].m
    end
end

# -----------------------------------------------------------------------------
# Unitful Support for Built-in LennardJonesParameters
# -----------------------------------------------------------------------------

function NBodySimulator.LennardJonesParameters(ϵ::Unitful.Energy, σ::Unitful.Length, R::Unitful.Length)
    LennardJonesParameters(austrip(ϵ), austrip(σ), austrip(R))
end


text:
using MultivariateTimeSeries
using Random,Test

let
    Random.seed!(0)
    R = rand(20,3)
    d = DataFrame(R)
    mts = MTS(d, [1,6,11,16])
    labels = rand(1:3, 4)

    @test nrow(mts.data) == 20
    @test mts.index == [1,6,11,16]
    @test eltype(mts.views) <: AbstractDataFrame
    @test length(mts) == 4

    @test size(mts[1]) == (5,3)
    @test size(mts) == (4,)
    @test ncol(mts) == 3
    @test names(mts) == ["x1", "x2", "x3"]
    @test length(collect(mts)) == length(mts.index)

    file = joinpath(dirname(@__FILE__), "mts.zip")
    write_data(file, mts, labels)

    @test has_labels(file)
    mts2, labels2 = read_data_labeled(file)
    mts3 = read_data(file)

    @test mts == mts2 == mts3
    @test labels == labels2
    v = vec(mts)

    mts4 = normalize01(mts)
    @test maximum(maximum(mts4)) <= 1.0
    @test minimum(minimum(mts4)) >= 0.0
end

let
    Random.seed!(0)
    v = [DataFrame(rand(rand(3:5), 3)) for i=1:5]
    mts = MTS(v)
    @test all(x==y for (x,y) in zip(mts, v))
end

let
    Random.seed!(0)
    N = 5
    v1 = [DataFrame(rand(rand(3:5), 3)) for i=1:N]
    v2 = [DataFrame(rand(rand(3:5), 3)) for i=1:N]
    v3 = [DataFrame(rand(rand(3:5), 3)) for i=1:N]
    m1 = MTS(v1)
    m2 = MTS(v2)
    m3 = MTS(v3)
    mts = vcat(m1, m2, m3)

    @test all([mts[i] == m1[i] for i=1:N])
    @test all([mts[i+N] == m2[i] for i=1:N])
    @test all([mts[i+2N] == m3[i] for i=1:N])
end

let
    Random.seed!(0)
    A = [1 2; 3 4; 5 6; 7 8]
    d = DataFrame(A)
    mts = MTS(d, [1, 3])
    inds = findfirst(d->d[:x1] > 2, mts)
    @test inds == [2, 1]
    mts1 = trim(mts, inds)
    @test map(nrow, mts1) == [2, 1]
    mts1 = transform(d->d[:,[:x1]], mts)
    @test map(ncol, mts1) == [1, 1]
end

let
    mts = test_mts()
    mts1 = zero(mts)
    @test mts1[1][1,1] == 0.0
end


text:
function JSONRPC.parse_params(::Type{Val{Symbol("\$/cancelRequest")}}, params)
    return CancelParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("\$/cancelRequest")},CancelParams}, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("\$/setTraceNotification")}}, params)
    return params
end

function process(r::JSONRPC.Request{Val{Symbol("\$/setTraceNotification")},Dict{String,Any}}, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/lint-package")}}, params)
    return
end

function process(r::JSONRPC.Request{Val{Symbol("julia/lint-package")},Nothing}, server)
    for (uri, f) in server.documents
    end
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/toggle-lint")}}, params)
    return TextDocumentIdentifier(params["textDocument"])
end

function process(r::JSONRPC.Request{Val{Symbol("julia/toggle-lint")},TextDocumentIdentifier}, server)
    doc = server.documents[URI2(r.uri)]
    doc._runlinter = !doc._runlinter
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/reload-modules")}}, params)
end

function process(r::JSONRPC.Request{Val{Symbol("julia/reload-modules")},Nothing}, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/toggleFileLint")}}, params)
    return params
end

function process(r::JSONRPC.Request{Val{Symbol("julia/toggleFileLint")}}, server)
    # path = get(r.params, "path", "")
    # uri = get(r.params, "external", "")
    # if isdir(uri2filepath(path))
    #     for doc in values(server.documents)
    #         uri2 = doc._uri
    #         server.debug_mode && @info "LINT: ignoring $path"
    #         if startswith(uri2, uri)
    #             toggle_file_lint(doc, server)
    #         end
    #     end
    # else
    #     if uri in map(i->i._uri, values(server.documents))
    #         server.debug_mode && @info "LINT: ignoring $path"
    #         doc = server.documents[URI2(uri)]
    #         toggle_file_lint(doc, server)
    #     end
    # end
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/toggle-log")}}, params)
end

function process(r::JSONRPC.Request{Val{Symbol("julia/toggle-log")},Nothing}, server)
    server.debug_mode = !server.debug_mode
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/getCurrentBlockOffsetRange")}}, params)
    return TextDocumentPositionParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("julia/getCurrentBlockOffsetRange")}}, server)
    if !haskey(server.documents, URI2(r.params.textDocument.uri))
        send(JSONRPC.Response(r.id, CancelParams(r.id)), server)
        return
    end
    tdpp = r.params
    doc = server.documents[URI2(tdpp.textDocument.uri)]
    offset = get_offset(doc, tdpp.position)
    x = getcst(doc)
    loc = 0
    p1, p2, p3 = 1, x.span, x.fullspan
    if x.typ === CSTParser.FileH
        (offset > x.fullspan || x.args === nothing) && return 1, x.span, x.fullspan
        for a in x.args
            if loc <= offset < loc + a.fullspan
                if a.typ === CSTParser.ModuleH
                    if loc + a.args[1].fullspan + a.args[2].fullspan < offset < loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].fullspan
                        loc0 = loc +  a.args[1].fullspan + a.args[2].fullspan
                        loc += a.args[1].fullspan + a.args[2].fullspan
                        for b in a.args[3].args
                            if loc <= offset < loc + b.fullspan
                                p1, p2, p3 = loc + 1, loc + b.span, loc + b.fullspan
                                break
                            end
                            loc += b.fullspan
                        end
                    else
                        p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                    end
                elseif a.typ === CSTParser.TopLevel
                    p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                    for b in a.args
                        if loc <= offset < loc + b.fullspan
                            p1, p2, p3 = loc + 1, loc + b.span, loc + b.fullspan
                        end
                        loc += b.fullspan
                    end
                else
                    p1, p2, p3 = loc + 1, loc + a.span, loc + a.fullspan
                end
            end
            loc += a.fullspan
        end
    end

    response = JSONRPC.Response(r.id, (isempty(doc._content) ? 0 : length(doc._content, 1, max(1, p1)), length(doc._content, 1, p2), length(doc._content, 1, p3)))

    send(response, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/activateenvironment")}}, params)
    return params
end

function process(r::JSONRPC.Request{Val{Symbol("julia/activateenvironment")}}, server)
    server.env_path = r.params
    server.symbol_server = StaticLint.SymbolServer.SymbolServerProcess(depot = server.depot_path, environment=server.env_path)
    @info "Restarted symbol server"
    server.packages = StaticLint.SymbolServer.getstore(server.symbol_server)
    @info "StaticLint store set"
    kill(server.symbol_server)

    for (uri, doc) in server.documents
        parse_all(doc, server)
    end
    @info "Finished reparsing everything"
end


text:
"""
    AbstractDiffBackend

An abstract type for diff backends. See [`FiniteDifferencesBackend`](@ref) for
an example.
"""
abstract type AbstractDiffBackend end

struct NoneDiffBackend <: AbstractDiffBackend end

"""
    _derivative(f, t[, backend::AbstractDiffBackend])

Compute the derivative of a callable `f` at time `t` computed using the given `backend`,
an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean derivatives, for Riemannian differentiation see
for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _derivative end

_derivative(f, t) = _derivative(f, t, default_differential_backend())

function _derivative!(f, X, t, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _derivative(f, t, backend))
end

"""
    _gradient(f, p[, backend::AbstractDiffBackend])

Compute the gradient of a callable `f` at point `p` computed using the given `backend`,
an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _gradient end

_gradient(f, p) = _gradient(f, p, default_differential_backend())

function _gradient!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _gradient(f, p, backend))
end

"""
    _jacobian(f, p[, backend::AbstractDiffBackend])

Compute the jacobian of a callable `f` at point `p` computed using the given `backend`,
an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _jacobian end

_jacobian(f, p) = _jacobian(f, p, default_differential_backend())

function _jacobian!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _jacobian(f, p, backend))
end

"""
    CurrentDiffBackend(backend::AbstractDiffBackend)

A mutable struct for storing the current differentiation backend in a global
constant [`_current_default_differential_backend`](@ref).

# See also

[`AbstractDiffBackend`](@ref), [`default_differential_backend`](@ref), [`set_default_differential_backend!`](@ref)
"""
mutable struct CurrentDiffBackend
    backend::AbstractDiffBackend
end

"""
    _current_default_differential_backend

The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default
differentiation backend.
"""
const _current_default_differential_backend = CurrentDiffBackend(NoneDiffBackend())
"""
    default_differential_backend() -> AbstractDiffBackend

Get the default differentiation backend.
"""
default_differential_backend() = _current_default_differential_backend.backend

"""
    set_default_differential_backend!(backend::AbstractDiffBackend)

Set current backend for differentiation to `backend`.
"""
function set_default_differential_backend!(backend::AbstractDiffBackend)
    _current_default_differential_backend.backend = backend
    return backend
end


text:
function is_prime(n::Int64)
    if n < 4
        return true
    end
    if n % 2 == 0
        return false
    end
    for a = 3:round(sqrt(n))
        if n % a == 0
            return false
        end
    end
    return true
end

function calculate_primes(n::Int64)
    primes = Int64[]
    p = 1
    for i = 1:n
        while ! is_prime(p)
            p += 1
        end
        push!(primes, p)
        p += 1
    end
    return primes
end

calculate_primes(1000000)


text:
using Test
using DECAES
using Logging

# Arbitrary default required parameters used during testing
const default_paramdict = Dict{Symbol,Any}(
    :TE => 8e-3,
    :nT2 => 42,
    :T2Range => (12e-3, 1.8),
    :SPWin => (12e-3, 37e-3),
    :MPWin => (37e-3, 650e-3),
    :Reg => "lcurve",
)

# Legacy settings which previously had defaults
const legacy_default_paramdict = Dict{Symbol,Any}(
    :TE => 10e-3,
    :nT2 => 40,
    :T2Range => (15e-3, 2.0),
    :SPWin => (14e-3, 40e-3),
    :MPWin => (40e-3, 200e-3),
    :Reg => "chi2",
    :Chi2Factor => 1.02,
)

# Write 4D image to disk
function write_image(filename, image)
    if endswith(filename, ".mat")
        DECAES.MAT.matwrite(filename, Dict("img" => image))
    else
        DECAES.NIfTI.niwrite(filename, DECAES.NIfTI.NIVolume(image))
    end
end

# Call main function on image file `image`
function run_main(image, args; make_settings_file::Bool)

    # Write input image to file for reading
    inputfilename = args[1]
    outputpath = args[3]
    inputfilebasename = joinpath(outputpath, "input")
    write_image(inputfilename, image)

    # Run main, possibly writing CLI args to settings file first
    if make_settings_file
        settings_file = joinpath(outputpath, "settings.txt")
        open(settings_file, "w") do file
            println(file, join(args, "\n"))
        end
        main(["@" * settings_file])
    else
        main(args)
    end

    # Check that only requested files were created
    t2maps_file, t2dist_file, t2parts_file, settings_file = inputfilebasename .* (".t2maps.mat", ".t2dist.mat", ".t2parts.mat", ".settings.txt")
    T2map, T2part = ("--T2map" ∈ args), ("--T2part" ∈ args)

    @test !xor(T2map,  isfile(t2maps_file))
    @test !xor(T2map,  isfile(t2dist_file))
    @test !xor(T2part, isfile(t2parts_file))
    @test !xor(make_settings_file, isfile(settings_file))

    t2maps  = T2map  ? DECAES.MAT.matread(t2maps_file) : nothing
    t2dist  = T2map  ? DECAES.MAT.matread(t2dist_file)["dist"] : nothing
    t2parts = T2part ? DECAES.MAT.matread(t2parts_file) : nothing

    return (; t2maps, t2dist, t2parts)
end

function construct_args(paramdict;
        argstype,
        inputfilename = nothing,
        outputpath = nothing,
        quiet::Bool = true,
        legacy::Bool = false,
        T2map::Bool = true,
        T2part::Bool = true,
    )

    if legacy
        # Add legacy default options to paramdict
        for (k,v) in legacy_default_paramdict
            if k ∉ keys(paramdict)
                paramdict[k] = v
            end
        end
    end

    if argstype === :cli
        #### CLI

        args = [inputfilename, rand(["--output", "-o"]), outputpath]
        T2map  && push!(args, "--T2map")
        T2part && push!(args, "--T2part")
        quiet  && push!(args, rand(["--quiet", "-q"]))
        legacy && push!(args, "--legacy")

        for (param, paramval) in paramdict
            # Default flag/value pairs for `nothing` values; `nothing` is always default if allowable, therefore no flag/val is passed
            if paramval !== nothing
                push!(args, "--" * string(param)) # CLI flags are prepended with "--"
                append!(args,
                    paramval isa Tuple ? [string(x) for x in paramval] : # Pass each arg separately
                    paramval isa Bool  ? [] : # No arg necessary, flag only
                    [string(paramval)] # Pass string
                )
            end
        end

        return args

    elseif argstype === :mat
        #### MATLAB

        @assert legacy
        t2map_args  = T2map  ? Dict{Symbol,Any}() : nothing
        t2part_args = T2part ? Dict{Symbol,Any}() : nothing

        # Only these params are used within MATLAB (possibly spelled differently)
        mat_t2map_params  = [:Chi2Factor, :MinRefAngle, :nRefAngles, :nT2, :SetRefConAngle, :Reg, :SaveRegParam, :SetFlipAngle, :T1, :T2Range, :TE, :Threshold, :vTEparam]
        mat_t2part_params = [:T2Range, :SPWin, :MPWin, :Sigmoid]

        for (param, paramval) in paramdict
            T2map  && (param ∈ mat_t2map_params)  && jl_to_mat_param!(t2map_args,  param, paramval)
            T2part && (param ∈ mat_t2part_params) && jl_to_mat_param!(t2part_args, param, paramval)
        end

        return t2map_args, t2part_args

    elseif argstype === :jl
        #### Julia

        t2map_args  = T2map  ? Dict{Symbol,Any}(:legacy => legacy) : nothing
        t2part_args = T2part ? Dict{Symbol,Any}(:legacy => legacy) : nothing

        t2map_fields = DECAES.fieldsof(T2mapOptions, Set)
        t2part_fields = DECAES.fieldsof(T2partOptions, Set)
        for (param, paramval) in paramdict
            T2map  && (param ∈ t2map_fields)  && setindex!(t2map_args, paramval, param)
            T2part && (param ∈ t2part_fields) && setindex!(t2part_args, paramval, param)
        end

        return t2map_args, t2part_args
    end
end

function jl_to_mat_param!(opts, param, paramval)

    # T2mapSEcorr parameters which aren't in the MATLAB API
    new_t2map_params = Set{Symbol}([
        :SaveResidualNorm,
        :SaveDecayCurve,
        :SaveNNLSBasis,
        :Silent,
    ])

    if param == :SaveRegParam # renamed parameter
        opts[:Save_regparam] = ifelse(paramval, "yes", "no")
    elseif param == :nRefAngles # renamed parameter
        opts[:nAngles] = paramval
    elseif param == :SetRefConAngle # renamed parameter
        opts[:RefCon] = paramval
    elseif param == :Reg # renamed value: "no" => "none", "lcurve" => "gcv"
        opts[:Reg] = paramval == "none" ? "no" : paramval == "gcv" ? "lcurve" : paramval
    elseif param == :SPWin # renamed parameter
        opts[:spwin] = paramval
    elseif param == :MPWin # renamed parameter
        opts[:mpwin] = paramval
    elseif param ∉ new_t2map_params # skip Julia-only parameters
        opts[param] = paramval
    end

    return opts
end

field_error_string(x, y) = "max val = $(maximum(abs, y)), max diff = $(maximum(abs, x.-y)), rel diff = $(maximum(abs, (x.-y)./y))"

function test_field!(allpassed, x, y, prefix = "failed:"; kwargs...)
    passed = isapprox(x, y; kwargs..., nans = true)
    allpassed[] &= passed
    !passed && println(prefix * " (" * field_error_string(x,y) * ")")
    @test passed
end

# Compare t2map results for approximately equality
function test_compare_t2map(out1, out2; kwargs...)
    maps1, dist1, maps2, dist2 = out1..., out2...
    allpassed = Ref(true)
    for s in keys(maps1)
        haskey(maps2, s) && test_field!(allpassed, maps1[s], maps2[s], "maps failed: $s"; kwargs...)
    end
    test_field!(allpassed, dist1, dist2, "dist failed"; kwargs...)
    return allpassed[]
end

# Compare t2part results for approximately equality
function test_compare_t2part(part1, part2; kwargs...)
    allpassed = Ref(true)
    for s in keys(part1)
        haskey(part2, s) && test_field!(allpassed, part1[s], part2[s], "parts failed: $s"; kwargs...)
    end
    return allpassed[]
end

# CLI parameter settings to loop over
#   -Each param value will be tested individually, with all other params set to default values
#   -Each list should contain some non-default/edge case values
const cli_params_perms = Any[
    (:MPWin            .=> [(38e-3, 180e-3)],),
    (:MinRefAngle      .=> [55.0],),
    (:SetRefConAngle   .=> [172.0],),
    (
        :Reg           .=> ["none", "chi2", "gcv", "lcurve"],
        :Chi2Factor    .=> [nothing, 1.025, nothing, nothing],
    ),
    (:SPWin            .=> [(13e-3, 37e-3)],),
    (:SaveResidualNorm .=> [false, true],),
    (:SaveDecayCurve   .=> [false, true],),
    (:SaveNNLSBasis    .=> [false, true],),
    (:SaveRegParam     .=> [false, true],),
    (:SetFlipAngle     .=> [nothing, 170.0],),
    (:Sigmoid          .=> [nothing, 1.0],),
    (:T1               .=> [0.95],),
    (:T2Range          .=> [(16e-3, 1.8)],),
    (:TE               .=> [8e-3, 11e-3],),
    (:Threshold        .=> [125_000.0, Inf],), # Include non-zero and infinite (i.e. either some or all voxels skipped)
    (:nRefAngles       .=> [9, 10],), # Include odd number
    (:nRefAnglesMin    .=> [4, 5],), # Include odd number
    (:nT2              .=> [2, 47],), # Include odd number
]

@testset "CLI" begin
    make_settings_perms = [false, true]
    file_suffix_perms = [".mat", ".nii", ".nii.gz"] # Note: no PAR/REC or XML/REC, since we can't write to them
    iters = (cli_params_perms, make_settings_perms, file_suffix_perms)
    nloop = max(length.(iters)...)
    repeat_until(x) = Iterators.take(Iterators.cycle(x), nloop)

    for (param_val_lists, make_settings_file, file_suffix) in zip(map(repeat_until, iters)...), param_val_pairs in zip(param_val_lists...), legacy in [false, true]
        paramdict = deepcopy(default_paramdict)
        for (param, paramval) in param_val_pairs
            paramdict[param] = paramval
        end

        image = DECAES.mock_image(nTE = rand([4,5,20,37]))
        settings_kwargs_jl = Dict{Symbol,Any}(:argstype => :jl, :quiet => rand([true,false]), :legacy => legacy, :T2map => true, :T2part => true)
        settings_kwargs_cli = Dict{Symbol,Any}(:argstype => :cli, :quiet => rand([true,false]), :legacy => legacy, :T2map => true, :T2part => true)
        jl_t2map_kwargs, jl_t2part_kwargs = construct_args(paramdict; settings_kwargs_jl...)

        # Run T2map and T2part through Julia API for comparison
        t2map, t2dist = DECAES.tee_capture(suppress_terminal = true, suppress_logfile = true) do io
            T2mapSEcorr(io, image; jl_t2map_kwargs...)
        end
        t2part = DECAES.tee_capture(suppress_terminal = true, suppress_logfile = true) do io
            T2partSEcorr(io, t2dist; jl_t2part_kwargs...)
        end

        # Run CLI with both --T2map and --T2part flags
        mktempdir() do path
            settings_kwargs_cli[:outputpath] = path
            settings_kwargs_cli[:inputfilename] = joinpath(path, "input" * file_suffix)
            cli_t2map_args = construct_args(paramdict; settings_kwargs_cli...)

            t2maps_cli, t2dist_cli, t2parts_cli = DECAES.redirect_to_devnull() do
                run_main(image, cli_t2map_args; make_settings_file = make_settings_file)
            end
            t2map_passed = test_compare_t2map((t2map, t2dist), (t2maps_cli, t2dist_cli); rtol = 1e-14)
            t2part_passed = test_compare_t2part(t2part, t2parts_cli; rtol = 1e-14)
            if !(t2map_passed && t2part_passed)
                println("\n ------------------------------- \n")
                println("CLI with --T2map and --T2part failed: $param_val_pairs")
                @show jl_t2map_kwargs
                @show jl_t2part_kwargs
                @show cli_t2map_args
                @show paramdict
                println("\n ------------------------------- \n")
            end
        end

        # Run CLI with --T2part flag only
        mktempdir() do path
            settings_kwargs_cli[:outputpath] = path
            settings_kwargs_cli[:inputfilename] = joinpath(path, "input" * file_suffix)
            settings_kwargs_cli[:T2map] = false
            cli_t2part_args = construct_args(paramdict; settings_kwargs_cli...)

            t2maps_cli, t2dist_cli, t2parts_cli = DECAES.redirect_to_devnull() do
                run_main(t2dist, cli_t2part_args; make_settings_file = make_settings_file)
            end
            t2part_passed = test_compare_t2part(t2part, t2parts_cli; rtol = 1e-14)
            if !t2part_passed
                println("\n ------------------------------- \n")
                println("CLI with --T2part only failed: $param_val_pairs")
                @show jl_t2map_kwargs
                @show jl_t2part_kwargs
                @show cli_t2part_args
                @show paramdict
                println("\n ------------------------------- \n")
            end
        end
    end
end

# ================================================================================
# UBC MWI Toolbox MATLAB compatibility tests
#   NOTE: For these tests to run, MATLAB must be installed on your default path.
#   Additionally, the MWI NNLS toolbox (https://github.com/ubcmri/ubcmwf)
#   folder "MWI_NNLS_toolbox_0319" (and subfolders) must be added to your
#   default MATLAB path.
# ================================================================================

# Helper functions
matlabify(x::AbstractString) = String(x)
matlabify(x::AbstractArray) = Float64.(x)
matlabify(x::Tuple) = [Float64.(x)...]
matlabify(x::Bool) = x
matlabify(x) = map(Float64, x)
matlabify(kwargs::Base.Iterators.Pairs) = Iterators.flatten([(string(k), matlabify(v)) for (k,v) in kwargs if v !== nothing])

mfile_exists(fname) = MATLAB.mxcall(:exist, 1, fname) == 2

mxT2mapSEcorr(image, maxCores = 6; kwargs...) =
    MATLAB.mxcall(:T2map_SEcorr_nechoes_2019, 2, image, maxCores, matlabify(kwargs)...)

mxT2partSEcorr(image; kwargs...) =
    MATLAB.mxcall(:T2part_SEcorr_2019, 1, image, matlabify(kwargs)...)

# Arbitrary non-default T2mapSEcorr options for testing
const mat_t2map_params_perms = Any[
    (:TE               .=> [9e-3],),
    (:T1               .=> [1.1],),
    (:Threshold        .=> [250.0],),
    (:nT2              .=> [10, 59],), # Include odd number
    (:T2Range          .=> [(8e-3, 1.0)],),
    (:SetRefConAngle   .=> [175.0],),
    (:MinRefAngle      .=> [60.0],),
    (:nRefAngles       .=> [7, 12],),
    (
        :Reg           .=> ["none", "chi2", "lcurve", "gcv"],
        :Chi2Factor    .=> [nothing, 1.03, nothing, nothing],
    ),
    (:SetFlipAngle     .=> [178.0],),
    (:SaveResidualNorm .=> [false, true],),
    (:SaveDecayCurve   .=> [false, true],),
    (:SaveRegParam     .=> [false, true],),
    (:SaveNNLSBasis    .=> [false, true],),
]

# Arbitrary non-default T2partSEcorr options for testing
const mat_t2part_params_perms = Any[
    (:T2Range    .=> [(11e-3, 1.5)],),
    (:SPWin      .=> [(12e-3, 28e-3)],),
    (:MPWin      .=> [(35e-3, 150e-3)],),
    (:Sigmoid    .=> [1.5],),
]

function matlab_tests()
    # Relative tolerance threshold for legacy algorithms to match MATLAB version
    default_rtol = 1e-10

    @testset "T2mapSEcorr" begin
        settings_kwargs_jl = Dict{Symbol,Any}(:argstype => :jl, :quiet => rand([true,false]), :legacy => true, :T2map => true, :T2part => true)
        settings_kwargs_mat = Dict{Symbol,Any}(:argstype => :mat, :quiet => rand([true,false]), :legacy => true, :T2map => true, :T2part => true)

        for param_val_lists in mat_t2map_params_perms, param_val_pairs in zip(param_val_lists...)
            rtol = default_rtol
            paramdict = deepcopy(legacy_default_paramdict)
            for (param, paramval) in param_val_pairs
                paramdict[param] = paramval
            end

            # The MATLAB flag "lcurve" for choosing the regularization parameter uses the Generalized Cross-Validation (GCV) method.
            # There are several issues with comparing the DECAES implementation with MATLAB:
            #   1.  GCV involves minimizing a functional GCV(mu), which is implemented using an internal call to `fminbnd` with a tolerance of 1e-3,
            #       and therefore the Julia outputs would only match to at best a tolerance of 1e-3
            #   2.  There is an error in the `G(mu,C_g,d_g)` subfunction in the MATLAB MWI Toolbox:
            #       -   Numerator should be ||A*x_mu - b||^2, not ||A*x_mu - b||
            #       -   See e.g. equation (1.4) in Fenu, C. et al., 2017, GCV for Tikhonov regularization by partial SVD (https://doi.org/10.1007/s10543-017-0662-0)
            #   3.  There is an error in methodology:
            #       -   Solving regularized ||Ax-b||^2 via [A; mu*I] \ [b; 0] is equivalent to minimizing
            #               ||Ax-b||^2 + mu^2 * ||x||^2
            #           as opposed to
            #               ||Ax-b||^2 + mu * ||x||^2
            #           as is used in the MWI toolbox
            #   4.  The MATLAB equivalent to --Reg="lcurve" has been replaced with the more accurate --Reg="gcv".
            #       In DECAES, --Reg="lcurve" implements a different method which utilizes the L-curve method; see:
            #           A. Cultrera and L. Callegaro, “A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems”
            #           IOPSciNotes, vol. 1, no. 2, p. 025004, Aug. 2020, doi: 10.1088/2633-1357/abad0d
            # This comparison is therefore skipped by default. If you have a version in which these errors are fixed,
            # the below line can be modified (with rtol set appropriately in accordance with your solver tolerance)
            if any(param === :Reg && paramval ∈ ("lcurve", "gcv") for (param, paramval) in param_val_pairs)
                continue
                # rtol = 1e-3
            end

            jl_t2map_kwargs,  _ = construct_args(paramdict; settings_kwargs_jl...)
            mat_t2map_kwargs, _ = construct_args(paramdict; settings_kwargs_mat...)

            # Run T2mapSEcorr
            image = DECAES.mock_image(nTE = rand([4,5,20,37]))
            t2map_out_jl = DECAES.tee_capture(suppress_terminal = true, suppress_logfile = true) do io
                T2mapSEcorr(io, image; jl_t2map_kwargs...)
            end
            t2map_out_mat = DECAES.redirect_to_devnull() do
                mxT2mapSEcorr(image; mat_t2map_kwargs...)
            end
            allpassed = test_compare_t2map(t2map_out_jl, t2map_out_mat; rtol = rtol)
            if !allpassed
                println("\n ------------------------------- \n")
                println("MATLAB T2mapSEcorr comparison failed: $param_val_pairs")
                @show jl_t2map_kwargs
                @show mat_t2map_kwargs
                @show paramdict
                println("\n ------------------------------- \n")
            end
        end
    end

    @testset "T2partSEcorr" begin
        settings_kwargs_jl = Dict{Symbol,Any}(:argstype => :jl, :quiet => rand([true,false]), :legacy => true, :T2map => false, :T2part => true)
        settings_kwargs_mat = Dict{Symbol,Any}(:argstype => :mat, :quiet => rand([true,false]), :legacy => true, :T2map => false, :T2part => true)

        for param_val_lists in mat_t2part_params_perms, param_val_pairs in zip(param_val_lists...)
            paramdict = deepcopy(legacy_default_paramdict)
            delete!(paramdict, :nT2) # inferred
            for (param, paramval) in param_val_pairs
                paramdict[param] = paramval
            end
            _, jl_t2part_kwargs  = construct_args(paramdict; settings_kwargs_jl...)
            _, mat_t2part_kwargs = construct_args(paramdict; settings_kwargs_mat...)

            # Run T2partSEcorr
            T2dist = DECAES.mock_T2_dist(nT2 = rand([5,19,20,40,60])) #rand([4,5,40,41]))
            t2part_jl = DECAES.tee_capture(suppress_terminal = true, suppress_logfile = true) do io
                T2partSEcorr(io, T2dist; jl_t2part_kwargs...)
            end
            t2part_mat = DECAES.redirect_to_devnull() do
                mxT2partSEcorr(T2dist; mat_t2part_kwargs...)
            end
            allpassed = test_compare_t2part(t2part_jl, t2part_mat; rtol = default_rtol)
            if !allpassed
                println("\n ------------------------------- \n")
                println("MATLAB T2partSEcorr comparison failed: $param_val_pairs")
                @show jl_t2part_kwargs
                @show mat_t2part_kwargs
                @show paramdict
                println("\n ------------------------------- \n")
            end
        end
    end
end

# Try loading MATLAB.jl and running tests
try
    # Check for environment flags
    mwi_toolbox_path = get(ENV, "DECAES_MWI_TOOLBOX_PATH", "")
    run_matlab_tests = get(ENV, "DECAES_RUN_MWI_TOOLBOX_TESTS", "") != "0"

    if run_matlab_tests
        @eval using MATLAB

        if !isempty(mwi_toolbox_path)
            mxcall(:addpath, 0, mwi_toolbox_path)
        end

        if mfile_exists("T2map_SEcorr_nechoes_2019") && mfile_exists("T2part_SEcorr_2019")
            matlab_tests()
        else
            @warn "Files T2map_SEcorr_nechoes_2019.m and T2part_SEcorr_2019.m were not found on the default MATLAB path. " *
                "Modify your default MATLAB path to include these files, or set the DECAES_MWI_TOOLBOX_PATH environment variable.\n\n" *
                "For example, on unix-like systems run" *
                "\n\n    export DECAES_MWI_TOOLBOX_PATH=/path/to/MWI_NNLS_toolbox_0319\n\n" *
                "before testing DECAES, or add a command such as" *
                "\n\n    addpath /path/to/MWI_NNLS_toolbox_0319\n\n" *
                "to your startup.m file in MATLAB."
        end
    end
catch e
    @warn "Failed to load Julia package MATLAB.jl; skipping UBCMWF MATLAB tests"
    @warn sprint(showerror, e, catch_backtrace())
end

nothing


text:
using BlogPlots
using Test

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


text:
using DataFrames, DataFramesMeta, Statistics
using PyPlot

# Top-level functions to process dataframe files (outputs of simulations)
# in order to create the csv data for the plots
# in Figure 8 and Figure 9.
# (the output .png files were not used in the paper, only for inspection)

# For Figure 9
function process_robustfrombayes_M(;savepath = homedir(),
                filepath = "/path/to/jldfiles_MULTI_GRIDROB/",
                learnerbase = "gridrobust_multi_tbayes_pc.jld2")
   jldfiles = searchdirendings(filepath, ".jld2")
   clist = clisthardcoded_robM()
   # THE FINALISTS ARE:
   #stochvalue = 5; pcstar = 0.004124626382901352
   #stochvalue = 0.14; pcstar = 0.004124626382901352
   #stochvalue = 5; pcstar = 0.04124626382901352
   stochvalue = 0.14; pcstar = 0.04124626382901352

   file = load(joinpath(filepath, learnerbase))
   learnernamebase = file["learnertype"]
   dfbase = file["df"]
   dfbase = @where(dfbase, :stochasticity .== stochvalue)

   dfresults = DataFrame()
   fig = figure(); ax = gca()
   for i in 1:length(jldfiles)
       file = load(joinpath(filepath, jldfiles[i]))
       learnername = file["learnertype"]
       df = file["df"]
       df = @where(df, :stochasticity .== stochvalue)
       indexminimum = argmin(@where(df,:changeprobability .== pcstar)[:meanterrors])
       thetastar = @where(df, :changeprobability .== pcstar)[:param][indexminimum]  # theta of min error for central pcstar. Ie optimal theta for central pc

       changeprobs = unique(df[:changeprobability])
       data = []
       for j in 1:length(changeprobs)
           error_thetastar_pc = @where(df,
                            :changeprobability .== changeprobs[j],
                            :param .== thetastar)[:meanterrors][1] # Error when using thetastar_pc_base of jth pc
           error_theta_pc_base = minimum(@where(dfbase,
                            :changeprobability .== changeprobs[j])[:meanterrors])
           err_diffs = error_thetastar_pc - error_theta_pc_base
           push!(data, err_diffs)
        end
        labelstring = makelabelstring(learnername, jldfiles[i])
        ax.plot(changeprobs, data, label = labelstring, color = clist[i])

        if i==1; dfresults[:x] = changeprobs; end
        dfresults[Symbol(labelstring)] = data
   end
   ax.legend(); ax.grid(true); xlabel("Pc");
   ylabel("Error diff from "* makelabelstring(learnernamebase, learnerbase))
   pcstarformatted = round(pcstar, digits=4)
   title("Robustness plot" * "stoch: " * string(stochvalue)
        * " , Pc: " * string(pcstarformatted))
   PyPlot.savefig(joinpath(savepath,
            "multi_robfrombayes_s"* string(stochvalue) * "_Pc" *string(pcstarformatted)*".png"))
    CSV.write(joinpath(savepath,
            "multi_robustnessfrombayes_s"* string(stochvalue) * "_Pc" *string(pcstarformatted)*".csv"),
            dfresults, delim = " ");
end
# For Figure 8
function process_robustfrombayes_G(;savepath = homedir(),
                filepath = "/path/to/jldfiles_GAUSSIAN_GRIDROB/",
                learnerbase = "gridrobust_gaussian_gbayes_pc.jld2")
   jldfiles = searchdirendings(filepath, ".jld2")
   clist = clisthardcoded_robG()
   # THE FINALISTS ARE:
   #sigmavalue = 5; pcstar = 0.004124626382901352
   #sigmavalue = 0.1; pcstar = 0.004124626382901352
   #sigmavalue = 5; pcstar = 0.04124626382901352
   sigmavalue = 0.1; pcstar = 0.04124626382901352

   file = load(joinpath(filepath, learnerbase))
   learnernamebase = file["learnertype"]
   dfbase = file["df"]
   dfbase = @where(dfbase, :sigma .== sigmavalue)

   dfresults = DataFrame()
   fig = figure(); ax = gca()
   for i in 1:length(jldfiles)
       file = load(joinpath(filepath, jldfiles[i]))
       learnername = file["learnertype"]
       df = file["df"]
       df = @where(df, :sigma .== sigmavalue)
       minerr = minimum(@where(df,
                        :changeprobability .== pcstar)[:meanmerrors])
        @show minerr
       indexminimum = argmin(@where(df,:changeprobability .== pcstar)[:meanmerrors])
       thetastar = @where(df, :changeprobability .== pcstar)[:param][indexminimum]  # theta of min error for central pcstar. Ie optimal theta for central pc
       @show thetastar

       changeprobs = unique(df[:changeprobability])
       data = []
       for j in 1:length(changeprobs)
           error_thetastar_pc = @where(df,
                            :changeprobability .== changeprobs[j],
                            :param .== thetastar)[:meanmerrors][1] # Error when using thetastar_pc_base of jth pc
           error_theta_pc_base = minimum(@where(dfbase,
                            :changeprobability .== changeprobs[j])[:meanmerrors])
           err_diffs = error_thetastar_pc - error_theta_pc_base
           push!(data, err_diffs)
        end
        labelstring = makelabelstring(learnername, jldfiles[i])
        ax.plot(changeprobs, data, label = labelstring, color = clist[i])

        if i==1; dfresults[:x] = changeprobs; end
        dfresults[Symbol(labelstring)] = data
   end
   ax.legend(); ax.grid(true); xlabel("Pc");
   ylabel("Error diff from "* makelabelstring(learnernamebase, learnerbase))
   ax.set_xscale("log")
   pcstarformatted = round(pcstar, digits=4)
   title("Robustness plot - Gaussian task - " * "sigma: " * string(sigmavalue)
        * " , Pc: " * string(pcstarformatted))
   PyPlot.savefig(joinpath(savepath,
            "gauss_robfrombayes_s"* string(sigmavalue) * "_Pc" *string(pcstarformatted)*".png"))
   CSV.write(joinpath(savepath, "gaussian_robustnessfrombayes_s"* string(sigmavalue) * "_Pc" *string(pcstarformatted)*".csv"), dfresults,
                        delim = " ");
end

function makelabelstring(learnername, jldfile)
    if occursin("NonstationaryModelLearning", learnername)
        substrings = split(learnername,".")
        labelstring = substrings[end]
    else
        labelstring = learnername
    end
    if occursin("1", jldfile)
        if occursin("10", jldfile)
            labelstring = labelstring*"10"
        else
            labelstring = labelstring*"1"
        end
    elseif occursin("20", jldfile)
        labelstring = labelstring*"20"
    end
    labelstring
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------- Not used eventually for the paper, only for inspection ---------------
function process_robust_M(;savepath = homedir(),
                filepath = "/path/to/jldfiles_MULTI_GRIDROB/")

   jldfiles = searchdirendings(filepath, ".jld2")
   clist = clisthardcoded_robM()
   # THE FINALISTS ARE:
   #stochvalue = 5; pcstar = 0.004124626382901352 # for log scale
   #stochvalue = 0.14; pcstar = 0.004124626382901352 # for log scale
   #stochvalue = 5; pcstar = 0.04124626382901352 # for log scale
   stochvalue = 0.14; pcstar = 0.04124626382901352 # for log scale

   dfresults = DataFrame()
   fig = figure(); ax = gca()
   for i in 1:length(jldfiles)
       file = load(joinpath(filepath, jldfiles[i]))
       learnername = file["learnertype"]
       df = file["df"]
       df = @where(df, :stochasticity .== stochvalue)
       error_thetastar_pcstar = minimum(@where(df,
                            :changeprobability .== pcstar)[:meanterrors])
       indexminimum = argmin(@where(df,
                            :changeprobability .== pcstar)[:meanterrors])
       thetastar = @where(df,
                    :changeprobability .== pcstar)[:param][indexminimum]
       changeprobs = unique(df[:changeprobability])
       data = []
       for j in 1:length(changeprobs)
           error_thetastar_pc = @where(df,
                            :changeprobability .== changeprobs[j],
                            :param .== thetastar)[:meanterrors][1]
           error_theta_pc = minimum(@where(df,
                            :changeprobability .== changeprobs[j])[:meanterrors])
           err_diffs = error_thetastar_pc - error_theta_pc #+ error_thetastar_pcstar
           push!(data, err_diffs)
        end
        labelstring = makelabelstring(learnername, jldfiles[i])
        ax.plot(changeprobs, data, label = labelstring, color = clist[i])
        if i==1; dfresults[:x] = changeprobs; end
        dfresults[Symbol(labelstring)] = data
   end
   ax.legend(); ax.grid(true); xlabel("Pc"); ylabel("Error diff")
   pcstarformatted = round(pcstar, digits=4)
   title("Robustness plot" * "stoch: " * string(stochvalue)
        * " , Pc: " * string(pcstarformatted))
   PyPlot.savefig(joinpath(savepath,
            "multi_rob_s"* string(stochvalue) * "_Pc" *string(pcstarformatted)*".png"))
    CSV.write(joinpath(savepath, "multi_robustness_s"* string(stochvalue) * "_Pc" *string(pcstarformatted)* ".csv"),
            dfresults, delim = " ");
end
function process_robust_G(;savepath = homedir(),
        filepath = "/path/to/jldfiles_GAUSSIAN_GRIDROB/"
        )
   jldfiles = searchdirendings(filepath, ".jld2")
   clist = clisthardcoded_robG()
   # THE FINALISTS ARE:
   sigmavalue = 5; pcstar = 0.004124626382901352
   #sigmavalue = 0.1; pcstar = 0.004124626382901352
   #sigmavalue = 5; pcstar = 0.04124626382901352
   # sigmavalue = 0.1; pcstar = 0.04124626382901352

   dfresults = DataFrame()
   fig = figure(); ax = gca()
   for i in 1:length(jldfiles)
       file = load(joinpath(filepath, jldfiles[i]))
       learnername = file["learnertype"]
       df = file["df"]
       df = @where(df, :sigma .== sigmavalue)
       error_thetastar_pcstar = minimum(@where(df,
                            :changeprobability .== pcstar)[:meanmerrors])
       indexminimum = argmin(@where(df,
                            :changeprobability .== pcstar)[:meanmerrors])
       thetastar = @where(df,
                    :changeprobability .== pcstar)[:param][indexminimum]
       changeprobs = unique(df[:changeprobability])
       data = []
       for j in 1:length(changeprobs)
           error_thetastar_pc = @where(df,
                            :changeprobability .== changeprobs[j],
                            :param .== thetastar)[:meanmerrors][1]
           error_theta_pc = minimum(@where(df,
                            :changeprobability .== changeprobs[j])[:meanmerrors])
           err_diffs = error_thetastar_pc - error_theta_pc #+ error_thetastar_pcstar
           push!(data, err_diffs)
        end
        labelstring = makelabelstring(learnername, jldfiles[i])
        ax.plot(changeprobs, data, label = labelstring, color = clist[i])
        if i==1; dfresults[:x] = changeprobs; end
        dfresults[Symbol(labelstring)] = data
   end
   ax.legend(); ax.grid(true); xlabel("Pc"); ylabel("Error diff")
   pcstarformatted = round(pcstar, digits=4)
   title("Robustness plot - Gaussian task - " * "sigma: " * string(sigmavalue)
        * " , Pc: " * string(pcstarformatted))
   PyPlot.savefig(joinpath(savepath,
            "gauss_rob_s"* string(sigmavalue) * "_Pc" *string(pcstarformatted)*".png"))
   CSV.write(joinpath(savepath, "gaussian_robustness_s"* string(sigmavalue) * "_Pc" *string(pcstarformatted)* "OLD_linear" * ".csv"),
                    dfresults, delim = " ");
end


text:
type BitVectorRepresentation <: Representation
  length::Int

  BitVectorRepresentation(length::Int) =
    new(length)
end

"""
Provides a definition for a bit vector representation.
"""
type BitVectorRepresentationDefinition <: RepresentationDefinition
  length::Int

  BitVectorRepresentationDefinition(length::Int) =
    new(length)
end

"""
Composes a bit vector representation from its provided definition.
"""
compose!(r::BitVectorRepresentationDefinition) =
  BitVectorRepresentation(r.length)

"""
Bit vector representations operate on a fixed-length vector of binary values,
implemented as a vector of Int objects (for now).

**Parameters:**
* `length::Int`, the length of this bit vector. Defaults to `100` if no value
is supplied.
"""
bit_vector() = bit_vector(100)
bit_vector(length::Int) = BitVectorRepresentationDefinition(length)

minimum_value(::BitVectorRepresentation) = 0
maximum_value(::BitVectorRepresentation) = 1
codon_type(::BitVectorRepresentation) = Int
chromosome(::BitVectorRepresentation) = Vector{Int}

"""
Generates a bit vector at random using its representation.
"""
generate(r::BitVectorRepresentation) = Base.rand(0:1, r.length)

# Converts a given bit vector chromosome into an int vector chromosome, destructively.
function convert!(r_from::BitVectorRepresentation,
  r_to::IntVectorRepresentation,
  c_from::IndividualStage{Vector{Int}},
  c_to::IndividualStage{Vector{Int}}
)
  c_to.value = Array(Int, r_to.length)
  for i in 1:r_to.length
    j = ((i - 1) * 8) + 1
    c_to.value[i] = parseint(join(c_from.value[j:j+7]), 2)
  end
end

function convert!(r_from::IntVectorRepresentation,
  r_to::BitVectorRepresentation,
  c_from::IndividualStage{Vector{Int}},
  c_to::IndividualStage{Vector{Int}}
)
  c_to.value = Array(Int, r_to.length)
  for i in 1:r_from.length
    offset = ((i - 1) * 8)
    binary = bin(c_from.value[i], 8)
    for p in 1:8
      j = offset + p
      c_to.value[j] = int(binary[p:p])
    end
  end
end


text:
# Obsoleted samplers

# These samplers are not actually used by anyone, and they have not been properly tested.
# The codes are still kept here in case we want to pick them up in future.
#

##### Draw Table #####

# Store an alias table
struct DiscreteDistributionTable <: Sampler{Univariate,Discrete}
    table::Vector{Vector{Int64}}
    bounds::Vector{Int64}
end

# TODO: Test if bit operations can speed up Base64 mod's and fld's
function DiscreteDistributionTable(probs::Vector{T}) where T <: Real
    # Cache the cardinality of the outcome set
    n = length(probs)

    # Convert all AbstractFloat's into integers
    vals = Vector{Int64}(undef, n)
    for i in 1:n
        vals[i] = round(Int, probs[i] * 64^9)
    end

    # Allocate digit table and digit sums as table bounds
    table = Vector{Vector{Int64}}(undef, 9)
    bounds = zeros(Int64, 9)

    # Special case for deterministic distributions
    for i in 1:n
        if vals[i] == 64^9
            table[1] = Vector{Int64}(undef, 64)
            for j in 1:64
                table[1][j] = i
            end
            bounds[1] = 64^9
            for j in 2:9
                table[j] = Vector{Int64}()
                bounds[j] = 64^9
            end
            return DiscreteDistributionTable(table, bounds)
        end
    end

    # Fill tables
    multiplier = 1
    for index in 9:-1:1
        counts = Vector{Int64}()
        for i in 1:n
            digit = mod(vals[i], 64)
            # vals[i] = fld(vals[i], 64)
            vals[i] >>= 6
            bounds[index] += digit
            for itr in 1:digit
                push!(counts, i)
            end
        end
        bounds[index] *= multiplier
        table[index] = counts
        # multiplier *= 64
        multiplier <<= 6
    end

    # Make bounds cumulative
    bounds = cumsum(bounds)

    return DiscreteDistributionTable(table, bounds)
end

function rand(table::DiscreteDistributionTable)
    # 64^9 - 1 == 0x003fffffffffffff
    i = rand(1:(64^9 - 1))
    # if i == 64^9
    #   return table.table[9][rand(1:length(table.table[9]))]
    # end
    bound = 1
    while i > table.bounds[bound] && bound < 9
        bound += 1
    end
    if bound > 1
        index = fld(i - table.bounds[bound - 1] - 1, 64^(9 - bound)) + 1
    else
        index = fld(i - 1, 64^(9 - bound)) + 1
    end
    return table.table[bound][index]
end

Base.show(io::IO, table::DiscreteDistributionTable) = @printf io "DiscreteDistributionTable"

##### Huffman Table ######

abstract type HuffmanNode{T} <: Sampler{Univariate,Discrete} end

struct HuffmanLeaf{T} <: HuffmanNode{T}
    value::T
    weight::UInt64
end

struct HuffmanBranch{T} <: HuffmanNode{T}
    left::HuffmanNode{T}
    right::HuffmanNode{T}
    weight::UInt64
end
HuffmanBranch(ha::HuffmanNode{T},hb::HuffmanNode{T}) where {T} = HuffmanBranch(ha, hb, ha.weight + hb.weight)

Base.isless(ha::HuffmanNode{T}, hb::HuffmanNode{T}) where {T} = isless(ha.weight,hb.weight)
Base.show(io::IO, t::HuffmanNode{T}) where {T} = show(io,typeof(t))

function Base.getindex(h::HuffmanBranch{T},u::UInt64) where T
    while isa(h,HuffmanBranch{T})
        if u < h.left.weight
            h = h.left
        else
            u -= h.left.weight
            h = h.right
        end
    end
    h.value
end

# build the huffman tree
# could be slightly more efficient using a Deque.
function huffman(values::AbstractVector{T},weights::AbstractVector{UInt64}) where T
    leafs = [HuffmanLeaf{T}(values[i],weights[i]) for i = 1:length(weights)]
    sort!(leafs; rev=true)

    branches = Vector{HuffmanBranch{T}}()

    while !isempty(leafs) || length(branches) > 1
        left = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        right = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        pushfirst!(branches,HuffmanBranch(left,right))
    end

    pop!(branches)
end

function rand(h::HuffmanNode{T}) where T
    w = h.weight
    # generate uniform UInt64 objects on the range 0:(w-1)
    # unfortunately we can't use Range objects, as they don't have sufficient length
    u = rand(UInt64)
    if (w & (w-1)) == 0
        # power of 2
        u = u & (w-1)
    else
        m = typemax(UInt64)
        lim = m - (rem(m,w)+1)
        while u > lim
            u = rand(UInt64)
        end
        u = rem(u,w)
    end
    h[u]
end


text:
function calc_rho( pw::PWGrid, Focc, psi::Array{ComplexF64,2} )
    Ω  = pw.Ω
    Ns = pw.Ns
    Npoints = pw.Npoints
    Nstates = size(psi)[2]
    #
    rho = zeros(ComplexF64,Npoints)

    # Transform to real space
    psiR = G_to_R(Ns,psi)

    # orthonormalization in real space
    ortho_gram_schmidt!( Nstates, psiR )
    psiR = sqrt(Npoints/Ω)*psiR

    for is = 1:Nstates
        for ip = 1:Npoints
            rho[ip] = rho[ip] + conj(psiR[ip,is])*psiR[ip,is]*Focc[is]
        end
    end
    return real(rho)
end


text:
println("Starting precompilation")
using OrdinaryDiffEq
println("Starting tests")
using Flux, DiffEqFlux
using Test
using Distributions
using Distances
using Zygote
using DistributionsAD
using LinearAlgebra

##callback to be used by all tests
function cb(p,l)
    @show l
    false
end

###
#test for default base distribution and deterministic trace CNF
###

nn = Chain(Dense(1, 1, tanh))
data_train = [Float32(rand(Beta(7,7))) for i in 1:100]
tspan = (0.0,10.0)
cnf_test = DeterministicCNF(nn,tspan,Tsit5())

function loss_adjoint(θ)
    logpx = [cnf_test(x,θ) for x in data_train]
    loss = -mean(logpx)
end

res = DiffEqFlux.sciml_train(loss_adjoint, 0.01.*cnf_test.p,
                                        ADAM(0.01), cb=cb,
                                        maxiters = 100)

θopt = res.minimizer
data_validate = [Float32(rand(Beta(7,7))) for i in 1:100]
actual_pdf = [pdf(Beta(7,7),r) for r in data_validate]
#use direct trace calculation for predictions
learned_pdf = [exp(cnf_test(r,θopt)) for r in data_validate]

@test totalvariation(learned_pdf, actual_pdf)/100 < 0.25

###
#test for alternative base distribution and deterministic trace CNF
###

nn = Chain(Dense(1, 3, tanh), Dense(3, 1, tanh))
data_train = [Float32(rand(Normal(6.0,0.7))) for i in 1:100]
tspan = (0.0,10.0)
cnf_test = DeterministicCNF(nn,tspan,Tsit5(),basedist=MvNormal([0.0],[2.0]))

function loss_adjoint(θ)
    logpx = [cnf_test(x,θ) for x in data_train]
    loss = -mean(logpx)
end

res = DiffEqFlux.sciml_train(loss_adjoint, 0.01.*cnf_test.p,
                                          ADAM(0.01), cb = cb,
                                          maxiters = 300)

θopt = res.minimizer
data_validate = [Float32(rand(Normal(6.0,0.7))) for i in 1:100]
actual_pdf = [pdf(Normal(6.0,0.7),r) for r in data_validate]
learned_pdf = [exp(cnf_test(r, θopt)) for r in data_validate]

@test totalvariation(learned_pdf, actual_pdf)/100 < 0.25

###
#test for multivariate distribution and deterministic trace CNF
###

nn = Chain(Dense(2, 2, tanh))
μ = ones(2)
Σ = 7*I + zeros(2,2)
mv_normal = MvNormal(μ, Σ)
data_train = [Float32.(rand(mv_normal)) for i in 1:100]
tspan = (0.0,10.0)
cnf_test = DeterministicCNF(nn,tspan,Tsit5())

function loss_adjoint(θ)
    logpx = [cnf_test(x,θ) for x in data_train]
    loss = -mean(logpx)
end

res = DiffEqFlux.sciml_train(loss_adjoint, cnf_test.p,
                                          ADAM(0.1), cb = cb,
                                          maxiters = 300)

θopt = res.minimizer
data_validate = [Float32.(rand(mv_normal)) for i in 1:100]
actual_pdf = [pdf(mv_normal,r) for r in data_validate]
learned_pdf = [exp(cnf_test(r, θopt)) for r in data_validate]

@test totalvariation(learned_pdf, actual_pdf)/100 < 0.25

###
#test for default multivariate distribution and FFJORD with regularizers
###

nn = Chain(Dense(1, 1, tanh))
data_train = [Float32(rand(Beta(7,7))) for i in 1:100]
tspan = (0.0,10.0)
ffjord_test = FFJORD(nn,tspan,Tsit5())

function loss_adjoint(θ)
    logpx = [-ffjord_test(x,θ,true)[1] + 0.1*(ffjord_test(x,θ,true)[2]+ffjord_test(x,θ,true)[3]) for x in data_train]
    loss = mean(logpx)
end
res = DiffEqFlux.sciml_train(loss_adjoint, 0.01.*ffjord_test.p,
                                        ADAM(0.01), cb=cb,
                                        maxiters = 100)

θopt = res.minimizer
data_validate = [Float32(rand(Beta(7,7))) for i in 1:100]
actual_pdf = [pdf(Beta(7,7),r) for r in data_validate]
#use direct trace calculation for predictions
learned_pdf = [exp(ffjord_test(r,θopt,false,false)[1]) for r in data_validate]

@test totalvariation(learned_pdf, actual_pdf)/100 < 0.25


text:
isdefined(Base, :__precompile__) && __precompile__()

module Einsum

using Base.Cartesian
export @einsum, @einsimd

macro einsum(ex)
    _einsum(ex)
end

macro einsimd(ex)
    _einsum(ex,true,true)
end

macro einsum_checkinbounds(ex)
    _einsum(ex,false)
end

function _einsum(ex::Expr, inbound=true, simd=false)

    # Get left hand side (lhs) and right hand side (rhs) of equation
    lhs = ex.args[1]
    rhs = ex.args[2]

    # Get info on the left-hand side
    lhs_idx,lhs_arr,lhs_dim = get_indices!(lhs)
    length(lhs_arr) != 1 && throw(ArgumentError("Left-hand side of ",
        "equation contains multiple arguments. Only a single referencing ",
        " expression (e.g. @einsum A[i] = ...) should be used,"))

    # Get info on the right-hand side
    rhs_idx,rhs_arr,rhs_dim = get_indices!(rhs)

    # remove duplicate indices on left-hand and right-hand side
    # and ensure that the array sizes match along these dimensions
    ###########################################################
    ex_check_dims = :()

    # remove duplicate indices on the right hand side
    for i in reverse(1:length(rhs_idx))
        duplicated = false
        di = rhs_dim[i]
        for j = 1:(i-1)
            if rhs_idx[j] == rhs_idx[i]
                # found a duplicate
                duplicated = true
                dj = rhs_dim[j]

                # add dimension check ensuring consistency
                ex_check_dims = quote
                    @assert $(esc(dj)) == $(esc(di))
                    $ex_check_dims
                end
            end
        end
        for j = 1:length(lhs_idx)
            if lhs_idx[j] == rhs_idx[i]
                dj = lhs_dim[j]
                if ex.head == :(:=)
                    # ex.head is :=
                    # infer the size of the lhs array
                    lhs_dim[j] = di
                else
                    # ex.head is =, +=, *=, etc.
                    lhs_dim[j] = :(min($dj,$di))
                end
                duplicated = true
            end
        end
        if duplicated
            deleteat!(rhs_idx,i)
            deleteat!(rhs_dim,i)
        end
    end

    # remove duplicate indices on the left hand side
    for i in reverse(1:length(lhs_idx))
        duplicated = false
        di = lhs_dim[i]

        # don't need to check rhs, already done above

        for j = 1:(i-1)
            if lhs_idx[j] == lhs_idx[i]
                # found a duplicate
                duplicated = true
                dj = lhs_dim[j]

                # add dimension check
                ex_check_dims = quote
                    @assert $(esc(dj)) == $(esc(di))
                    $ex_check_dims
                end
            end
        end
        if duplicated
            deleteat!(lhs_idx,i)
            deleteat!(lhs_dim,i)
        end
    end

    # Create output array if specified by user
    ex_get_type = :(nothing)
    ex_create_arrays = :(nothing)
    ex_assignment_op = :(=)

    if ex.head == :(:=)

        # infer type of allocated array
        #    e.g. rhs_arr = [:A,:B]
        #    then the following line produces :(promote_type(eltype(A),eltype(B)))
        rhs_type = Expr(:call,:promote_type, [ Expr(:call,:eltype,arr) for arr in rhs_arr ]...)

        ex_get_type = :($(esc(:(local T = $rhs_type))))
        if length(lhs_dim) > 0
            ex_create_arrays = :($(esc(:($(lhs_arr[1]) = Array{$rhs_type}($(lhs_dim...))))))
        else
            ex_create_arrays = :($(esc(:($(lhs_arr[1]) = zero($rhs_type)))))
        end
    else
        ex_get_type = :($(esc(:(local T = eltype($(lhs_arr[1]))))))
        ex_create_arrays = :(nothing)
        ex_assignment_op = ex.head
    end

    # Copy equation, ex is the Expr we'll build up and return.
    unquote_offsets!(ex)

    if length(rhs_idx) > 0
        # There are indices on rhs that do not appear in lhs.
        # We sum over these variables.

        # Innermost expression has form s += rhs
        ex.args[1] = :s
        ex.head = :(+=)
        ex = esc(ex)

        # Nest loops to iterate over the summed out variables
        ex = nest_loops(ex,rhs_idx,rhs_dim,simd)

        lhs_assignment = Expr(ex_assignment_op, lhs, :s)
        # Prepend with s = 0, and append with assignment
        # to the left hand side of the equation.
        ex = quote
            $(esc(:(local s = zero(T))))
            $ex
            $(esc(lhs_assignment))
        end
    else
        # We do not sum over any indices
        # ex.head = :(=)
        ex.head = ex_assignment_op
        ex = :($(esc(ex)))
    end

    # Next loops to iterate over the destination variables
    ex = nest_loops(ex,lhs_idx,lhs_dim)

    # Assemble full expression and return
    if inbound
        return quote
            $ex_create_arrays
            let
                @inbounds begin
                    $ex_check_dims
                    $ex_get_type
                    $ex
                end
            end
        end
    else
        return quote
            $ex_create_arrays
            let
                $ex_check_dims
                $ex_get_type
                $ex
            end
        end
    end
end

function nest_loops(ex::Expr,idx::Vector{Symbol},dim::Vector{Expr},simd=false)
    if simd && !isempty(idx)
        # innermost index and dimension
        i = idx[1]
        d = dim[1]

        # Add @simd to the innermost loop.
        ex = quote
            local $(esc(i))
            @simd for $(esc(i)) = 1:$(esc(d))
                $(ex)
            end
        end
        start_ = 2
    else
        start_ = 1
    end

    # Add remaining for loops
    for j = start_:length(idx)
        # index and dimension we are looping over
        i = idx[j]
        d = dim[j]

        # add for loop around expression
        ex = quote
            local $(esc(i))
            for $(esc(i)) = 1:$(esc(d))
                $(ex)
            end
        end
    end
    return ex
end

function get_indices!(
        ex::Symbol,
        idx_store=Symbol[],
        arr_store=Symbol[ex],
        dim_store=Expr[]
    )
    return idx_store,arr_store,dim_store
end

@inline get_indices!(ex::Number,args...) = (args...)

function get_indices!(
        ex::Expr,
        idx_store=Symbol[],
        arr_store=Symbol[],
        dim_store=Expr[]
    )

    if ex.head == :ref
        # e.g. A[i,j,k] #
        push!(arr_store, ex.args[1])

        # iterate over indices (e.g. i,j,k)
        for (i,arg) in enumerate(ex.args[2:end])

            if typeof(arg) == Symbol
                # e.g. A[i]
                #    First, push :i to index list
                #    Second, push size(A,1) to dimension list
                push!(idx_store,arg)
                push!(dim_store,:(size($(ex.args[1]),$i)))

            elseif typeof(arg) <: Number
                # e.g. A[5]
                #    Do nothing, since we don't iterate over this dimension
                continue
            else
                # e.g. A[i+:offset] or A[i+5]
                #    arg is an Expr in this case
                #    We restrict it to be a Symbol (e.g. :i) followed by either
                #        a number or quoted expression.
                #    As before, push :i to index list
                #    Need to add/subtract off the offset to dimension list
                arg.head == :quote && continue
                @assert arg.head == :call
                @assert length(arg.args) == 3
                op = arg.args[1]
                sym = arg.args[2]
                offT = typeof(arg.args[3])
                if offT <: Integer
                    off = arg.args[3]::Integer
                elseif offT <: Expr && arg.args[3].head == :quote
                    off = arg.args[3].args[1]::Symbol
                elseif offT == QuoteNode
                    off = arg.args[3].value::Symbol
                else
                    throw(ArgumentError("improper expression inside reference on rhs"))
                end
                @assert typeof(sym) == Symbol

                # push :i to indices we're iterating over
                push!(idx_store, sym)

                # need to invert + or - to determine iteration range
                if op == :+
                    push!(dim_store,:( (size($(ex.args[1]),$i) - $off )))
                elseif op == :-
                    push!(dim_store,:( (size($(ex.args[1]),$i) + $off )))
                else
                    throw(ArgumentError("operations inside ref on rhs are limited to + or -"))
                end
            end
        end
    else
        # e.g. 2*A[i,j] or transpose(A[i,j])
        @assert ex.head == :call
        for arg in ex.args[2:end]
            get_indices!(arg,idx_store,arr_store,dim_store)
        end
    end
    idx_store,arr_store,dim_store
end

function unquote_offsets!(ex::Expr, inside_ref=false)
    inside_ref = inside_ref || ex.head == :ref
    for i = 1:length(ex.args)
        if isa(ex.args[i], Expr)
            if ex.args[i].head == :quote && inside_ref
                ex.args[i] = :($(ex.args[i].args[1]))
            else
                unquote_offsets!(ex.args[i], inside_ref)
            end
        end
    end
    return ex
end

# end module
############
end


text:
module EUPower

function SayHi()
    print("Hello")
    return 1
end

end


text:
include("QuadLoss.jl")
include("LogisticLoss.jl")
include("MultinomialLoss.jl")
include("OvALoss.jl")
include("MultinomialOrdinalLoss.jl")
include("OrdisticLoss.jl")
include("BvSLoss.jl")


text:
# Data loading for example
"""
    load_brexit()

Returns a `DataFrame` with quarterly GDP data on 30 OECD countries, borrowed from the analysis of
the effect of Brexit on UK GDP undertaken by the [Centre for European Reform](
https://www.cer.eu/insights/cost-brexit-june-2019)

With the Brexit referendum having taken place on 23 June 2016, the `TreatmentPanel` for Brexit
should be specified as
    ```
    TreatmentPanel(df, "United Kingdom" => Date(2016, 7, 1);
        outcome = :realgdp, id_var = :country, t_var = :quarter)
    ```
"""
function load_brexit()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","brexit.csv"), DataFrame; stringtype = String)
end

"""
    load_brexit_panel()

Returns a `BalancedPanel` object based on the data set available through load_brexit(), with
treatment assignment specified as `"United Kingdom" => Date(2016, 7, 1)`

"""
function load_brexit_panel()
    df_brexit = load_brexit()
    BalancedPanel(df_brexit, "United Kingdom" => Date(2016, 7, 1);
        id_var = :country, t_var = :quarter, outcome_var = :realgdp)
end

"""
    load_germany()

    Returns a `DataFrame` with annual observations of per-capita GDP and covariates between 1960 and
    2003. The data is used in Abadie, Diamond and Hainmueller (2015) to study the effect of German
    unification on per-capita GDP.

    As German reunification occurred in 1990, the `TreatmentPanel` for use in synthetic control
    estimation should be specified as:

    ```
    TreatmentPanel(df, "West Germany" => 1990,
        id_var = "country", t_var = "year", outcome_var = "gdp")
    ```
"""
function load_germany()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","germany.csv"), DataFrame; stringtype = String)
end

"""
    load_germany_panel()

    Returns a `BalancedPanel` object based on the data available from `load_germany()`, with
    treatment assignment specified as `"West Germany" => 1990`
"""
function load_germany_panel()
    df_germany = load_germany()
    BalancedPanel(df_germany, "West Germany" => 1990;
        id_var = "country", t_var = "year", outcome_var = "gdp")
end

"""
    load_basque()

    Returns a `DataFrame` with annual observations of per-capita GDP and covariates between 1955 and
    1997. The data is used in Abadie and Gardeazabal (2003) to study the effect of ETA terrorism
          on per-capita GDP in the Basque country.

    To undertake analysis simliar to that in Abadie and Gardeazabal (2003), the `TreatmentPanel` for
    use in synthetic control estimation should be specified as:

    ```
    TreatmentPanel(df, "Basque Country (Pais Vasco)" => 1970,
        id_var = "regionname", t_var = "year", outcome_var = "gdpcap")
    ```
"""
function load_basque()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","basque.csv"), DataFrame;
    stringtype = String, missingstring = "NA")
end

"""
    load_basque_panel()

    Returns a `BalancedPanel` object based on the data available from `load_basque()`, with
    treatment assignment specified as `"Basque Country (Pais Vasco)" => 1970`
"""
function load_basque_panel()
    df_basque = load_basque()
    BalancedPanel(df_basque, "Basque Country (Pais Vasco)" => 1970;
        id_var = "regionname", t_var = "year", outcome_var = "gdpcap")
end

"""
    load_smoking()

    Returns a `DataFrame` with annual observations of cigarette sales and covariates for 39 US
    states between 1970 and 2000. The data is used in Abadie, Diamond and Hainmueller (2010) to
    study the effect of Proposition 99 on cigarette sales in California

    To undertake analysis simliar to that in Abadie, Diamond and Hainmueller (2010), the
    `TreatmentPanel` for use in synthetic control estimation should be specified as:

    ```
    TreatmentPanel(df, 3 => 1989,
        id_var = "state", t_var = "year", outcome_var = "cigsale")
    ```
"""
function load_smoking()
    df_smoking = CSV.read(joinpath(dirname(@__FILE__),"..","data","smoking.csv"), DataFrame;
    stringtype = String, missingstring = "NA")
    select!(df_smoking, Not(:Column1))
end

"""
    load_smoking_panel()

    Returns a `BalancedPanel` object based on the data available from `load_smoking()`, with
    treatment assignment specified as `3 => 1989`
"""
function load_smoking_panel()
    df_smoking = load_smoking()
    BalancedPanel(df_smoking, 3 => 1989;
        id_var = :state, t_var = :year, outcome_var = :cigsale)
end

max_stars_repo_name: nilshg/SynthControl

max_issues_repo_name: nilshg/SynthControl

max_forks_repo_name: nilshg/SynthControl


text:
using PoincareInvariants
using GeometricIntegrators
using GeometricIntegrators.Utils
using Test


text:
@testset "Hubbard Model Attractive" begin
    # constructors
    m = HubbardModelAttractive(dims=1, L=8);
    @test m.L == 8 && m.dims == 1
    @test typeof(m) == MonteCarlo.HubbardModelAttractive{MonteCarlo.Chain}
    d = Dict(:dims=>2,:L=>3)
    m = HubbardModelAttractive(d)
    @test typeof(m) == MonteCarlo.HubbardModelAttractive{MonteCarlo.SquareLattice}
    @test m.L == 3 && m.dims == 2

end