text:
"""
    TemperatureSensor(; name)

Absolute temperature sensor in kelvin.

This is an ideal absolute temperature sensor which returns the temperature of the connected port in kelvin as an output
signal. The sensor itself has no thermal interaction with whatever it is connected to. Furthermore, no thermocouple-like
lags are associated with this sensor model.
"""
function TemperatureSensor(; name)
    @named port = HeatPort()
    @variables T(t) # [K] Absolute temperature
    eqs = [
        T ~ port.T
        port.Q_flow ~ 0
    ]
    ODESystem(eqs, t, [T], [], systems=[port], name=name)
end

"""
    RelativeTemperatureSensor(; name)

Relative Temperature sensor.

The relative temperature `port_a.T - port_b.T` is determined between the two ports of this component and is provided as
output signal in kelvin.
"""
function RelativeTemperatureSensor(; name)
    @named port_a = HeatPort()
    @named port_b = HeatPort()
    @variables T(t) # [K] Relative temperature a.T - b.T
    eqs = [
        T ~ port_a.T - port_b.T
        port_a.Q_flow ~ 0
        port_b.Q_flow ~ 0
    ]
    ODESystem(eqs, t, [T], [], systems=[port_a, port_b], name=name)
end

"""
    HeatFlowSensor(; name)

Heat flow rate sensor.

This model is capable of monitoring the heat flow rate flowing through this component. The sensed value of heat flow rate
is the amount that passes through this sensor while keeping the temperature drop across the sensor zero. This is an ideal
model so it does not absorb any energy and it has no direct effect on the thermal response of a system it is included in.
The output signal is positive, if the heat flows from `port_a` to `port_b`.
"""
function HeatFlowSensor(; name)
    @named port_a = HeatPort()
    @named port_b = HeatPort()
    @variables Q_flow(t) # [W] Heat flow from port a to port b
    eqs = [
        port_a.T ~ port_b.T
        port_a.Q_flow + port_b.Q_flow ~ 0
        Q_flow ~ port_a.Q_flow
    ]
    ODESystem(eqs, t, [Q_flow], [], systems=[port_a, port_b], name=name)
end


text:
"""
myEstimator

An abstract type for estimators.
"""
abstract type myEstimator end

# Methods ======================================================================
"""
coef(fit::myEstimator)

A method to retrieve the coefficient from a myEstimator object.
"""
function coef(fit::myEstimator)
    return fit.β
end #COEF.MYESTIMATOR

"""
predict(fit::myEstimator)

A method to calculate predictions of a myEstimator object.
"""
function predict(fit::myEstimator, data = nothing)
  # Calculate and return predictions
  isnothing(data) ? fitted = fit.X * fit.β : fitted = data * fit.β
  return(fitted)
end #PREDICT.MYESTIMATOR

ext: jl


text:
lsv = reval("ls()")
@test length(lsv) == 0
@test isa(lsv, RObject{StrSxp})

lsd = reval("ls(\"package:datasets\")")
@test isa(lsv, RObject{StrSxp})
@test length(lsd) > 50
@test "airmiles" in rcopy(lsd)

langsexp = rlang(:solve, RObject([1 2; 0 4]))
@test length(langsexp) == 2
@test rcopy(reval(langsexp)) == [1 -0.5; 0 0.25]
@test rcopy(langsexp[1]) == :solve
langsexp[1] = RObject(:det)
langsexp[2] = RObject([1 2; 0 0])
@test rcopy(reval(langsexp))[1] == 0

globalEnv[:x] = RObject([1,2,3])
@test rcopy(globalEnv[:x]) == [1,2,3]
globalEnv[:y] = RObject([4,5,6])
@test rcopy(rcall(Symbol("+"),:x,:y)) == [5,7,9]

x = 1:10
@rput x
@rget x
@test isa(x,Vector{Int32})
@test all(x .== 1:10)

y = "foo"
@rput x y::StrSxp
@rget x y::Array{Compat.UTF8String}
@test isa(y,Vector{Compat.UTF8String})
@test y[1] == "foo"

@test sprint(io -> rprint(io,RObject([1,2,3]))) == "[1] 1 2 3\n"

@test_approx_eq rcopy(rcall(:besselI, 1.0, 2.0)) besseli(2.0,1.0)
@test_approx_eq rcopy(rcall(:besselI, 1.0, 2.0, var"expon.scaled"=true)) besselix(2.0,1.0)

# callbacks
function testfn(x,y;a=3,b=4)
    [x;y;a;b]
end

r = rcall(testfn, 1, 2)
@test isa(r,RObject{IntSxp})
@test rcopy(r) == [1,2,3,4]

r = rcall(testfn, 1, 2,b=6)
@test isa(r,RObject{IntSxp})
@test rcopy(r) == [1,2,3,6]

r = rcall(:optimize,sin,[-2,0])
@test_approx_eq_eps r[:minimum][1] -pi/2 eps()^0.25
r = rcall(:optimize,sin,[0,2],maximum=true)
@test_approx_eq_eps r[:maximum][1] pi/2 eps()^0.25

nullfn() = nothing
@test isa(rcall(nullfn), RObject{NilSxp})

# graphics
RCall.rgui_init()
f = tempname()
rcall(:png,f)
rcall(:plot,1:10)
rcall(Symbol("dev.off"))
@test isfile(f)
@test !RCall.rgui_start(true)
@test_throws ErrorException RCall.rgui_start()
@test RCall.rgui_stop()

# S4 rprint
@test contains(sprint(io ->
rprint(io, reval("""
   setClass("Foo", representation(x = "numeric"))
   foo <- new("Foo", x = 20)
"""))), "An object of class")

# S3 rprint
@test contains(sprint(io ->
rprint(io, reval("""
   print.Bar <- function(x) print("hello")
   bar <- 1
   class(bar) <- "Bar"
   bar
"""))), "hello")

# operators
a = reval("a=c(1,2,3)")
b = reval("b=c(4,5,6)")
@test rcopy(a+b)==rcopy("a+b")
@test rcopy(a-b)==rcopy("a-b")
@test rcopy(a*b)==rcopy("a*b")
@test rcopy(a/b)==rcopy("a/b")
@test rcopy(a^b)==rcopy("a^b")


text:
## Set the number of cores to use
Num_Core = 25
if nprocs() != Num_Core
    addprocs(Num_Core - 1)
end
include("../../forc-helpers/prepare.jl")
include("../run_all.jl")
level = 4
for h in [3, 6, 12]
    for J in [301, 373, 493]
        h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2 = prepare(level, h, J)
        run_LAS2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    end
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
using MaximumGeneratorProfitTCRDD
using Documenter

makedocs(;
    modules=[MaximumGeneratorProfitTCRDD],
    authors="Invenia Technical Computing Corporation",
    repo="https://gitlab.invenia.ca/invenia/research/MaximumGeneratorProfitTCRDD.jl/blob/{commit}{path}#L{line}",
    sitename="MaximumGeneratorProfitTCRDD.jl",
    format=Documenter.HTML(;
        prettyurls=false,
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
        "Contributing" => "contributing.md",
    ],
    strict=true,
    checkdocs=:exports,
)

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# transfer_operators.jl : tranform quantities from fine to coarse grid and back
Base.broadcastable(scaling::UniformScaling) = Ref(scaling)

# restriction
R̃(dir::Int,op::TransferKind,n::Int...) = kron(sparse.(I,n[end:-1:dir+1].-1,n[end:-1:dir+1].-1)...,R₁(op,n[dir]),sparse.(I,n[dir-1:-1:1].-1,n[dir-1:-1:1].-1)...)

# interpolation
P̃(dir::Int,op::TransferKind,n::Int...) = kron(sparse.(I,n[end:-1:dir+1].-1,n[end:-1:dir+1].-1)...,P₁(op,n[dir]),sparse.(I,n[dir-1:-1:1].-1,n[dir-1:-1:1].-1)...)


text:
using TimedAutomata, Base.Iterators
using GraphMakie
##

##
neuron_sequential = Neuron((((((),:A),),:B),),:C);
neuron_sequential_ta = TA(neuron_sequential)
neuron_sequential_p = prune(neuron_sequential_ta)

#zone = TAConstraint(neuron_sequential_p.clocks, zeros(TABound{Int},length(neuron_sequential_p.clocks)+1,length(neuron_sequential_p.clocks)+1))
zone = @constraint clk_1-clk_2==0 && clk_2-clk_3==0 && clk_3-clk_1==0 Int
TimedAutomata.set_clock_list!(zone, neuron_sequential_ta.clocks)
neuron_sequential_tts = TTS(neuron_sequential_p; initial_zone=zone)
@assert all(neuron_sequential_tts.states[end][1].components[1])
sequential_language_full=language(neuron_sequential_tts; final_states=neuron_sequential_tts.states[end:end], ignore_messages=(_)->false)
sequential_language_no_silent=language(neuron_sequential_tts; final_states=neuron_sequential_tts.states[end:end])
sequential_language_input=language(neuron_sequential_tts; final_states=neuron_sequential_tts.states[end:end], forbid_messages=(m)->m.symbol∈(:↻,), ignore_messages=(m)->m.direction!=MSG_IN)
sequential_language_no_reset_no_spike=language(neuron_sequential_tts; final_states=neuron_sequential_tts.states[end:end], forbid_messages=(m)->m.symbol∈(:spike,:↻))

##
neuron_parallel = Neuron((((),:A),((),:B)),:C);
neuron_parallel_ta = TA(neuron_parallel)
neuron_parallel_p = prune(neuron_parallel_ta)
neuron_parallel_tts = TTS(neuron_parallel_p; initial_zone=zone)
id = findfirst(((l,_),)->all(l.components[1]), neuron_parallel_tts.states)
parallel_language_full=language(neuron_parallel_tts; final_states=neuron_parallel_tts.states[id:id], ignore_messages=(_)->false)
parallel_language_no_silent=language(neuron_parallel_tts; final_states=neuron_parallel_tts.states[id:id])
parallel_language_input=language(neuron_parallel_tts; final_states=neuron_parallel_tts.states[id:id], forbid_messages=(m)->m.symbol∈(:↻,), ignore_messages=(m)->m.direction!=MSG_IN)
parallel_language_no_reset_no_spike=language(neuron_parallel_tts; final_states=neuron_parallel_tts.states[id:id], forbid_messages=(m)->m.symbol∈(:spike,:↻))

## Plotting

f,ax,_ = graphplot(neuron_sequential_ta);
hidedecorations!(ax)
display(f)
print("Press key to continue: ")
read(stdin, 1)
save(joinpath("examples","neurons", "figures","sequential_unpruned.png"), f)

f,ax,_ = graphplot(neuron_sequential_p);
hidedecorations!(ax)
display(f)
print("Press key to continue: ")
read(stdin, 1)
save(joinpath("examples","neurons", "figures","sequential_pruned.png"), f)

f,ax,_ = graphplot(neuron_parallel_p);
hidedecorations!(ax)
display(f)
print("Press key to continue: ")
read(stdin, 1)
save(joinpath("examples","neurons", "figures","parallel_pruned.png"), f)


text:
type DirichiletLift

    N::Int
    nd::Int
    nv::Int
    idx::Array{Int,1}
    vidx::Array{Int,1}
    A::Array{Float64,2}
    function DirichiletLift(n, is_dirichilet)
        nd = 0
        for i = 1:n
            nd = is_dirichilet[i]?nd+1:nd
        end
        nv = n - nd
        idx = zeros(Int, nd)
        vidx = zeros(Int, n-nd)
        icnt = 1
        ivar = 1
        for i = 1:n
            if is_dirichilet[i]
                idx[icnt] = i
                icnt += 1
            else
                vidx[ivar] = i
                ivar += 1
            end
        end

        A = zeros(nv, nd)

        new(n, nd, nv, idx, vidx, A)
    end

end

function set_dirichilet_matrix(dlift::DirichiletLift, M)
    A = dlift.A
    N = dlift.N
    nd = dlift.nd
    nv = dlift.nv
    idx = dlift.idx
    vidx = dlift.vidx

    for i = 1:nv
        iv = vidx[i]
        for k = 1:nd
            A[i,k] = M[iv, idx[k]]
        end
    end

end

function lift(dlift::DirichiletLift, f)

    idx = dlift.idx
    nv = dlift.nv
    nd = dlift.nd
    vidx = dlift.vidx
    A = dlift.A

    for i = 1:nv
        iv = vidx[i]
        for k = 1:nd
            f[iv] -= A[i,k] * f[idx[k]]
        end
    end
end

function lift(dlift::DirichiletLift, f, u)

    idx = dlift.idx
    nv = dlift.nv
    nd = dlift.nd
    vidx = dlift.vidx
    A = dlift.A

    for i = 1:nv
        iv = vidx[i]
        for k = 1:nd
            f[iv] -= A[i,k] * u[idx[k]]
        end
    end
end


text:
# setting up needed Julia packages
using Pkg
Pkg.add(["HTTP","JSON2"])
using HTTP, JSON2

# define NOAA NDBC data access web address & data access token
base_ncdc_url = "https://www.ncdc.noaa.gov/cdo-web/api/v2/";
token = "FMAPxDJFXSsMKsJWUDiOhZkINCfYGZqX"; # need to request your own from https://www.ncdc.noaa.gov/cdo-web/token

# use HTTP.get to obtain location categories, see https://www.ncdc.noaa.gov/cdo-web/webservices/v2#gettingStarted
endpoint = "locationcategories";
r = HTTP.get(base_ncdc_url * endpoint, ["token" => token])
@pretty String(r.body)

endpoint = "datasets";
id = "/GHCND"
rds = HTTP.get(base_ncdc_url * endpoint * id, ["token" => token])
@pretty String(rds.body)

endpoint = "data";
datasetid = "GHCND"
datatypeid = "TOBS"
stationid = Dict(["williamsburg" => "GHCND:USC00449151", "newportnews" => "GHCND:USW00093741"]) # Williamsburg, Newport News
startdate = "2020-04-01"
enddate = "2020-04-03"
rd = HTTP.get(base_ncdc_url * endpoint * "?datasetid=" * datasetid * "&stationid=" * stationid["williamsburg"] * "&datatypeid=" * datatypeid * "&startdate=" * startdate * "&enddate=" * enddate, ["token" => token])
@pretty String(rd.body)


text:
# TODO probably deserves separate struct, or deeper looking into internet
function mixed_radix_conv(indexes::Vector{Int}, dims::Vector{Int})
    @assert length(indexes) == length(dims)
    ind = 1 # because Julia enumerates from 1
    d_prod = 1
    for (i, d) = reverse(collect(zip(indexes .- 1, dims)))
        ind += d_prod*i
        d_prod *= d
    end
    ind
end

function mixed_radix_conv(index::Int, dims::Vector{Int})
    @assert 1 <= index <= prod(dims)
    index -= 1
    n = length(dims)
    res = ones(Int, n)
    for (ind, d) = enumerate(reverse(dims))
        res[ind] += index % d
        index = div(index, d)
    end
    reverse(res)
end


text:
using DataFrames
import SMC

@everywhere struct Problem
    x::Vector{Float64}
    y::Vector{Float64}
    prior_std::Float64
    noise_std::Float64
end

@everywhere begin
    include("../aide.jl")
    include("exact_closed_form.jl")
    include("smc.jl")
    include("bbvi.jl")
    include("mcmc.jl")
end

include("schema.jl")

function time_sampler(sampler::Sampler)
    const NUM_TIMING_RUNS = 100
    times = []
    for i=1:NUM_TIMING_RUNS
        tic()
        simulate(sampler)
        push!(times, toq())
    end
    median(times)
end

function generate_prior_data(problem::Problem, num_replicates::Int)
    println("generating prior data...")
    exact_sampler = ExactLinregSampler(problem)
    prior_sampler = PriorSampler(problem.x, problem.y, problem.prior_std, problem.noise_std)
    kls = aide(exact_sampler, prior_sampler, 1, 1, num_replicates, num_replicates)

    data = DataFrame()
    data[COL_SAMPLER_NAME] = ["prior"]
    data[COL_AIDE_ESTIMATE] = [mean(kls)]
    data[COL_AIDE_STDERR] = [std(kls) / sqrt(num_replicates)]
    data[COL_MEDIAN_RUNTIME] = [time_sampler(prior_sampler)]
    return data
end

function generate_posterior_data(problem::Problem, num_replicates::Int)
    println("generating posterior data...")
    exact_sampler = ExactLinregSampler(problem)
    kls = aide(exact_sampler, exact_sampler, 1, 1, num_replicates, num_replicates)

    data = DataFrame()
    data[COL_SAMPLER_NAME] = ["posterior"]
    data[COL_AIDE_ESTIMATE] = [mean(kls)]
    data[COL_AIDE_STDERR] = [std(kls) / sqrt(num_replicates)]
    data[COL_MEDIAN_RUNTIME] = [time_sampler(exact_sampler)]
    return data
end

function generate_smc_data(problem::Problem, num_replicates::Int)
    println("generating SMC data...")
    exact_sampler = ExactLinregSampler(problem)

    k_column = Int[]
    num_rejuvenation_column = Int[]
    num_particles_column = Int[]
    aide_est_column = Float64[]
    aide_stderr_column = Float64[]
    median_runtime_column = Float64[]

    all_ks = [1, 10, 10^3]
    all_num_rejuvenation = [2] #[0, 1, 2, 4, 8] # TODO add back to make the rejuv. plot for appendix
    all_num_particles = [1, 2, 3, 6, 10, 20, 30, 60, 100]
    for k in all_ks
        for num_rejuvenation in all_num_rejuvenation
            for num_particles in all_num_particles
                println("k=$k, num_rejuvenation=$num_rejuvenation, num_particles=$num_particles")
                push!(k_column, k)
                push!(num_rejuvenation_column, num_rejuvenation)
                push!(num_particles_column, num_particles)
                smc = LinregSMCSampler(problem.x, problem.y, problem.prior_std, problem.noise_std, num_rejuvenation, num_particles)

                kls = aide(exact_sampler, smc, 1, k, num_replicates, num_replicates)
                push!(aide_est_column, mean(kls))
                push!(aide_stderr_column, std(kls) / sqrt(num_replicates))
                push!(median_runtime_column, time_sampler(smc))
            end
        end
    end

    data = DataFrame()
    data[COL_SAMPLER_NAME] = fill("smc", length(aide_est_column))
    data[COL_NUM_METAINFERENCE] = k_column
    data[COL_NUM_REJUVENATION] = num_rejuvenation_column
    data[COL_NUM_PARTICLES] = num_particles_column
    data[COL_AIDE_ESTIMATE] = aide_est_column
    data[COL_AIDE_STDERR] = aide_stderr_column
    data[COL_MEDIAN_RUNTIME] = median_runtime_column
    return data
end

function generate_bbvi_data(problem::Problem, num_replicates::Int)
    println("generating BBVI data...")
    exact_sampler = ExactLinregSampler(problem)

    all_num_vi_iters = convert(Array{Int,1},floor.(10.^[0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]))
    bbvi_num_samples = 10
    bbvi_a = 10000.0
    bbvi_b = 0.75

    num_sgd_iters_column = Int[]
    aide_est_column = Float64[]
    aide_stderr_column = Float64[]
    median_runtime_column = Float64[]
    for num_iters in all_num_vi_iters
        println("num_iters=$num_iters")
        push!(num_sgd_iters_column, num_iters)
        params = BBVIParams(0.9, bbvi_num_samples, bbvi_a, bbvi_b, num_iters, problem.prior_std, problem.noise_std) # 1000 # was 0.7
        vi_sampler = BBVISampler(params, problem.x, problem.y)
        kls = aide(exact_sampler, vi_sampler, 1, 1, num_replicates, num_replicates)
        push!(aide_est_column, mean(kls))
        push!(aide_stderr_column, std(kls) / sqrt(num_replicates))
        times = []
        for i=1:num_replicates
            tic()
            simulate(vi_sampler)
            push!(times, toq())
        end
        push!(median_runtime_column, median(times))
    end

    data = DataFrame()
    data[COL_SAMPLER_NAME] = fill("smc", length(aide_est_column))
    data[COL_NUM_METAINFERENCE] = fill(1, length(aide_est_column))
    data[COL_NUM_SGD_ITERS] = num_sgd_iters_column
    data[COL_AIDE_ESTIMATE] = aide_est_column
    data[COL_AIDE_STDERR] = aide_stderr_column
    data[COL_MEDIAN_RUNTIME] = median_runtime_column
    return data
end

function generate_mcmc_data(problem::Problem, num_replicates::Int)
    println("generating MCMC data...")
    exact_sampler = ExactLinregSampler(problem)

    k_column = Int[]
    num_mcmc_iters_column = Int[]
    aide_est_column = Float64[]
    aide_stderr_column = Float64[]
    median_runtime_column = Float64[]

    all_num_mh_iters = [0, 1, 2, 4, 8, 16, 32, 64, 128]
    all_ks = [1, 10, 10^3]
    for num_iters in all_num_mh_iters
        mh_sampler = LinregMHSampler(problem.x, problem.y, problem.prior_std, problem.noise_std, num_iters)
        for k in all_ks
            println("num_iters=$num_iters, k=$k")
            push!(k_column, k)
            push!(num_mcmc_iters_column, num_iters)
            kls = aide(exact_sampler, mh_sampler, 1, k, num_replicates, num_replicates)
            push!(aide_est_column, mean(kls))
            push!(aide_stderr_column, std(kls) / sqrt(num_replicates))
            push!(median_runtime_column, time_sampler(mh_sampler))
        end
    end

    data = DataFrame()
    data[COL_SAMPLER_NAME] = fill("mcmc", length(aide_est_column))
    data[COL_NUM_METAINFERENCE] = k_column
    data[COL_NUM_MCMC_ITERS] = num_mcmc_iters_column
    data[COL_AIDE_ESTIMATE] = aide_est_column
    data[COL_AIDE_STDERR] = aide_stderr_column
    data[COL_MEDIAN_RUNTIME] = median_runtime_column
    return data
end

problem = Problem([-2.0, 2.0], [-2.0, 2.0], 1.0, 1.0)
num_replicates =  10000
data_dir = "data"
if !Base.Filesystem.isdir(data_dir)
    Base.Filesystem.mkdir(data_dir)
end

prior_data = generate_prior_data(problem, num_replicates)
writetable("$data_dir/prior_data.csv", prior_data)

posterior_data = generate_posterior_data(problem, num_replicates)
writetable("$data_dir/posterior_data.csv", posterior_data)

smc_data = generate_smc_data(problem, num_replicates)
writetable("$data_dir/smc_data.csv", smc_data)

bbvi_data = generate_bbvi_data(problem, num_replicates)
writetable("$data_dir/bbvi_data.csv", bbvi_data)

mcmc_data = generate_mcmc_data(problem, num_replicates)
writetable("$data_dir/mcmc_data.csv", mcmc_data)


text:
using MPI
using MPIAsyncPools
using Test

function parse_line(line::AbstractString; prefix="test")
    if !startswith(line, prefix)
        return -1, ""
    end
    line = line[length(prefix)+1:end]
    segments = split(line, ":")
    length(segments) == 2 || error("Unexpected colon in $line")
    rank = parse(Int, segments[1])
    output = strip(segments[2])
    rank, output
end

mpitestexec(n::Integer, file::String) = mpiexec(cmd -> read(`$cmd -n $n --mca btl ^openib julia --project $file`, String))

@testset "MPIAsyncPools.jl" begin
    n = 3
    for line in split(mpitestexec(n, "kmap1.jl"), "\n")
        rank, output = parse_line(line)
        if rank == -1
            continue
        end
        println((rank, line))
    end

    n = 3
    for line in split(mpitestexec(n, "kmap2.jl"), "\n")
        rank, output = parse_line(line)
        if rank == -1
            continue
        end
        println((rank, line))
    end

    n = 10
    for line in split(mpitestexec(n, "kmap2.jl"), "\n")
        rank, output = parse_line(line)
        if rank == -1
            continue
        end
        println((rank, line))
    end

    println("No error messages indicates that all tests passed")
end


text:
@testset "pendulum" begin

    env = PendulumEnv(; rng = StableRNG(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    env = PendulumNonInteractiveEnv(; rng = StableRNG(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
using DiffEqBase, ClimaTimeSteppers, LinearAlgebra, Test

include("utils.jl")

dts = 0.5 .^ (4:7)

for (prob, sol, tscale) in [
    (linear_prob, linear_sol, 1)
    (sincos_prob, sincos_sol, 1)
]

    @test convergence_order(prob, sol, LSRKEulerMethod(), dts.*tscale)     ≈ 1 atol=0.1
    @test convergence_order(prob, sol, LSRK54CarpenterKennedy(), dts.*tscale)       ≈ 4 atol=0.05
    @test convergence_order(prob, sol, LSRK144NiegemannDiehlBusch(),  dts.*tscale)   ≈ 4 atol=0.05

    @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
    @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
    @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
    @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05
end

for (prob, sol, tscale) in [
    (linear_prob_wfactt, linear_sol, 1)
]
    @test convergence_order(prob, sol, SSPKnoth(linsolve=linsolve_direct), dts.*tscale) ≈ 2 atol=0.05

end

# ForwardEulerODEFunction
for (prob, sol, tscale) in [
    (linear_prob_fe, linear_sol, 1)
    (sincos_prob_fe, sincos_sol, 1)
]
    @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
    @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
    @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
    @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05

end

if ArrayType == Array
for (prob, sol) in [
    imex_autonomous_prob => imex_autonomous_sol,
    #imex_nonautonomous_prob => imex_nonautonomous_sol,
]
    # IMEX
    @test convergence_order(prob, sol, ARK1ForwardBackwardEuler(linsolve=DirectSolver), dts)       ≈ 1 atol=0.1
    @test convergence_order(prob, sol, ARK2ImplicitExplicitMidpoint(linsolve=DirectSolver), dts)   ≈ 2 atol=0.05
    @test convergence_order(prob, sol, ARK2GiraldoKellyConstantinescu(linsolve=DirectSolver), dts) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, ARK2GiraldoKellyConstantinescu(linsolve=DirectSolver; paperversion=true), dts) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, ARK437L2SA1KennedyCarpenter(linsolve=DirectSolver), dts)    ≈ 4 atol=0.05
    @test convergence_order(prob, sol, ARK548L2SA2KennedyCarpenter(linsolve=DirectSolver), dts)    ≈ 5 atol=0.05
end
end
for (prob, sol) in [
    imex_autonomous_prob => imex_autonomous_sol,
    imex_nonautonomous_prob => imex_nonautonomous_sol,
    # kpr_multirate_prob => kpr_sol,
]
    # Multirate
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),LSRK54CarpenterKennedy()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 4 atol=0.05
    # MIS
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),MIS2()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),MIS3C()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),MIS4()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 3 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),MIS4a()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 3 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),TVDMISA()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),TVDMISB()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05

    # Wicker Skamarock
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),WSRK2()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05
    @test convergence_order(prob, sol, Multirate(LSRK54CarpenterKennedy(),WSRK3()), dts;
        fast_dt = 0.5^12, adjustfinal=true) ≈ 2 atol=0.05

    end


text:
function build(type::Type{VLAgentModel}; kwargs...)

    try

        # get the args as a dictionary -
        args = Dict(kwargs)

        # setup default values for the key words args -
        agent_trade_logic = get(args, :agent_trade_logic, _default_agent_trade_logic)
        agent_update_logic = get(args, :agent_update_logic, _default_agent_update_logic!)
        number_of_iterations = get(args, :number_of_iterations, 100)
        number_of_assets = get(args, :number_of_assets, 1)
        wealth = get(args, :agent_wealth, 1000.0)
        initial_position_array = get(args, :agent_position_array, zeros(number_of_assets))

        # build a blank agent -
        agent_model = VLAgentModel()
        agent_model.agent_trade_logic = agent_trade_logic
        agent_model.agent_update_logic = agent_update_logic
        agent_model.agent_id = uuid4()

        # what is the agents initial wealth -
        wealth_array = Array{Float64,1}()
        push!(wealth_array, wealth)
        agent_model.agent_wealth_array = wealth_array

        # initialize the initial position array -
        agent_position_array = zeros(number_of_iterations, number_of_assets)
        for asset_index = 1:number_of_assets
            agent_position_array[1,asset_index] = initial_position_array[asset_index]
        end
        agent_model.agent_position_array = agent_position_array

        # return -
        return agent_model
    catch error
        rethrow(error)
    end
end

function build(type::Type{VLGameWorld}; kwargs...)

    try

        # get the args as a dictionary -
        args = Dict(kwargs)

        # check for defaults -
        number_of_iterations = get(args, :number_of_iterations, 100)
        number_of_assets = get(args, :number_of_assets, 1)
        exchange_logic = get(args, :market_exchange_logic, _default_exchange_logic!)
        liquidity_parameter_array = get(args, :liquidity_parameter_array, 0.01 * ones(number_of_assets))

        # Create world array -
        game_world = VLGameWorld()
        game_world.number_of_iterations = number_of_iterations
        game_world.exchange_logic = exchange_logic
        game_world.liquidity_parameter_array = liquidity_parameter_array
        game_world.asset_price_array = ones(number_of_iterations, number_of_assets)

        # return -
        return game_world
    catch error
        rethrow(error)
    end
end


text:
using Base.Cartesian

"""
    @reduce(n::Int, op::Symbol, ex::Expr)

        Wrapping macro for @ncall. Is equivalent and generates
        op(ex_1, ..., ex_n)

    Example:

        @reduce 2 (+) u -> v_u^2
        > Generates +(v_1^2, v_2^2)
"""
macro reduce(n::Int, op::Symbol, ex::Expr)
    esc(:(@ncall($n, $op, $ex)))
end

"""
    VerletList(size::Int, capacity::Int, cutoff::Float64, offset::Vector{Int}, list::Vector{Int})
    VerletList(size::Int)

        Hold information regarding the neighbouring particles of each particle
        in the system.
"""
mutable struct VerletList
    size::Int              # number of particles the list refers to
    capacity::Int          # capacity of list
    cutoff::Float64        # cutoff (nm)
    offset::Vector{Int}    # pointer particle->list position
    list::Vector{Int}      # nblist
end
VerletList(size::Int) = VerletList(size, size, -1.0, zeros(Int, size), zeros(Int, size))

"""
    resize!(vlist::VerletList, n::Int)

        Resize a Verlet list to hold 'n' particles. If 'n' is higher than the
        current length of the list, appends 0's.

    Example:

        resize!(vl, 8)
"""
Base.resize!(vlist::VerletList, n::Int) = begin
    resize!(vlist.list, n)
    vlist.capacity = n
    vlist
end

"""
    get_coords(coords::Matrix{T}, u::Int, i::Int) where {T <: AbstractFloat}

        Returns, from the coordinates at array position (atom number) 'i', the
        dimendion number 'u' (i.e. dimension X is 1, Y is 2 and Z is 3). Allows
        for correct coordinates retrieval independent of AoS vs SoA format.

    Example:

        get_coords(coords, u, i)
"""
function get_coords(coords::Matrix{T}, u::Int, i::Int) where {T <: AbstractFloat}
    if argmin(collect(size(coords))) == 1
        return coords[u, i]
    else
        return coords[i, u]
    end
end

"""
    update_serial!(vlist::VerletList, coords::Matrix{T}) where {T <: AbstractFloat}

        Updates the Verlet list (using a serial approach) according to the
        defined 'vlist.cutoff' and the given coordinates 'coords' (can be in
        both AoS or SoA formats).

    Example:

        update_serial!(vl, xyz)
"""
function update_serial!(vlist::VerletList, coords::Matrix{T}) where {T<:AbstractFloat}

    @assert vlist.size == max(size(coords)...) "incompatible sizes"
    # println(size(coords))

    # cutoff squared
    cutsq = convert(T, vlist.cutoff*vlist.cutoff)

    offset = 1
    natoms = vlist.size
    @inbounds for i = 1:natoms

        vlist.offset[i] = offset
        # @nexprs 3 u -> xi_u = coords[i,u]
        @nexprs 3 u -> xi_u = coords[i, u]

        for j = (i+1):natoms

            @nexprs 3 u -> vij_u = coords[j, u] - xi_u
            # @nexprs 3 u -> vij_u = coords[j,u] - xi_u
            dij_sq = @reduce 3 (+) u -> vij_u*vij_u

            if dij_sq < cutsq
                vlist.list[offset] = j
                offset += 1
                if offset == vlist.capacity
                    resize!(vlist, vlist.capacity + natoms)
                end
            end
        end

        vlist.list[offset] = -1
        offset += 1
        if (i < natoms) && (offset == vlist.capacity)
            resize!(vlist, vlist.capacity + natoms)
        end

    end

    vlist
end

# TODO: update_simd4!
# TODO: update_simd8!


text:
using MCMCChains
using StatsBase
using Test
using Random

# Tests for missing values.
function testdiff(cdf1, cdf2)
    m1 = convert(Array, cdf1)
    m2 = convert(Array, cdf2)
    return all(((x, y),) -> isapprox(x, y; atol=1e-2), zip(m1, m2))
end

@testset "utils" begin
    x = [1, missing, 3, 2]
    @test MCMCChains.cummean(x) == [1., 1., 2., 2.]
end

@testset "stats" begin
    Random.seed!(1234)
    chn = Chains(randn(1000, 2, 2))

    # Call describe without missing values.
    describe(devnull, chn; showall=true)
    s1, s2 = summarystats(chn), quantile(chn)

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, 2) .* missing, dims = 1))

    # Call describe with missing values.
    describe(devnull, chn_m; showall=true)
    m1, m2 = summarystats(chn_m), quantile(chn_m)

    @test testdiff(s1, m1)
end

@testset "diagnostic functions" begin
    Random.seed!(1234)
    nchains = 2
    chn = Chains(randn(5000, 2, nchains))

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, nchains) .* missing, dims = 1))

    # Currently we only check if diag. functions throw an assertion if a value is missing.
    @test_throws AssertionError gelmandiag(chn_m)

    gw_1 = gewekediag(chn)
    gw_2 = gewekediag(chn_m)
    hd_1 = heideldiag(chn)
    hd_2 = heideldiag(chn_m)
    rf_1 = rafterydiag(chn)
    rf_2 = rafterydiag(chn_m)

    @testset "diagnostics missing tests" for i in 1:nchains
        @test testdiff(gw_1, gw_2)
        @test testdiff(hd_1, hd_2)
        @test testdiff(rf_1, rf_2)
    end

    @test_throws AssertionError discretediag(chn_m)
end


text:
using ComputationalResources
using Base.Test

for r in (CPU1(5), CPUThreads(5), ArrayFireLibs(5), CUDALibs(5), OpenCLLibs(5))
    for T in (CPU1, CPUThreads, ArrayFireLibs, CUDALibs, OpenCLLibs)
        r1 = T(r)
        @test r1.settings === 5
    end
end

push!(LOAD_PATH, joinpath(dirname(@__FILE__), "packages"))

import Dummy1

@test Dummy1.foo(CPU1())
@test_throws MethodError Dummy1.foo(ArrayFireLibs())

addresource(ArrayFireLibs)
@test haveresource(ArrayFireLibs)

reload("Dummy1")

@test Dummy1.foo(CPU1())
@test Dummy1.foo(ArrayFireLibs())

rmresource(ArrayFireLibs)

reload("Dummy1")

@test Dummy1.foo(CPU1())
@test_throws MethodError Dummy1.foo(ArrayFireLibs())

pop!(LOAD_PATH)

@test isa(CPUThreads(), AbstractCPU{Void})
@test isa(CUDALibs(), AbstractResource{Void})
@test isa(OpenCLLibs(), AbstractResource{Void})

nothing


text:
import SDtoolbox: CJspeed, CANTERA_MECH_FILETYPE

P₁ = 1e5
T₁ = 300
ϕ = 1

X₁ = "H2:$(42*ϕ), O2:21,N2:79"
mech = "gri30."*CANTERA_MECH_FILETYPE

# reference value from online NASA CEA
V_CJ_ref = 1964.6

# compare results for all CJ algorithms
algorithms = "CEA", "umin", "aeq"
V_CJ = Float64[]

for method in algorithms
    push!(V_CJ,CJspeed(P₁,T₁,X₁,mech,method=method))
end

rel_err = maximum(abs.(1 .- V_CJ / (sum(V_CJ)/3) ))
ref_err = maximum(abs.(V_CJ_ref .- V_CJ))/V_CJ_ref


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

#This script simulates the Fisher-KPP equation and fits
#a neural PDE to the data with the growth (aka reaction) term replaced
#by a feed-forward neural network and the diffusion term with a CNN

using PyPlot, Printf
using LinearAlgebra
using Flux, DiffEqFlux, Optim, DiffEqSensitivity
using BSON: @save, @load
using Flux: @epochs
using OrdinaryDiffEq

#parameter
D = 0.01; #diffusion
r = 1.0; #reaction rate

#domain
X = 1.0; T = 5;
dx = 0.04; dt = T/10;
x = collect(0:dx:X);
t = collect(0:dt:T);
Nx = Int64(X/dx+1);
Nt = Int64(T/dt+1);

#initial conditions
Amp = 1.0;
Delta = 0.2
#IC-1
rho0 = Amp*(tanh.((x .- (0.5 - Delta/2))/(Delta/10)) - tanh.((x .- (0.5 + Delta/2))/(Delta/10)))/2
#IC-2
#rho0 = Amp*(1 .- tanh.((x .- 0.2)/(Delta/6)))/2.

save_folder = "data"

if isdir(save_folder)
    rm(save_folder, recursive=true)
end
mkdir(save_folder)

close("all")
figure()
plot(x, rho0)
title("Initial Condition")
gcf()

########################
# Generate training data
########################
reaction(u) = r * u .* (1 .- u)
lap = diagm(0 => -2.0 * ones(Nx), 1=> ones(Nx-1), -1 => ones(Nx-1)) ./ dx^2
#Periodic BC
lap[1,end] = 1.0/dx^2
lap[end,1] = 1.0/dx^2
#Neumann BC
#lap[1,2] = 2.0/dx^2
#lap[end,end-1] = 2.0/dx^2

function rc_ode(rho, p, t)
    #finite difference
    D * lap * rho + reaction.(rho)
end

prob = ODEProblem(rc_ode, rho0, (0.0, T), saveat=dt)
sol = solve(prob, Tsit5());
ode_data = Array(sol);

figure(figsize=(8,3))

subplot(121)
pcolor(x,t,ode_data')
xlabel("x"); ylabel("t");
colorbar()

subplot(122)
for i in 1:2:Nt
    plot(x, ode_data[:,i], label="t=$(sol.t[i])")
end
xlabel("x"); ylabel(L"$\rho$")
legend(frameon=false, fontsize=7, bbox_to_anchor=(1, 1), loc="upper left", ncol=1)
tight_layout()
savefig(@sprintf("%s/training_data.pdf", save_folder))
gcf()

########################
# Define the neural PDE
########################
n_weights = 10

#for the reaction term
rx_nn = Chain(Dense(1, n_weights, tanh),
                Dense(n_weights, 2*n_weights, tanh),
                Dense(2*n_weights, n_weights, tanh),
                Dense(n_weights, 1),
                x -> x[1])

#conv with bias with initial values as 1/dx^2
w_err = 0.0
init_w = reshape([1.1 -2.5 1.0], (3, 1, 1, 1))
diff_cnn_ = Conv(init_w, [0.], pad=(0,0,0,0))

#initialize D0 close to D/dx^2
D0 = [6.5]

p1,re1 = Flux.destructure(rx_nn)
p2,re2 = Flux.destructure(diff_cnn_)
p = [p1;p2;D0]
full_restructure(p) = re1(p[1:length(p1)]), re2(p[(length(p1)+1):end-1]), p[end]

function nn_ode(u,p,t)
    rx_nn = re1(p[1:length(p1)])

    u_cnn_1   = [p[end-4] * u[end] + p[end-3] * u[1] + p[end-2] * u[2]]
    u_cnn     = [p[end-4] * u[i-1] + p[end-3] * u[i] + p[end-2] * u[i+1] for i in 2:Nx-1]
    u_cnn_end = [p[end-4] * u[end-1] + p[end-3] * u[end] + p[end-2] * u[1]]

    # Equivalent using Flux, but slower!
    #CNN term with periodic BC
    #diff_cnn_ = Conv(reshape(p[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
    #u_cnn = reshape(diff_cnn_(reshape(u, (Nx, 1, 1, 1))), (Nx-2,))
    #u_cnn_1 = reshape(diff_cnn_(reshape(vcat(u[end:end], u[1:1], u[2:2]), (3, 1, 1, 1))), (1,))
    #u_cnn_end = reshape(diff_cnn_(reshape(vcat(u[end-1:end-1], u[end:end], u[1:1]), (3, 1, 1, 1))), (1,))

    [rx_nn([u[i]])[1] for i in 1:Nx] + p[end] * vcat(u_cnn_1, u_cnn, u_cnn_end)
end

########################
# Soving the neural PDE and setting up loss function
########################
prob_nn = ODEProblem(nn_ode, rho0, (0.0, T), p)
sol_nn = concrete_solve(prob_nn,Tsit5(), rho0, p)

function predict_rd(θ)
  # No ReverseDiff if using Flux
  Array(concrete_solve(prob_nn,Tsit5(),rho0,θ,saveat=dt,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

#match data and force the weights of the CNN to add up to zero
function loss_rd(p)
    pred = predict_rd(p)
    sum(abs2, ode_data .- pred) + 10^2 * abs(sum(p[end-4 : end-2])), pred
end

########################
# Training
########################

#Optimizer
opt = ADAM(0.001)

global count = 0
global save_count = 0
save_freq = 50

train_arr = Float64[]
diff_arr = Float64[]
w1_arr = Float64[]
w2_arr = Float64[]
w3_arr = Float64[]

#callback function to observe training
cb = function (p,l,pred)
    rx_nn, diff_cnn_, D0 = full_restructure(p)
    push!(train_arr, l)
    push!(diff_arr, p[end])

    weight = diff_cnn_.weight[:]
    push!(w1_arr, weight[1])
    push!(w2_arr, weight[2])
    push!(w3_arr, weight[3])

    println(@sprintf("Loss: %0.4f\tD0: %0.4f Weights:(%0.4f,\t %0.4f, \t%0.4f) \t Sum: %0.4f"
            ,l, D0[1], weight[1], weight[2], weight[3], sum(weight)))

    global count

    if count==0
        fig = figure(figsize=(8,2.5));
        ttl = fig.suptitle(@sprintf("Epoch = %d", count), y=1.05)
        global ttl
        subplot(131)
        pcolormesh(x,t,ode_data')
        xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
        colorbar()

        subplot(132)
        img = pcolormesh(x,t,pred')
        global img
        xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
        colorbar(); clim([0, 1]);

        ax = subplot(133); global ax
        u = collect(0:0.01:1)
        rx_line = plot(u, rx_nn.([[elem] for elem in u]), label="NN")[1];
        global rx_line
        plot(u, reaction.(u), label="True")
        title("Reaction Term")
        legend(loc="upper right", frameon=false, fontsize=8);
        ylim([0, r*0.25+0.2])

        subplots_adjust(top=0.8)
        tight_layout()
    end

    if count>0
        println("updating figure")
        img.set_array(pred[1:end-1, 1:end-1][:])
        ttl.set_text(@sprintf("Epoch = %d", count))

        u = collect(0:0.01:1)
        rx_pred = rx_nn.([[elem] for elem in u])
        rx_line.set_ydata(rx_pred)
        u = collect(0:0.01:1)

        min_lim = min(minimum(rx_pred), minimum(reaction.(u)))-0.1
        max_lim = max(maximum(rx_pred), maximum(reaction.(u)))+0.1

        ax.set_ylim([min_lim, max_lim])
    end

    global save_count
    if count%save_freq == 0
        println("saved figure")
        savefig(@sprintf("%s/pred_%05d.png", save_folder, save_count), dpi=200, bbox_inches="tight")
        save_count += 1
    end

    display(gcf())
    count += 1

    false
end

#train
res1 = DiffEqFlux.sciml_train(loss_rd, p, ADAM(0.001), cb=cb, maxiters = 100)
res2 = DiffEqFlux.sciml_train(loss_rd, res1.minimizer, ADAM(0.001), cb=cb, maxiters = 300)
res3 = DiffEqFlux.sciml_train(loss_rd, res2.minimizer, BFGS(), cb=cb, maxiters = 1000)

pstar = res3.minimizer

## Save trained model
@save @sprintf("%s/model.bson", save_folder) pstar

########################
# Plot for paper
########################
@load @sprintf("%s/model.bson", save_folder) pstar
#re-defintions for newly loaded data

diff_cnn_ = Conv(reshape(pstar[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
diff_cnn(x) = diff_cnn_(x) .- diff_cnn_.bias
D0 = res3.minimizer[end]

fig = figure(figsize=(4,4))

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 10
rcParams["text.usetex"] = true
rcParams["font.family"] = "serif"
rcParams["font.sans-serif"] = "Helvetica"
rcParams["axes.titlesize"] = 10

subplot(221)
pcolormesh(x,t,ode_data', rasterized=true)
xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
yticks([0, 1, 2, 3, 4, 5])

ax = subplot(222)
cur_pred = predict_rd(pstar)[1]
img = pcolormesh(x,t,cur_pred', rasterized=true)
global img
xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
yticks([0, 1, 2, 3, 4, 5])
cax = fig.add_axes([.48,.62,.02,.29])
colb = fig.colorbar(img, cax=cax)
colb.ax.set_title(L"$\rho$")
clim([0, 1]);
colb.set_ticks([0, 1])

subplot(223)
plot(Flux.data(w1_arr ./ w3_arr) .- 1, label=L"$w_1/w_3 - 1$")
plot(Flux.data(w1_arr .+ w2_arr .+ w3_arr), label=L"$w_1 + w_2 + w_3$")
axhline(0.0, linestyle="--", color="k")
xlabel("Epochs"); title("CNN Weights")
xticks([0, 1500, 3000]); yticks([-0.4, -0.3,-0.2, -0.1, 0.0, 0.1])
legend(loc="lower right", frameon=false, fontsize=6)

subplot(224)
u = collect(0:0.01:1)
plot(u, rx_nn.([[elem] for elem in u]), label="UPDE")[1];
plot(u, reaction.(u), linestyle="--", label="True")
xlabel(L"$\rho$")
title("Reaction Term")
legend(loc="lower center", frameon=false, fontsize=6);
ylim([0, 0.3])

tight_layout(h_pad=1)
gcf()
savefig(@sprintf("%s/fisher_kpp.pdf", save_folder))

#plot loss vs epochs and save
figure(figsize=(6,3))
plot(log.(train_arr), "k.", markersize=1)
xlabel("Epochs"); ylabel("Log(loss)")
tight_layout()
savefig(@sprintf("%s/loss_vs_epoch.pdf", save_folder))
gcf()


text:
"""
    shuffle!(df::DataFrame)

Shuffles a dataframe in place.
"""
function Base.shuffle!(df::DataFrame)
    permutation = shuffle(collect(1:size(df)[1]))
    tdf = copyColumns(df)
    for i in 1:length(permutation)
        df[i, :] = tdf[permutation[i], :]
    end
    return df
end
export shuffle!

"""
    convertNulls!{T}(A::Array{T, 1}, newvalue::T)

Converts all null values to a particular value.
"""
function convertNulls!(A::Vector{T}, newvalue::T) where T<:AbstractFloat
    A[isnan.(A) .| isnull(A)] .= newvalue
    A
end
export convertNulls!

"""
    convertNulls!(df::DataFrame, cols::Vector{Symbol}, newvalue::Any)

Convert all null values in columns of a DataFrame to a particular value.

There is also a method for passing a single column symbol, not as a vector.
"""
function convertNulls!(df::DataFrame, cols::Vector{Symbol}, newvalue::Any)
    for col in cols
        df[col] = convertNulls(df[col], newvalue)
    end
    return
end
convertNulls!(df::DataFrame, col::Symbol, newvalue) = convertNulls!(df, [col], newvalue)
export convertNulls!

"""
    copyColumns(df::DataFrame)

The default copy method for dataframes only copies one level deep, so basically it stores
an array of columns.  If you assign elements of individual (column) arrays then, it can
make changes to references of those arrays that exist elsewhere.

This method instead creates a new dataframe out of copies of the (column) arrays.

This is not named copy due to the fact that there is already an explicit copy(::DataFrame)
implementation in dataframes.

Note that deepcopy is recursive, so this is *NOT* the same thing as deepcopy(df), which
copies literally everything.
"""
function copyColumns(df::DataFrame)
    ndf = DataFrame()
    for col in names(df)
        ndf[col] = copy(df[col])
    end
    return ndf
end
export copyColumns

"""
    applyCatConstraints(dict, df[, kwargs])

Returns a copy of the dataframe `df` with categorical constraints applied.  `dict` should
be a dictionary with keys equal to column names in `df` and values equal to the categorical
values that column is allowed to take on.  For example, to select gauge bosons we can
pass `Dict(:PID=>[i for i in 21:24; -24])`.  Alternatively, the values in the dictionary
can be functions which return boolean values, in which case the returned dataframe will
be the one with column values for which the functions return true.

Note that this requires that the dictionary values are either `Vector` or `Function`
(though one can of course mix the two types).

Alternatively, instead of passing a `Dict` one can pass keywords, for example
`applyCatConstraints(df, PID=[i for i in 21:24; -24])`.
"""
function applyCatConstraints(dict::Dict, df::DataFrame)
    constr = Bool[true for i in 1:size(df)[1]]
    for (col, values) in dict
        constr &= if typeof(values) <: Vector
            convert(BitArray, map(x -> x ∈ values, df[col]))
        elseif typeof(values) <: Function
            convert(BitArray, map(values, df[col]))
        else
            throw(ArgumentError("Constraints must be either vectors or functions."))
        end
    end
    return df[constr, :]
end

function applyCatConstraints(df::DataFrame; kwargs...)
    dct = Dict(kwargs)
    applyCatConstraints(dct, df)
end
export applyCatConstraints

"""
    randomData(dtypes...; nrows=10^4)

Creates a random dataframe with columns of types specified by `dtypes`.  This is useful
for testing various dataframe related functionality.
"""
function randomData(dtypes::DataType...; nrows::Integer=10^4,
                    names::Vector{Symbol}=Symbol[])::DataFrame
    df = DataFrame()
    for (idx, dtype) in enumerate(dtypes)
        col = Symbol(string("col_", idx))
        if dtype <: Real
            df[col] = rand(dtype, nrows)
        elseif dtype <: AbstractString
            df[col] = [randstring(rand(8:16)) for i in 1:nrows]
        elseif dtype <: Dates.TimeType
            df[col] = [dtype(now()) + Dates.Day(i) for i in 1:nrows]
        elseif dtype <: Symbol
            df[col] = [Symbol(randstring(rand(4:12))) for i in 1:nrows]
        end
    end
    if length(names) > 0
        names!(df, names)
    end
    return df
end
export randomData

nans2nulls!(X::Array{Union{T,Null}}) where T = (X[find(isnan(X))] .= null; X)
nans2nulls!(data::DataFrame, col::Symbol) = nans2nulls!(data[col])

"""
    nans2nulls(col)
    nans2nulls(df, colname)

Converts all `NaN`s appearing in the column to `null`.
"""
nans2nulls(col::AbstractVector{T}) where T = nans2nulls!(convert(Array{Union{T,Null}}, col))

nans2nulls(df::DataFrame, col::Symbol) = nans2nulls(df[col])
export nans2nulls

"""
    getCategoryVector(A, vals[, dtype])

Get a vector which is 1 for each `a ∈ A` that satisfies `a ∈ vals`, and 0 otherwise.
Any nulls will be mapped to 0.

Optionally, one can specify the datatype of the output vector.
"""
function getCategoryVector(A::Vector{T}, vals::Vector{T}, ::Type{U}=Int64) where {T,U}
    # this is for efficiency
    valsdict = Dict{T, Void}(v=>nothing for v ∈ vals)
    Vector{U}([a ∈ keys(valsdict) for a ∈ A])
end

function getCategoryVector(A::Vector{Union{T,Null}}, vals::Vector{T}, ::Type{U}=Int64) where {T,U}
    # this is for efficiency
    valsdict = Dict{T, Void}(v=>nothing for v ∈ vals)
    o = map(a -> a ∈ keys(valsdict), A)
    convert(Array{U}, convert(Array, o, 0))
end

function getCategoryVector(A::Vector{T}, val::T, ::Type{U}=Int64) where {T,U}
    getCategoryVector(A, [val], U)
end

function getCategoryVector(A::Vector{Union{T,Null}}, val::T, ::Type{U}=Int64) where {T,U}
    getCategoryVector(A, [val], U)
end

function getCategoryVector(df::AbstractDataFrame, col::Symbol, vals::Vector, ::Type{U}=Int64) where U
    getCategoryVector(df[col], vals, U)
end

function getCategoryVector(df::AbstractDataFrame, col::Symbol, val, ::Type{U}=Int64) where U
    getCategoryVector(df[col], [val], U)
end
export getCategoryVector

"""
    getMatrixDict([T,] df, keycols, datacols)

Gets a dictionary the keys of which are the keys of a groupby of `df` by the columns
`keycols` and the values of which are the matrices produced by taking `sdf[datacols]`
of each `SubDataFrame` `sdf` in the groupby.  Note that the keys are always tuples
even if `keycols` only has one element.

If a type `T` is provided, the output matrices will be of type `Matrix{T}`.
"""
function getMatrixDict(df::DataFrame, keycols::Vector{Symbol}, datacols::Vector{Symbol})
    keycoltypes = getUnwrappedColumnElTypes(df, keycols)
    dict = Dict{Tuple{keycoltypes...},Matrix}()
    for sdf ∈ groupby(df, keycols)
        key = tuple(convert(Array{Any}, sdf[1, keycols])...)
        dict[key] = convert(Array, sdf[datacols])
    end
    dict
end

function getMatrixDict(::Type{T}, gdf::GroupedDataFrame, keycols::Vector{Symbol},
                       datacols::Vector{Symbol}) where  T
    keycoltypes = getUnwrappedColumnElTypes(gdf.parent, keycols)
    dict = Dict{Tuple{keycoltypes...},Matrix{T}}()
    for sdf ∈ gdf
        key = tuple(convert(Array{Any}, sdf[1, keycols])...)
        dict[key] = convert(Array{T}, sdf[datacols])
    end
    dict
end

function getMatrixDict(::Type{T}, gdf::GroupedDataFrame, keycols::Vector{Symbol},
                       Xcols::Vector{Symbol}, ycols::Vector{Symbol}) where T
    keycoltypes = getUnwrappedColumnElTypes(gdf.parent, keycols)
    Xdict = Dict{Tuple{keycoltypes...},Matrix{T}}()
    ydict = Dict{Tuple{keycoltypes...},Matrix{T}}()
    for sdf ∈ gdf
        key = tuple(convert(Array{Any}, sdf[1, keycols])...)
        Xdict[key] = convert(Array{T}, sdf[Xcols])
        ydict[key] = convert(Array{T}, sdf[ycols])
    end
    Xdict, ydict
end

function getMatrixDict(::Type{T}, df::DataFrame, keycols::Vector{Symbol},
                       datacols::Vector{Symbol}) where T
    getMatrixDict(T, groupby(df, keycols), keycols, datacols)
end

# this version is used by grouped dataframe
function getMatrixDict(::Type{T}, df::DataFrame, keycols::Vector{Symbol},
                       Xcols::Vector{Symbol}, ycols::Vector{Symbol}) where T
    getMatrixDict(T, groupby(df, keycols), keycols, Xcols, ycols)
end

export getMatrixDict

#=========================================================================================
Convert an array of Julia objects to a datatable
=========================================================================================#
"""
    struct_to_table(s)
    struct_to_table(v)

Converts a Julia compound type s into a single row of a `DataFrame` with column names
equal to field names and row values equal to field values.

Alternatively, once could pass a vector of compound type instances to create a `DataFrame`
with each row corresponding to an instance of the type.
"""
function struct_to_table(t::T) where T
    fnames = fieldnames(T)
    dat = Vector{Any}(length(fnames))
    for i ∈ 1:length(fnames)
        dat[i] = [getfield(t, fnames[i])]
    end
    DataFrame(dat, fnames)
end

function struct_to_table(v::AbstractVector{T}) where T
    fnames = fieldnames(T)
    dtypes = (fieldtype(T, n) for n ∈ fnames)
    dat = Any[Vector{dtype}(length(v)) for dtype ∈ dtypes]

    for i ∈ 1:length(v), j ∈ 1:length(fnames)
        dat[j][i] = getfield(v[i], fnames[j])
    end

    DataFrame(dat, fnames)
end
export struct_to_table

"""
    Dict(df, keycol, valcol)

Create a dictionary out of the specified columns of a dataframe.
"""
function Dict(df::DataFrame, keycol::Symbol, valcol::Symbol)::Dict
    Dict(df[keycol], df[valcol])
end
export Dict

"""
    datasize(df)

Returns the total size of data stored by a dataframe in megabytes.
"""
datasize(df::DataFrame)::Float64 = sum(sizeof.(df.columns))/10^6
export datasize


text:
using DataFrames

struct DBPerson
    database::SQLite.DB
    pkey::Int
    name::String
    email::String
end

Base.show(io::IO, per::DBPerson) = print(io, per.name)

function NeXLUncertainties.asa(::Type{DataFrame}, people::AbstractArray{DBPerson})
    return DataFrame(
        PKey = [ person.pkey for person in people ],
        Name = [ person.name for person in people ],
        EMail = [ person.email for person in people ]
    )
end

function Base.write(db::SQLite.DB, ::Type{DBPerson}, name::String, email::String)::Int
    (find(db, DBPerson, email) == -1) || error("A person with e-mail address $email already exists.")
    stmt1 = SQLite.Stmt(db, "INSERT INTO PERSON ( NAME, EMAIL ) VALUES ( ?, ? );")
    q = DBInterface.execute(stmt1, ( name, lowercase(email) ))
    return DBInterface.lastrowid(q)
end

function Base.read(db::SQLite.DB, ::Type{DBPerson}, pkey::Int)::DBPerson
    stmt1 = SQLite.Stmt(db, "SELECT PKEY, NAME, EMAIL FROM PERSON WHERE PKEY=?;")
    q = DBInterface.execute(stmt1, ( pkey, ))
    if SQLite.done(q)
        error("No known person with key '$(pkey)'.")
    end
    r = SQLite.Row(q)
    return DBPerson(db, r[:PKEY], r[:NAME], r[:EMAIL])
end

function Base.findall(db::SQLite.DB, ::Type{DBPerson})::Vector{DBPerson}
    stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM PERSON;")
    q = DBInterface.execute(stmt1)
    return [ read(db, DBPerson, r[:PKEY]) for r in q]
end

function find(db::SQLite.DB, ::Type{DBPerson}, email::String)::Int
    stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM PERSON WHERE EMAIL=?;")
    q = DBInterface.execute(stmt1, ( lowercase(email), ))
    return SQLite.done(q) ? -1 : SQLite.Row(q)[:PKEY]
end

function Base.read(db::SQLite.DB, ::Type{DBPerson}, email::String)::DBPerson
    stmt1 = SQLite.Stmt(db, "SELECT PKEY, NAME, EMAIL FROM PERSON WHERE EMAIL=?;")
    q = DBInterface.execute(stmt1, ( email, ))
    if SQLite.done(q)
        error("No known person with e-mail '$(email)'.")
    end
    r = SQLite.Row(q)
    return DBPerson(db, r[:PKEY], r[:NAME], r[:EMAIL])
end


text:
"""
    abstract type AbstractHyperRectangle{N,T} <: AbstractElement{ReferenceHyperCube{N},SVector{N,T}}

Axis-aligned hyperrectangle in `N` dimensions with coordinates of type
`SVector{N,T}`.
"""
abstract type AbstractHyperRectangle{N,T} <: AbstractElement{ReferenceHyperCube{N},SVector{N,T}} end

Base.eltype(::AbstractHyperRectangle{N,T}) where {N,T}     = T
ambient_dimension(::AbstractHyperRectangle{N}) where {N}   = N
geometric_dimension(::AbstractHyperRectangle{N}) where {N} = N
low_corner(r::AbstractHyperRectangle)  = abstractmethod(r)
high_corner(r::AbstractHyperRectangle) = abstractmethod(r)
diameter(r::AbstractHyperRectangle)    = abstractmethod(r)
center(r::AbstractHyperRectangle)      = abstractmethod(r)
radius(r::AbstractHyperRectangle)      = abstractmethod(r)
vertices(r::AbstractHyperRectangle)    = abstractmethod(r)

Base.in(point,h::AbstractHyperRectangle) = all(low_corner(h) .<= point .<= high_corner(h))

function vertices(rec::AbstractHyperRectangle{2})
    lc = low_corner(rec)
    hc = high_corner(rec)
    return SVector(
        SVector(lc[1],lc[2]),
        SVector(hc[1],lc[2]),
        SVector(hc[1],hc[2]),
        SVector(lc[1],hc[2])
    )
end
function vertices(rec::AbstractHyperRectangle{3})
    lc = low_corner(rec)
    hc = high_corner(rec)
    return SVector(
        # lower face
        SVector(lc[1],lc[2],lc[3]),
        SVector(hc[1],lc[2],lc[3]),
        SVector(hc[1],hc[2],lc[3]),
        SVector(lc[1],hc[2],lc[3]),
        # upper face
        SVector(lc[1],lc[2],hc[3]),
        SVector(hc[1],lc[2],hc[3]),
        SVector(hc[1],hc[2],hc[3]),
        SVector(lc[1],hc[2],hc[3])
    )
end

## implement the AbstractElement interface for convenience
function (el::AbstractHyperRectangle)(u)
    @assert u in domain(el)
    lc  = low_corner(el)
    hc  = high_corner(el)
    # map from reference domain to
    v = @. lc + (hc - lc)*u
    return v
end

function jacobian(el::AbstractHyperRectangle,u)
    @assert u in domain(el)
    lc  = low_corner(el)
    hc  = high_corner(el)
    return SDiagonal(hc - lc)
end

"""
    distance(r1::AbstractHyperRectangle,r2::AbstractHyperRectangle)

The (minimal) Euclidean distance between a point `x ∈ r1` and `y ∈ r2`.
"""
function distance(rec1::AbstractHyperRectangle{N},rec2::AbstractHyperRectangle{N}) where {N}
    d2 = 0
    rec1_low_corner  = low_corner(rec1)
    rec1_high_corner = high_corner(rec1)
    rec2_low_corner  = low_corner(rec2)
    rec2_high_corner = high_corner(rec2)
    for i=1:N
        d2 += max(0,rec1_low_corner[i] - rec2_high_corner[i])^2 +
              max(0,rec2_low_corner[i] - rec1_high_corner[i])^2
    end
    return sqrt(d2)
end

"""
    distance(x::SVector,r::HyperRectangle)

The (minimal) Euclidean distance between the point `x` and any point `y ∈ r`.
"""
function distance(pt::SVector{N},rec::AbstractHyperRectangle{N}) where {N}
    d2 = 0
    rec_low_corner  = low_corner(rec)
    rec_high_corner = high_corner(rec)
    for i=1:N
        d2 += max(0,pt[i] - rec_high_corner[i])^2 +
              max(0,rec_low_corner[i] - pt[i])^2
    end
    return sqrt(d2)
end
distance(rec::AbstractHyperRectangle{N},pt::SVector{N}) where {N} = distance(pt,rec)

######
# HyperRectangle
######

"""
    struct HyperRectangle{N,T}

Axis-aligned hyperrectangle in `N` dimensions given by `low_corner::SVector{N,T}` and
`high_corner::SVector{N,T}`.
"""
struct HyperRectangle{N,T} <: AbstractHyperRectangle{N,T}
    low_corner::SVector{N,T}
    high_corner::SVector{N,T}
end
HyperRectangle(l::Tuple,h::Tuple)     = HyperRectangle(SVector(l),SVector(h))
HyperRectangle(l::SVector,h::SVector) = HyperRectangle(promote(l,h)...)
# 1d case
HyperRectangle(a::Number,b::Number) = HyperRectangle(SVector(a),SVector(b))

low_corner(r::HyperRectangle)  = r.low_corner
high_corner(r::HyperRectangle) = r.high_corner
center(r::HyperRectangle)   = (low_corner(r) + high_corner(r)) / 2
diameter(r::HyperRectangle) = norm(high_corner(r) .- low_corner(r),2)
radius(r::HyperRectangle)   = diameter(r) / 2
Base.isapprox(h1::HyperRectangle,h2::HyperRectangle;kwargs...) = isapprox(h1.low_corner,h2.low_corner;kwargs...) && isapprox(h1.high_corner,h2.high_corner;kwargs...)

######
# HyperCube
######

"""
    struct HyperCube{N,T}

Axis-aligned hypercube in `N` dimensions given by `low_corner::SVector{N,T}` and
`side::T`.
"""
struct HyperCube{N,T} <: AbstractHyperRectangle{N,T}
    low_corner::SVector{N,T}
    side::T
end
# 1d case
HyperCube(low_corner::Number,side::Number) = HyperCube(SVector(low_corner),side)

side(r::HyperCube) = r.side
low_corner(r::HyperCube)  = r.low_corner
high_corner(r::HyperCube) = low_corner(r) .+ side(r)
center(r::HyperCube) = low_corner(r) .+ side(r)/2
diameter(r::HyperCube{N}) where N = side(r)*sqrt(N)
radius(r::HyperCube) = diameter(r) / 2
Base.isapprox(h1::HyperCube,h2::HyperCube;kwargs...) = isapprox(h1.low_corner,h2.low_corner;kwargs...) && isapprox(h1.side,h2.side;kwargs...)

######
# Utils
######

function HyperRectangle(els,cube=false)
    isempty(els)  && (error("data cannot be empty") )
    lb  = first(els) |> coords
    ub  = first(els) |> coords
    for el in els
        pt = coords(el)
        lb  = min.(lb,pt)
        ub  = max.(ub,pt)
    end
    if cube # fit a square/cube instead
        w  = maximum(ub-lb)
        xc = (ub + lb) / 2
        lb = xc .- w/2
        ub = xc .+ w/2
        # TODO: return HyperCube instead
    end
    return HyperRectangle(lb,ub)
end

"""
    split(rec::AbstractHyperRectangle,[axis]::Int,[place])

Split a hyperrectangle in two along the `axis` direction at the  position
`place`. Returns a tuple with the two resulting hyperrectangles.

When no `place` is given, defaults to splitting in the middle of the axis.

When no axis and no place is given, defaults to splitting along the largest axis.
"""
function Base.split(rec::AbstractHyperRectangle,axis,place)
    rec_low_corner  = low_corner(rec)
    rec_high_corner = high_corner(rec)
    N            = ambient_dimension(rec)
    high_corner1 = svector(n-> n==axis ? place : rec_high_corner[n], N)
    low_corner2  = svector(n-> n==axis ? place : rec_low_corner[n], N)
    rec1         = HyperRectangle(rec_low_corner, high_corner1)
    rec2         = HyperRectangle(low_corner2,rec_high_corner)
    return (rec1, rec2)
end
function Base.split(rec::AbstractHyperRectangle,axis)
    place = (high_corner(rec)[axis] + low_corner(rec)[axis])/2
    return split(rec,axis,place)
end
function Base.split(rec::AbstractHyperRectangle)
    axis = argmax(high_corner(rec) .- low_corner(rec))
    return split(rec,axis)
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# Hyper parameters
ωs = [1.0, 0.0, 0.0]
ωp = 1e-7
# PO₄ mean and variance of observations fom WOA18
μDIPobs3D, σ²DIPobs3D = WorldOceanAtlasTools.fit_to_grid(grd, "PO₄")
μDIPobs, σ²DIPobs = μDIPobs3D[iwet], σ²DIPobs3D[iwet]
μx = (μDIPobs, missing, missing)
σ²x = (σ²DIPobs, missing, missing)
δconvert(x) = @. exp(x)
cs = (δconvert, identity, identity)
# generate mismatch functions
f, ∇ₓf = f_and_∇ₓf(ωs, μx, σ²x, v, ωp, TestParameters)

#===========================================
Tests
===========================================#

# TODO add tests for other types of mismatch functions!

@testset "Mismatch / cost functions" begin
    nt = length(T_all)
    nx = nt * nb
    np = length(p)
    x = p.xgeo * ones(nx)
    λ = p2λ(p)
    fλ = f(x, λ)
    ∇ₓfλ = ∇ₓf(x, λ)
    @testset "Mismatch function" begin
        @test fλ isa Float64
        @test fλ > 0
        @test fλ ≈ f(x, p)
    end
    @testset "Derivative w.r.t the state, x" begin
        @test size(∇ₓfλ) == (1, nx)
        @test ∇ₓfλ ≈ ∇ₓf(x, p)
    end
end


text:
using SafeTestsets

@safetestset "Runge-Kutta Integrators                                                         " begin include("rk_integrators_tests.jl") end
@safetestset "VPRK Integrators                                                                " begin include("vprk_integrators_tests.jl") end
@safetestset "SPARK Integrators                                                               " begin include("spark_integrators_tests.jl") end
@safetestset "Splitting Integrators                                                           " begin include("splitting_integrators_tests.jl") end
@safetestset "Galerkin Integrators                                                            " begin include("galerkin_integrators_tests.jl") end
@safetestset "Stochastic Integrators                                                          " begin include("stochastic_integrators_tests.jl") end


text:
using CircleFinding
using Test

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

max_stars_repo_name: joeynelson/CircleFinding

max_issues_repo_name: joeynelson/CircleFinding

max_forks_repo_name: joeynelson/CircleFinding


text:
using BenchmarkConfigSweeps

BenchmarkConfigSweeps.run(
    joinpath(@__DIR__, "build"),
    joinpath(@__DIR__, "benchmarks.jl"),
    BenchmarkConfigSweeps.env.(
        "OPENBLAS_NUM_THREADS" .=> 1:max(1, min(16, Sys.CPU_THREADS ÷ 2)),
    ),
)


text:
versioninfo()

import Pkg; Pkg.add("Arpack")
import Pkg; Pkg.add("ArgParse")

using LinearAlgebra
using Arpack
using SparseArrays
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--N"
            help = "set Nsize (should be >=4)"
            arg_type = Int
            default = 8
        "--Sz"
            help = "set Sz"
            arg_type = Int
            default = 0
    end
    return parse_args(s)
end

function snoob(x)
    next = 0
    if(x>0)
        smallest = x & -(x)
        ripple = x + smallest
        ones = xor(x,ripple)
        ones = div((ones >> 2),smallest)
        next = ripple | ones
    end
    return next
end

function count_bit(n)
    count = 0
    while (n>0)
        count += n & 1
        n >>= 1
    end
    return count
end

function init_parameters(N,Sz)
    Nup = div(N,2) + Sz
    Nhilbert = binomial(N,Nup)
    ihfbit = 1 << (div(N,2))
    irght = ihfbit-1
    ilft = xor(((1<<N)-1),irght)
    iup = (1<<(N-Nup))-1
    return Nup, Nhilbert, ihfbit, irght, ilft, iup
end

function make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)
    list_1 = zeros(Int,Nhilbert)
    list_ja = zeros(Int,ihfbit)
    list_jb = zeros(Int,ihfbit)
    ii = iup
    ja = 0
    jb = 0
    ia_old = ii & irght
    ib_old = div((ii & ilft),ihfbit)
    list_1[1] = ii
    list_ja[ia_old+1] = ja
    list_jb[ib_old+1] = jb
    ii = snoob(ii)
    for i in 1:Nhilbert-1
        ia = ii & irght
        ib = div((ii & ilft),ihfbit)
        if ib == ib_old
            ja += 1
        else
            jb += ja+1
            ja = 0
        end
        list_1[i+1] = ii
        list_ja[ia+1] = ja
        list_jb[ib+1] = jb
        ia_old = ia
        ib_old = ib
        ii = snoob(ii)
    end
    return list_1, list_ja, list_jb
end

function get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)
    ia = ii & irght
    ib = div((ii & ilft),ihfbit)
    ja = list_ja[ia+1]
    jb = list_jb[ib+1]
    return ja+jb
end

function make_hamiltonian(J1,D1,N,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    listki = ones(Int,(N+1)*Nhilbert)
    loc = ones(Int,(N+1)*Nhilbert)
    elemnt = zeros((N+1)*Nhilbert)
    for k in 0:N
        for i in 0:Nhilbert-1
            listki[k*Nhilbert+i+1] = i+1
        end
    end
    for k in 0:N-1
        isite1 = k
        isite2 = mod((k+1),N)
        is1 = 1<<isite1
        is2 = 1<<isite2
        is0 = is1 + is2
        wght = -2.0*J1[k+1]
        diag = wght*0.5*D1[k+1]
        for i in 0:Nhilbert-1
            ii = list_1[i+1]
            ibit = ii & is0
            if (ibit==0 || ibit==is0)
                elemnt[N*Nhilbert+i+1] -= diag
                loc[N*Nhilbert+i+1] = i+1
            else
                elemnt[N*Nhilbert+i+1] += diag
                loc[N*Nhilbert+i+1] = i+1
                iexchg = xor(ii,is0)
                newcfg = get_ja_plus_jb(iexchg,irght,ilft,ihfbit,list_ja,list_jb)
                elemnt[k*Nhilbert+i+1] = -wght
                loc[k*Nhilbert+i+1] = newcfg+1
            end
        end
    end
    HamCSC = dropzeros(sparse(listki,loc,elemnt))
    return HamCSC
end

function main()
    parsed_args = parse_commandline()
#    println("Parsed args:")
#    for (arg,val) in parsed_args
#        println("  $arg  =>  $val")
#    end

#    N = 16 # should be N>=4
#    Sz = 0

    N = parsed_args["N"]
    Sz = parsed_args["Sz"]

    Nup, Nhilbert, ihfbit, irght, ilft, iup = init_parameters(N,Sz)

    binirght = string(irght,base=2,pad=N)
    binilft = string(ilft,base=2,pad=N)
    biniup = string(iup,base=2,pad=N)

    println("N=",N)
    println("Sz=",Sz)
    println("Nup=",Nup)
    println("Nhilbert=",Nhilbert)
    println("ihfbit=",ihfbit)
    println("irght,binirght=",irght," ",binirght)
    println("ilft,binilft=",ilft," ",binilft)
    println("iup,biniup=",iup," ",biniup)

    @time list_1, list_ja, list_jb = make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)

    #println("list_1=",list_1)
    #println("list_ja=",list_ja)
    #println("list_jb=",list_jb)

    #println("")
    #println("i ii binii ja+jb")
    #for i in 0:Nhilbert-1
    #    ii = list_1[i+1]
    #    binii = string(ii,base=2,pad=N)
    #    ind = get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)
    #    println(i," ",ii," ",binii," ",ind)
    #end

    J1 = ones(N) # J_{ij}>0: AF
    D1 = ones(N) # D_{ij}>0: AF

    @time Ham = make_hamiltonian(J1,D1,N,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)
    #println(Ham)

    @time ene,vec = eigs(Ham, nev=5, which=:LM)
    #println("# GS energy:",ene[1])
    println("# energy:",ene[1]," ",ene[2]," ",ene[3]," ",ene[4]," ",ene[5])
    #vec_sgn = sign(argmax(vec[:,1]))
    #println("# GS wave function:")
    #for i in 0:Nhilbert-1
    #    ii = list_1[i+1]
    #    binii = string(ii,base=2,pad=N)
    #    println(i," ",vec[i+1,1]*vec_sgn," ",binii)
    #end
end

main()


text:
__precompile__()

module MP3

deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(deps)? include(deps) : error("MP3 is not properly installed. Please run: Pkg.build(\"MP3\")")

# package code goes here
using SampledSignals
using FileIO
using FixedPointNumbers
using SIUnits
using SIUnits.ShortUnits

# methods to override
import SampledSignals: nchannels, nframes, samplerate, unsafe_read!, unsafe_write
import FileIO: load, save

# re-export
export load, save
export Hz, kHz, s

# types used for fixed-point 16-bit and 32-bit encoding
typealias PCM16Sample Fixed{Int16, 15}
typealias PCM32Sample Fixed{Int32, 31}

# shortcut to SIQuantity type of Hz
typealias Hertz SIUnits.SIQuantity{Int64,0,0,-1,0,0,0,0,0,0}

type MP3INFO
    nframes::Int64
    nchannels::Int32
    samplerate::Int32
    datatype::DataType
end

"""create an MP3INFO object from given audio buffer"""
function MP3INFO{T}(buf::SampleBuf{T})
    MP3INFO(size(buf, 1), size(buf, 2), buf.samplerate.val, T)
end

include("readers.jl")
include("writers.jl")

function __init__()
    initialize_readers()
    initialize_writers()

    # MP3 files with ID3v1 (or no tags) and ID3v2 tags have different headers
    magic = (UInt8[0xff, 0xfb], UInt8[0x49, 0x44, 0x33])
    add_format(format"MP3", magic, [".mp3"], [:MP3])
end

end # module


text:
module TestDecisionTree

# using Revise
using Test
import CategoricalArrays
using MLJ

# load code to be tested:
import MLJModels
import DecisionTree # MLJModels.DecisionTree_ now available via lazy-loading
using MLJModels.DecisionTree_

# get some test data:
task = load_iris();
X, y = X_and_y(task) # a table and categorical vector

baretree = DecisionTreeClassifier(target_type=String)

baretree.max_depth = 1
fitresult, cache, report = MLJ.fit(baretree, 1, X, y);
baretree.max_depth = -1 # no max depth
fitresult, cache, report = MLJ.update(baretree, 1, fitresult, cache, X, y);
@test fitresult isa MLJ.fitresult_type(baretree)
# in this case decision tree is a perfect predictor:
yhat = MLJ.predict_mode(baretree, fitresult, X);
@test yhat == y

# but pruning upsets this:
baretree.post_prune = true
baretree.merge_purity_threshold=0.1
fitresult, cache, report = MLJ.update(baretree, 2, fitresult, cache, X, y)
yhat = MLJ.predict_mode(baretree, fitresult, X);
@test yhat != y
yhat = MLJ.predict(baretree, fitresult, X);
cross_entropy(y, yhat)

# check preservation of levels:
yyhat = predict_mode(baretree, fitresult, MLJ.selectrows(X, 1:3))
@test CategoricalArrays.levels(yyhat) == CategoricalArrays.levels(y)

info(baretree)

# testing machine interface:
tree = machine(baretree, X, y)
fit!(tree)
yyhat = predict_mode(tree, MLJ.selectrows(X, 1:3))

end
true


text:
# Remember to set JULIA_NUM_THREADS=8
#using Pkg
#Pkg.add("Distributions")
#Pkg.add("ProgressMeter")
#Pkg.add("HDF5")

#using ProfileView
@everywhere begin
    using Distributions
    using LinearAlgebra
    using Base.Threads
    using ProgressMeter
    using HDF5

    const L = 4
    const N = 16*(L^3)

    const nsamples = 32#56
    const equilibration_mc_steps = 100000#40000

    const N_T = 17
    const N_J = 29
    const J_grid = range(-1, 0.4, length=N_J) # [-0.1]
    const T_grid = 10 .^(range(1, -3, length=N_T))

    normal_d = Normal(0,1)
    uniform_d = Uniform(0,1)

    const a1 = [1,0,0]
    const a2 = [0,1,0]
    const a3 = [0,0,1]

    const coords = 1/8. .* [[1,1,1],[5,5,1],[5,1,5],[1,5,5],
        [1,3,3],[5,3,7],[5,7,3],[1,7,7],
        [3,1,3],[3,5,7],[7,5,3],[7,1,7],
        [3,3,1],[3,7,5],[7,3,5],[7,7,1]]

    function create_random_spin_configuration()
        spins = rand(normal_d, (N,3))
        spins ./ sqrt.(sum(spins.*spins, dims=2))
    end

    coordinates = zeros(Float64, L, L, L, 16, 3)
    for lx in 1:L
        for ly in 1:L
            for lz in 1:L
                for i in 1:16
                    coordinates[lx, ly, lz, i, :] = coords[i] + (lx-1)*a1 + (ly-1)*a2 + (lz-1)*a3
                end
            end
        end
    end

    const x_local = 1/sqrt(6) .* [[-2,1,1], [-2,-1,-1], [2,1,-1], [2,-1,1]]
    const y_local = 1/sqrt(2) .* [[0,-1,1], [0,1,-1], [0,-1,-1], [0,1,1]]
    const z_local = 1/sqrt(3) .* [[1,1,1], [1,-1,-1], [-1,1,-1], [-1,-1,1]]

    function pbc(l)
        if l < 1
            l = L
        end
        if l == (L+1)
            l = 1
        end
        l
    end

    function calculate_indices(indices, inverse_indices)
        k = 1
        for lx in 1:L
            for ly in 1:L
                for lz in 1:L
                    for i in 1:16
                        indices[k,:] = [lx, ly, lz, i]
                        inverse_indices[lx, ly, lz, i] = k
                        k += 1
                    end
                end
            end
        end
    end

    indices = zeros(Int64, N, 4)
    inverse_indices = zeros(Int64, L, L, L, 16)
    calculate_indices(indices, inverse_indices)

    function create_bonds()
        bonds = zeros(Int64, 6, N)
        bond_distance = norm(coords[5]-coords[1])
        for site_i in 1:N
            lx, ly, lz, i = indices[site_i,:]
            coords_i = coordinates[lx, ly, lz, i, :]
            n = 1
            for j in 1:16
                for x in -1:1
                    for y in -1:1
                        for z in -1:1
                            d = norm(coords_i - (coords[j] + a1 * (lx+x-1) + a2 * (ly+y-1) + a3 * (lz+z-1)))
                            if abs(d - bond_distance) < 0.00001
                                bonds[n, site_i] = inverse_indices[pbc(lx+x), pbc(ly+y), pbc(lz+z), j]
                                n = n + 1
                            end
                        end
                    end
                end
            end
        end
        return bonds
    end

    function step!(spins, bonds, J, beta)
        site_i = rand(1:N)
        S1 = spins[site_i,:]
        S2 = rand(normal_d, 3)
        normalize!(S2)

        Sj = sum(spins[bonds[:, site_i],:], dims=1)

        dE = 2.0*J*(
                    (S1[1]-S2[1])*Sj[1] +
                    (S1[2]-S2[2])*Sj[2]
                ) - (S1[3]-S2[3])*Sj[3]

        if (dE <= 0) || (rand(uniform_d) < exp(-dE*beta))
            spins[site_i,:] = S2
        end
        return nothing
    end

    function mcs!(spins, bonds, J, beta)
        for i in 1:N
            step!(spins, bonds, J, beta)
        end
        return nothing
    end

    function local_to_global(spins_local)
        spins_global = zeros(Float64, N, 3)
        for k in 1:N
            lx, ly, lz, atom_i = indices[k,:]
            Sx, Sy, Sz = spins_local[k,:]
            i = (atom_i-1) ÷ 4
            spins_global[k, :] = x_local[i+1,:][1] * Sx + y_local[i+1, :][1] * Sy + z_local[i+1, :][1] * Sz
        end
        return spins_global
    end

    function mc_run(params)
        sample_step, J = params
        spins = create_random_spin_configuration()
        bonds = create_bonds()
        results = zeros(Float64, 2, N_T, N, 3)
        for (T_i, T) in enumerate(T_grid)
            beta = 1/T
            for mc_step in 1:equilibration_mc_steps
                mcs!(spins, bonds, J, beta)
            end
            results[1, T_i, :, :] = spins
            results[2, T_i, :, :] = local_to_global(spins)
        end
        return results
    end
end

function run()
    snapshots_filename = string("snapshots_", L, "_", equilibration_mc_steps, "_", nsamples, ".hdf5")
    h5open(snapshots_filename, "w") do fid
        write(fid, "coordinates", coordinates)
    end
    h5write(snapshots_filename, "J", collect(J_grid))
    h5write(snapshots_filename, "T", collect(T_grid))
    h5write(snapshots_filename, "L", L)

    results_local = zeros(Float64, N_J, N_T, nsamples, N, 3)
    results_global = zeros(Float64, N_J, N_T, nsamples, N, 3)
    #results_global = zeros(Float64, nsamples, N_J, N_T, N, 3)

    params = collect(Iterators.product(1:nsamples, J_grid))
    params_i = collect(Iterators.product(1:nsamples, 1:length(J_grid)))

    results = @showprogress 1 "Computing..." pmap(mc_run, params)
    for nsample in 1:nsamples
        for J_i in 1:length(J_grid)
            results_local[J_i,:,nsample, :, :] = results[nsample,J_i][1, :, :, :]
            results_global[J_i,:,nsample, :, :] = results[nsample,J_i][2, :, :, :]
        end
    end

    h5write(snapshots_filename, "spins_local", results_local)
    h5write(snapshots_filename, "spins_global", results_global)
end

println(N, " spins")
run()
#=
function many_step(spins, bonds, n)
    for i in 1:n
        step!(spins, bonds, -0.4, 1)
    end
    return nothing
end

function test_run()
    spins = create_random_spin_configuration()
    bonds = create_bonds()
    #@profview many_step(spins, bonds, 1)
    @time many_step(spins, bonds, 4000000)
end

test_run()
=#
#c = Condition()
#wait(c)


text:
# Plots power spectrum components of approximate Gaussian PSD:
using PyPlot

#nc = 4
#amp = [2.953823582534889, -0.04961614737314587, -1.6900719772140012, -0.21413545794774125]
#b =[1.589510614489265, 3.4729345651560872, 1.8698492189074047, 6.363451205766626]
#omega =[0.3329119766576696, 4.5999414581901386,1.7115890219170917, 2.205939216106191]
#amp = [2.953823582534889, -0.04961614737314587, -0.21413545794774125, -1.690719772140012]
#b =[1.589510614489265, 3.4729345651560872, 6.363451205766626, 1.8698492189074047]
#omega =[0.3329119766576696, 4.5999414581901386, 2.205939216106191,1.1115890219170917]

nc = 3
amp = [2.670143,-0.458705,-1.211438]
b = [1.500258,4.125558,1.744602]
omega = [0.346731,1.390912,1.798436]

nf = 100
f = 2*pi*logspace(-2,2,nf)
psdtot =  zeros(nf)
for i=1:nc
  psd = amp[i]*b[i]*sqrt(2./pi)*(b[i]^2+omega[i]^2+f.^2)./(b[i]^4+(omega[i]^2-f.^2).^2+2.*b[i]^2.*(omega[i]^2+f[i]^2))
  loglog(f,abs(psd))
  psdtot += psd
  loglog(f,psdtot,linestyle="dashed")
  read(STDIN,Char)
end

#plot(f,psdtot)

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
printstyled("  sync!\n", color=:light_green)

# Method of time calculation in old sync!
function sync_times(t::AbstractArray{Int64, 2}, fs::Float64, t_min::Int64, t_max::Int64)
  t_x = t_expand(t, fs)
  ii = vcat(findall(t_x.<t_min), findall(t_x.>t_max))
  return ii
end
sync_times(t::AbstractArray{Int64, 2}, fs::Float64, t_min::Int64) = get_sync_inds(t, fs, t_min, endtime(t, fs))

# Based on code in sync!
function prune_x!(x::SeisIO.FloatArray, x_del::Array{UnitRange, 1})
  nr = size(x_del, 1)
  for i in nr:-1:1
    if !isempty(x_del[i])
      deleteat!(x, x_del[i])
    end
  end
  return nothing
end

printstyled("    sync_t\n", color=:light_green)

# Cases to consider
printstyled("      nothing removed\n", color=:light_green)
# most basic case: a well-formed time matrix
Δ = 20000
ts = 1583455810004000
nx = 40000
fs = sμ/Δ

t = [1 ts; nx 0]
te = endtime(t, Δ)
t_min = ts
t_max = te
(xi, W) = sync_t(t, Δ, t_min, t_max)

@test xi == [1 nx]
@test W == t_win(t, Δ)

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

printstyled("      empty window\n", color=:light_green)
t = [1 ts-nx*Δ; nx 0]
(xi, W) = sync_t(t, Δ, t_min, t_max)

@test isempty(xi)
@test isempty(W)

printstyled("      simple truncations\n", color=:light_green)
# truncate start
t = [1 ts; nx 0]
nclip = 4
(xi, W) = sync_t(t, Δ, ts+nclip*Δ, te)

@test xi == [nclip+1 nx]

# truncate end
t_min = ts
t_max = te-nclip*Δ
(xi, W) = sync_t(t, Δ, t_min, t_max)

@test xi == [1 nx-nclip]

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# a time matrix with one gap equal to n samples
n = 1
gi = 10
t = [1 ts; gi n*Δ; nx 0]
(xi, W) = sync_t(t, Δ, ts, te)

@test xi == [1 gi-1; gi nx-n]

n = 120
gi = 10
t = [1 ts; gi n*Δ; nx 0]
(xi, W) = sync_t(t, Δ, ts, te)

@test xi == [1 gi-1; gi nx-n]

# two gaps, n and m samples
m = 33
gj = 2000
t = [1 ts; gi n*Δ; gj m*Δ; nx 0]
(xi, W) = sync_t(t, Δ, ts, te)

@test xi == [1 gi-1; gi gj-1; gj nx-m-n]

t_min = ts
t_max = te
i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# negative gap
t = [1 ts; gi -n*Δ; nx 0]
(xi, W) = sync_t(t, Δ, ts, te)

@test xi == [1 gi-1; n+1 nx]

# here, the number of samples before ts in window 2 is n-gi, or 110; the first
# n-gi samples of this window are samples gi to n-gi+gi = n; so the window
# starts at n+1, or 121

t_min = ts
t_max = te
i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# truncate other (windows not in chronological order)
n = 120
gi = 10
nx = 100

# start of middle window
W = ts .+ Δ.*([   1       gi
               -n+1        0
               gi+1       nx ] .- 1)
t_min = ts - div(n,2)*Δ
t_max = last(W)
t = w_time(W, Δ)
xi0 = x_inds(t)
(xi, W) = sync_t(t, Δ, t_min, t_max)

@test xi-xi0 == [0 0; div(n,2) 0; 0 0]

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# end of middle window
W = ts .+ Δ.*([   1       gi
               nx+1       nx+n
               gi+1       nx ] .- 1)
t_min = ts
t_max = ts + (nx - 1 + div(n,2))*Δ
t = w_time(W, Δ)
xi0 = x_inds(t)
(xi, W) = sync_t(t, Δ, t_min, t_max)

@test xi-xi0 == [0 0; 0 -div(n,2); 0 0]

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

printstyled("      segment removed\n", color=:light_green)
# first window gets emptied
gi = 10
nx = 100
gl = 6
W = ts .+ Δ.*[  -gl       -1
                  1     gi-1
              gi+gl       nx]
t_min = ts
t_max = last(W)
t = w_time(W, Δ)
xi0 = x_inds(t)
xi, W1 = sync_t(t, Δ, ts, last(W))

@test size(xi, 1) == size(W1, 1) == 2
@test sum(diff(xi, dims=2).+1) == nx-gl

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# second window gets emptied (no gap, windows 1 and 2 are out of order)
W = ts .+ Δ.*[    0     gi-1
                -gl       -1
                 gi  nx-gl-1]
t = w_time(W, Δ)
xi = x_inds(t)
xi1, W1 = sync_t(t, Δ, ts, last(W))

@test size(xi1, 1) == size(W1, 1) == 2
@test xi1 == xi[[1,3],:]

# last window gets emptied
W = ts .+ Δ.*[    0       gi
               gi+2       nx
               nx+2    nx+gi]
t_min = ts
t_max = W[2,2]
t = w_time(W, Δ)
xi = x_inds(t)
xi1, W1 = sync_t(t, Δ, t_min, t_max)

@test size(xi1, 1) == size(W1, 1) == 2
@test xi1 == xi[[1,2],:]

printstyled("      multiple segs removed\n", color=:light_green)
# first + last emptied
gi = 10
W = ts .+ Δ.*[    0     gi-1
               gi+2       nx
               nx+2    nx+gi]
t_min = ts + gi*Δ
t_max = ts + nx*Δ
t = w_time(W, Δ)
xi0 = x_inds(t)
xi, W1 = sync_t(t, Δ, t_min, t_max)

@test size(xi, 1) == size(W1, 1) == 1
@test xi == xi0[[2],:]

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# first + last emptied, middle truncated by os
os = 10
t_max = ts + (nx-os)*Δ
t = w_time(W, Δ)
xi0 = x_inds(t)
xi, W1 = sync_t(t, Δ, t_min, t_max)

@test size(xi, 1) == size(W1, 1) == 1
@test xi0[[2],:] .- xi == [0 os]

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# first emptied, last truncated, middle is one point
gi = 666
nx = gi + 100
W = ts .+ Δ.*[    0     gi-1
               gi+2       nx
               nx+2    nx+gi]
t_min = ts + nx*Δ
t_max = ts + (nx+div(gi,2))*Δ
t = w_time(W, Δ)
xi0 = x_inds(t)
xi, W1 = sync_t(t, Δ, t_min, t_max)

@test size(xi, 1) == size(W1, 1) == 2
@test W1[1,1] == W1[1,2]
@test div(W1[2,2] - W1[2,1], Δ) + 2 == div(gi, 2)

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# second and last emptied
gi = 10
W = ts .+ Δ.*[    0     gi-1
               gi+2       nx
               nx+2    nx+gi]
t_min = ts + (gi+1)*Δ
t_max = ts + nx*Δ
t = w_time(W, Δ)
xi = x_inds(t)
xi1, W1 = sync_t(t, Δ, t_min, t_max)

@test size(xi1, 1) == size(W1, 1) == 1
@test xi1 == xi[[2],:]

printstyled("      complicated cases\n", color=:light_green)
# complex case with no completely emptied windows; formed as W and not :t
nx = 40000
W = ts .+ Δ.*[1     10
              12    12
              14    14
              23    40
              51    100
              -1200 0
              101   9500
              12000 41318]

t_min = ts - 300Δ
t_max = last(W) - 2000Δ
t = w_time(W, Δ)
xi0 = x_inds(t)
xi, W = sync_t(t, Δ, t_min, t_max)

j = sortperm(W[:,1])
W1 = W[j,:]
t1 = w_time(W1, Δ)
@test t1[1,2] == t_min
@test endtime(t1, Δ) == t_max

# check our work; sorted matrix should be from t_min to t_max
t2 = w_time(W, Δ)
@test diff(x_inds(t2), dims=2) == diff(xi, dims=2)

i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi, t[end, 1])
test_del_ranges(x_del, i2)

# last window emptied; previous window truncated; first window truncated
nx = 1000
gi = 17
os = 50
W = ts .+ Δ.*[ -1*os      gi
               gi+2    nx+gi
               nx+gi+2   2nx ]
t_min = ts
t_max = W[2,2] - Δ*gi
t = w_time(W, Δ)
xi = x_inds(t)
xi1, W1 = sync_t(t, Δ, t_min, t_max)
t1 = w_time(W1, Δ)

@test size(xi1, 1) == size(W1, 1) == 2
@test xi1[2,2] == nx+os
@test W[[1,2],:] != W1
@test xi1[[1,2],:] != xi
@test t1[1,2] == t_min
@test endtime(t1, Δ) == t_max

nx_true = t[end,1]
i2 = sync_times(t, fs, t_min, t_max)
x_del = get_del_ranges(xi1, nx_true)
test_del_ranges(x_del, i2)

# Final step, putting it all together
x = randn(Float32, nx_true)
x1 = deepcopy(x)
x2 = deepcopy(x)
prune_x!(x1, x_del)
deleteat!(x2, i2)
@test x1 == x2

fs = 100.0
Δ = round(Int64, sμ/fs)
t0 = 1582934400000000
ns = 10
Δ3 = 1.0

# Test three seismic and three irregularly sampled channels
S = randSeisData(3, s=1.0)

# seismic channels will be ns seconds of data sampled at 100 Hz with preset gaps
nx = round(Int64, ns * fs)
for i = 1:3
  S.x[i] = rand(nx)       # Set S.x to be 1000 samples each
  S.fs[i] = fs            # Set S.fs to be uniformly 100 Hz
end
S.t[1] = [1 t0; nx 0]
S.t[2] = [1 t0-ns*1000000; nx 0]
S.t[3] = [1 t0; 2 round(Int64, Δ3*sμ); nx 0]

# irregular channels will be 100 points each, randomly sampled between the start and the end
append!(S, randSeisData(3, c=1.0))
n_irr = div(nx, 10)
ni = div(n_irr, ns)
for i = 4:6
  S.x[i] = rand(n_irr)
  S.fs[i] = 0.0
  S.t[i] = zeros(n_irr, 2)
  tv = zeros(Int64, n_irr)

  #= the definition of rmax here prevents rng giving us a value above tmax
    in test (6), which breaks the test
  =#
  for j = 1:ns
    si = 1 + (j-1)*ni
    ei = j*ni
    rmin = t0 + (j-1)*1000000
    rmax = t0 - 1 + (j*1000000) - (j == ns ? 2Δ : 0)
    r = UnitRange{Int64}(rmin:rmax)
    tv[si:ei] .= rand(r, ni)
  end
  sort!(tv)
  S.t[i][:,1] .= 1:n_irr
  S.t[i][:,2] .= tv
end
sx = Lx(S)

# Set derived variables
ts, te = get_edge_times(S)
ts₀ = minimum(ts)
ts₁ = maximum(ts)

# TEST 1 =====================================================================
# Sync to last start; don't sync ends
printstyled(stdout,"    sync to last start; don't sync ends\n", color=:light_green)
T = sync(S)
basic_checks(T)
ts_new, te_new = get_edge_times(T)

# Expectation: T[2] is deleted
@test T.n == 5
@test S.id[1] == T.id[1]
@test findid(S.id[2], T.id) == 0
for i = 3:6
  @test(S.id[i] == T.id[i-1])
end

# Expectation: the latest time in S.t is now the earliest start time
@test minimum(ts_new .≥ ts₁) == true
@test maximum(ts_new .== ts₁) == true

# TEST 2 =====================================================================
# Sync to last start; don't sync ends

# change trace 2 to begin only 1s earlier; resync
S.t[2] = [1 t0-1000000; nx 0]
ts, te = get_edge_times(S)
ts₀ = minimum(ts)
ts₁ = maximum(ts)

T = sync(S)
ts_new, te_new = get_edge_times(T)

# Expectation: T.x[2] is at least 100 pts shorter and exactly 100 shorter than T.x[1]
@test T.n == 6
tx = Lx(T)
@test tx[2] ≤ nx-100
@test tx[2] == tx[1] - 100
@test tx[2] ≤ tx[3]
@test minimum(ts_new .≥ ts₁) == true
@test maximum(ts_new .== ts₁) == true
basic_checks(T)

# TEST 3 =====================================================================
# Sync to last start; sync to first end
printstyled(stdout,"    sync to last start; sync to first end\n", color=:light_green)
te₀ = minimum(te)
te₁ = maximum(te)

T = sync(S, t="first")
basic_checks(T)
ts_new, te_new = get_edge_times(T)

# Expectation: T.x[1] and T.x[2] are same length, T.x[3] is shorter
tx = Lx(T)
@test tx[1] == tx[2]
@test tx[2] > tx[3]
@test minimum(te_new .≤ te₀) == true
@test maximum(te_new .== te₀) == true
@test minimum(sx.-tx.≥0) == true                # vx is always longer

# TEST 4 =====================================================================
# Sync to first start; sync to first end
printstyled(stdout,"    sync to first start; sync to last end\n", color=:light_green)
U = sync(S, s="first", t="first")
basic_checks(U)
ts_new, te_new = get_edge_times(U)

# Expectation:
ux = Lx(U)
@test minimum(ux.-tx.≥0) == true                # vx is always longer
@test tx[1]+100 ≤ ux[1]                        # U.x[1] gains >100 samples
@test minimum(ts_new .≥ ts₀) == true           # Should still be true
@test 2 ∈ findall((ts_new .== ts₀).== true)   # Defined earliest start time
@test ux[2]==ux[1]                            # Due to how we set the gap
@test ux[1]-ux[3] == 100
@test minimum(te_new .≤ te₀) == true
@test maximum(te_new .== te₀) == true
for i = 4:6
  @test ux[i] ≥ tx[i]
end

# TEST 5 =====================================================================
# Sync to first start; sync to last end
V = sync(S, s="first", t="last")
basic_checks(V)
ts_new, te_new = get_edge_times(V)
vx = Lx(V)

# Expectation:
@test minimum(vx.-sx.≥0) == true                # vx is always longer
t1 = t_expand(V.t[1], V.fs[1])
t2 = t_expand(S.t[1], S.fs[1])
j = setdiff(t1, t2)
@test vx[1]-sx[1] == length(j)
t2 = t_expand(U.t[1], U.fs[1])
j = setdiff(t1, t2)
@test vx[1]-ux[1] == length(j)
@test minimum(ts_new .≥ ts₀) == true           # Should still be true
@test 2 ∈ findall((ts_new .== ts₀).== true)   # Defined earliest start time
@test vx[2] == vx[1]                           # Due to how we set the gap
t2 = t_expand(V.t[2], V.fs[2])
t3 = t_expand(V.t[3], V.fs[3])
j = setdiff(t2, t3)
@test vx[2] - vx[3] == length(j)
@test minimum(te_new .≤ te₁) == true
@test maximum(te_new .== te₁) == true
for i = 4:6
  @test vx[i] ≥ tx[i]
  @test vx[i] == sx[i]
end

# TEST 6 =====================================================================
# Sync to s = DateTime (first trace start), t = "none"
# trace 3 should be 100 samples shorter
# so should trace 2
printstyled(stdout,"    sync to DateTime; don't sync ends\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + Δ*(nx-1)
ds₆ = u2d(ts₆*μs)
de₆ = u2d(te₆*μs)
W = sync(S, s=ds₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                  # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

# Repeat with an end time
# te₆ = S.t[1][1,2] + round(Int64, sμ*(nx-1)/fs)
# de₆ = u2d(te₆*μs)
W = sync(S, s=ds₆, t=de₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[3]-wx[3] == 100                  # Trace 3 is 100 samples shorter
for i in [1,2,4,5,6]
  @test sx[i] == wx[i]    # No change in other trace lengths; 2 gets padded
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

# TEST 7 =====================================================================
# Sync to DateTime 1s before first trace, 0.01s after; is it 101 pts longer?

# Repeat with an increased window range
ds₆ -= Second(1)
de₆ += Millisecond(10)
X = sync(S, s=ds₆, t=de₆)
basic_checks(X)
ts_new, te_new = get_edge_times(X)
xx = Lx(X)

# Expectations:
@test xx[1]-sx[1] == 101    # Trace 1 is 101 samples longer
@test xx[2]-sx[2] == 101    # Trace 2 is also 101 samples longer
@test xx[3]-sx[3] == 1      # Trace 3 loses 100 samples at the end, but gains 100 at end and 1 at start ... net gain +1
for i = 4:6
  @test xx[i] == sx[i]
end

# Sync to DateTime 2s after first trace, t="first"
ds₆ += Second(3)
X = sync(S, s=ds₆, t=de₆)
basic_checks(X)
ts_new, te_new = get_edge_times(X)
xx = Lx(X)

# Expectations:
@test xx[1] == 801          # Trace 1 should be 801 samples
@test xx[1] == xx[2]        # Should be the same due to padding
# In fact trace 2 should have 101 points appended
found = findall([occursin("appended 101", i) for i in X.notes[2]])
@test isempty(found) == false

# TEST 8 =====================================================================
# A few simple combinations; do these work?
printstyled(stdout,"    sync start to DateTime; sync first end\n", color=:light_green)
X = sync(S, s=ds₆, t="first"); basic_checks(X)

printstyled(stdout,"    sync start to DateTime; sync to last end\n", color=:light_green)
X = sync(S, s=ds₆, t="last"); basic_checks(X)

printstyled(stdout,"    sync start to string time; sync to last end\n", color=:light_green)
ss = string(ds₆)
Y = sync(S, s=ds₆, t="last", v=3); basic_checks(Y)

# Expect: X != Y due to notes, but all other fields equal
for f in datafields
  if f != :notes
    @test isequal(getfield(X,f), getfield(Y,f))
  end
end

# TEST 9 =====================================================================
# Do we actually prune campaign data when all times are out of range?
printstyled(stdout,"    prune all irregular data when all times are out of range\n", color=:light_green)

ss = string(ds₆)
Z = deepcopy(S)
t1 = deepcopy(Z.t[5])
t1 = hcat(t1[:,1:1], vcat(0, diff(t1[:,2:2], dims=1)))
Z.t[5] = deepcopy(t1)

redirect_stdout(out) do
  sync!(Z, v=3); basic_checks(Z)
end

# Expect: Z[5] is gone
for i in [1,2,3,4,6]
  @test any(Z.id.==S.id[i])
end
@test (any(Z.id.==S.id[5]) == false)

# ===========================================================================
# method extenson to SeisChannel
printstyled(stdout,"    SeisChannel method extension\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + Δ*(nx-1)
ds₆ = u2d(ts₆*μs)
de₆ =  u2d(te₆*μs)
C = deepcopy(S[1])
sync!(C, s=ds₆)
W = SeisData(C)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Repeat with an end time
C = deepcopy(S[1])
W = SeisData(sync(C, s=ds₆, t=de₆))
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# ===========================================================================
# method extenson to SeisEvent
printstyled(stdout,"    SeisEvent method extension\n", color=:light_green)
ts₆ = S.t[1][1,2]
te₆ = S.t[1][1,2] + Δ*(nx-1)
ds₆ = u2d(ts₆*μs)
de₆ =  u2d(te₆*μs)
Ev = SeisEvent(hdr = randSeisHdr(), data = deepcopy(S))
sync!(Ev.data, s=ds₆)
W = Ev.data
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                      # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time

Ev = SeisEvent(hdr = randSeisHdr(), data = deepcopy(S))
W = sync(Ev.data, s=ds₆)
basic_checks(W)
ts_new, te_new = get_edge_times(W)
wx = Lx(W)

# Expectations:
@test sx[2]-wx[2] == 100                      # Trace 2 is 100 samples shorter
for i in [1,3,4,5,6]
  @test sx[i] == wx[i]                        # No change in other trace lengths
end
@test minimum(ts_new .≥ ts₆) == true           # We start at ts₆, not before
@test findfirst(ts_new .== ts₆).== 1          # Defined start time


text:
### Random Variable subtypes

## Abstract parameter

abstract Parameter{S<:ValueSupport, F<:VariateForm} <: Variable{Random}

typealias DiscreteParameter{F<:VariateForm} Parameter{Discrete, F}
typealias ContinuousParameter{F<:VariateForm} Parameter{Continuous, F}

typealias UnivariateParameter{S<:ValueSupport} Parameter{S, Univariate}
typealias MultivariateParameter{S<:ValueSupport} Parameter{S, Multivariate}
typealias MatrixvariateParameter{S<:ValueSupport} Parameter{S, Matrixvariate}

typealias ParameterVector{P<:Parameter} Vector{P}

### Code generation of parameter fields

function codegen_setfield(parameter::Parameter, field::Symbol, f::Function)
  @gensym setfield
  quote
    function $setfield(_state::$(default_state_type(parameter)), _states::VariableStateVector)
      setfield!($parameter, $(QuoteNode(field)), $f(_state, _states))
    end
  end
end

function codegen_closure(parameter::Parameter, f::Function)
  @gensym closure
  quote
    function $closure(_state::$(default_state_type(parameter)))
      $(f)(_state, $(parameter).states)
    end
  end
end

function codegen_target_closure_via_distribution(parameter::Parameter, distribution::Symbol, f::Function, field::Symbol)
  @gensym target_closure_via_distribution
  quote
    function $target_closure_via_distribution(_state::$(default_state_type(parameter)))
      setfield!(_state, $(QuoteNode(field)), $(f)(getfield($parameter, $(QuoteNode(distribution))), _state.value))
    end
  end
end

function codegen_sumtarget_closure(
  parameter::Parameter, plfield::Symbol, ppfield::Symbol, stfield::Symbol, slfield::Symbol, spfield::Symbol
)
  body = []

  push!(body, :(getfield($parameter, $(QuoteNode(plfield)))(_state)))
  push!(body, :(getfield($parameter, $(QuoteNode(ppfield)))(_state)))
  push!(body, :(setfield!(
    _state,
    $(QuoteNode(stfield)),
    getfield(_state, $(QuoteNode(slfield)))+getfield(_state, $(QuoteNode(spfield)))))
  )

  @gensym sumtarget_closure

  quote
    function $sumtarget_closure(_state::$(default_state_type(parameter)))
      $(body...)
    end
  end
end

function codegen_uptotarget_closures(parameter::Parameter, fields::Vector{Symbol})
  body = []
  local f::Symbol

  for i in 1:length(fields)
    f = fields[i]
    push!(body, :(getfield($parameter, $(QuoteNode(f)))(_state)))
  end

  @gensym uptotarget_closures

  quote
    function $uptotarget_closures(_state::$(default_state_type(parameter)))
      $(body...)
    end
  end
end

setpdf!{S<:ValueSupport, F<:VariateForm}(parameter::Parameter{S, F}, state::ParameterState{S, F}) =
  parameter.setpdf(state, parameter.states)

setprior!{S<:ValueSupport, F<:VariateForm}(parameter::Parameter{S, F}, state::ParameterState{S, F}) =
  parameter.setprior(state, parameter.states)

value_support{S<:ValueSupport, F<:VariateForm}(::Type{Parameter{S, F}}) = S
variate_form{S<:ValueSupport, F<:VariateForm}(::Type{Parameter{S, F}}) = F

function check_support{S<:ValueSupport, F<:VariateForm}(parameter::Parameter{S, F}, state::ParameterState{S, F})
  if value_support(parameter) != value_support(state)
    warn("Value support of parameter ($(value_support(parameter))) and of ($(value_support(state))) not in agreement")
  end

  if variate_form(parameter) != variate_form(state)
    error("Variate form of parameter ($(variate_form(parameter))) and of ($(variate_form(state))) not in agreement")
  end
end

dotshape(variable::Parameter) = "circle"