text:
using Parquet
using Test
using Dates

function test_load(file::String)
    p = ParFile(file)
    println("loaded $file")
    @test isa(p.meta, Parquet.FileMetaData)

    rg = rowgroups(p)
    @test length(rg) > 0
    println("\tfound $(length(rg)) row groups")

    cc = columns(p, 1)
    println("\tfound $(length(cc)) column chunks in row group 1")

    cnames = colnames(rg[1])
    @test length(cnames) == length(cc)
    println("\tcolumns: $cnames")

    pg = pages(p, 1, 1)
    println("\tfound $(length(pg)) pages in column chunk 1 of row group 1")

    bpg = bytes(pg[1])
    println("\tread $(length(bpg)) bytes from page 1, column chunk 1 of row group 1")

    cc1 = columns(p, rg[1], cnames[1:2:3])
    @test length(cc) > length(cc1)
    @test length(rowgroups(p, [cnames[1], ["xyz"]])) == 0

    iob = IOBuffer()
    show(iob, p)
    sb = take!(iob)
    @test !isempty(sb)
    println("\t" * String(sb))

    println("\tsuccess")
end

function test_load_all_pages()
    testfolder = joinpath(@__DIR__, "parquet-compatibility")
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            testfile = joinpath(testfolder, "parquet-testdata", "impala", "1.1.1-$encformat", "$fname.impala.parquet")
            test_load(testfile)
        end
    end

    testfolder = joinpath(@__DIR__, "julia-parquet-compatibility")
    for encformat in ("ZSTD", "SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            testfile = joinpath(testfolder, "Parquet_Files", "$(encformat)_pandas_pyarrow_$(fname).parquet")
            test_load(testfile)
        end
    end
end

function test_load_boolean_and_ts()
    println("testing booleans and timestamps...")
    p = ParFile(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"))

    rg = rowgroups(p)
    @test length(rg) == 1
    cc = columns(p, 1)
    @test length(cc) == 11
    cnames = colnames(rg[1])
    @test length(cnames) == length(cc)
    @test cnames[2] == ["bool_col"]
    pg = pages(p, 1, 1)
    @test length(pg) == 2

    rc = RecordCursor(p; rows=1:2, colnames=colnames(p))
    @test length(rc) == 2
    @test eltype(rc) == NamedTuple{(:id, :bool_col, :tinyint_col, :smallint_col, :int_col, :bigint_col, :float_col, :double_col, :date_string_col, :string_col, :timestamp_col),Tuple{Union{Missing, Int32},Union{Missing, Bool},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int64},Union{Missing, Float32},Union{Missing, Float64},Union{Missing, Array{UInt8,1}},Union{Missing, Array{UInt8,1}},Union{Missing, Int128}}}

    values = Any[]
    for rec in rc
       push!(values, rec)
    end

    @test [v.bool_col for v in values] == [true,false]
    @test [logical_timestamp(v.timestamp_col) for v in values] == [DateTime("2009-04-01T12:00:00"), DateTime("2009-04-01T12:01:00")]
    @test [logical_timestamp(v.timestamp_col; offset=Dates.Second(30)) for v in values] == [DateTime("2009-04-01T12:00:30"), DateTime("2009-04-01T12:01:30")]
    @test [logical_string(v.date_string_col) for v in values] == ["04/01/09", "04/01/09"]
    #dlm,headers=readdlm("booltest/alltypes.csv", ','; header=true)
    #@test [v.bool_col for v in values] == dlm[:,2]  # skipping for now as this needs additional dependency on DelimitedFiles
end

function test_logical_type_mapping()
    println("testing default logical type mapping...")

    p = ParFile(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"); map_logical_types=true)
    rc = RecordCursor(p)
    @test eltype(rc) == NamedTuple{(:id, :bool_col, :tinyint_col, :smallint_col, :int_col, :bigint_col, :float_col, :double_col, :date_string_col, :string_col, :timestamp_col),Tuple{Union{Missing, Int32},Union{Missing, Bool},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int64},Union{Missing, Float32},Union{Missing, Float64},Union{Missing, String},Union{Missing, String},Union{Missing, DateTime}}}

    values = collect(rc)
    @test [v.bool_col for v in values] == [true,false]
    @test [v.timestamp_col for v in values] == [DateTime("2009-04-01T12:00:00"), DateTime("2009-04-01T12:01:00")]
    @test [v.date_string_col for v in values] == ["04/01/09", "04/01/09"]

    p = ParFile(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"); map_logical_types=Dict(["timestamp_col"]=>(DateTime,v->logical_timestamp(v; offset=Dates.Second(30)))))
    rc = RecordCursor(p)
    @test eltype(rc) == NamedTuple{(:id, :bool_col, :tinyint_col, :smallint_col, :int_col, :bigint_col, :float_col, :double_col, :date_string_col, :string_col, :timestamp_col),Tuple{Union{Missing, Int32},Union{Missing, Bool},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int64},Union{Missing, Float32},Union{Missing, Float64},Union{Missing, Array{UInt8,1}},Union{Missing, Array{UInt8,1}},Union{Missing, DateTime}}}

    values = collect(rc)
    @test [v.timestamp_col for v in values] == [DateTime("2009-04-01T12:00:30"), DateTime("2009-04-01T12:01:30")]

    cc = BatchedColumnsCursor(p)
    values, _state = iterate(cc)
    @test values.timestamp_col == [DateTime("2009-04-01T12:00:30"), DateTime("2009-04-01T12:01:30")]
end

function test_load_nested()
    println("testing nested columns...")
    p = ParFile(joinpath(@__DIR__, "nested", "nested1.parquet"))

    @test nrows(p) == 100
    @test ncols(p) == 5

    rc = RecordCursor(p)
    @test length(rc) == 100
    @test eltype(rc) == NamedTuple{(:_adobe_corpnew,),Tuple{NamedTuple{(:id, :vocab, :frequency, :max_len, :reduced_max_len),Tuple{Union{Missing, Int32},Union{Missing, Array{UInt8,1}},Union{Missing, Int32},Union{Missing, Float64},Union{Missing, Int32}}}}}

    values = Any[]
    for rec in rc
        push!(values, rec)
    end

    v = values[1]._adobe_corpnew
    @test v.frequency == 3
    @test v.id == 1375
    @test v.max_len == 64192.0
    @test v.reduced_max_len == 64
    @test String(copy(v.vocab)) == "10385911_a"

    v = values[100]._adobe_corpnew
    @test v.frequency == 61322
    @test v.id == 724
    @test v.max_len == 64192.0
    @test v.reduced_max_len == 64
    @test String(copy(v.vocab)) == "12400277_a"

    p = ParFile(joinpath(@__DIR__, "nested", "nested.parq"))

    @test nrows(p) == 10
    @test ncols(p) == 1

    rc = RecordCursor(p)
    @test length(rc) == 10
    @test eltype(rc) == NamedTuple{(:nest,),Tuple{Union{Missing, NamedTuple{(:thing,),Tuple{Union{Missing, NamedTuple{(:list,),Tuple{Array{NamedTuple{(:element,),Tuple{Union{Missing, Array{UInt8,1}}}},1}}}}}}}}}

    values = collect(rc)
    v = first(values)
    @test length(v.nest.thing.list) == 1
    @test v.nest.thing.list[1].element == UInt8[0x68,0x69]
    v = last(values)
    @test length(v.nest.thing.list) == 1
    @test v.nest.thing.list[1].element == UInt8[0x77,0x6f,0x72,0x6c,0x64]
end

test_load_all_pages()
test_load_boolean_and_ts()
test_logical_type_mapping()
test_load_nested()


text:
# # Bayesian Neural Network

# We borrow this tutorial from the [official Turing Docs](https://turing.ml/dev/tutorials/03-bayesian-neural-network/). We will show how the explicit parameterization of Lux enables first-class composability with packages which expect flattened out parameter vectors.

# We will use [Turing.jl](https://turing.ml) with [Lux.jl](https://lux.csail.mit.edu/stable) to implement implementing a classification algorithm. Lets start by importing the relevant libraries

## Import libraries.
using Lux
using Pkg #hide
Pkg.activate(joinpath(dirname(pathof(Lux)), "..", "examples")) #hide
using Turing, Plots, Random, ReverseDiff, NNlib, Functors

## Hide sampling progress.
Turing.setprogress!(false);

## Use reverse_diff due to the number of parameters in neural networks.
Turing.setadbackend(:reversediff)

# ## Generating data

# Our goal here is to use a Bayesian neural network to classify points in an artificial dataset. The code below generates data points arranged in a box-like pattern and displays a graph of the dataset we'll be working with.

## Number of points to generate.
N = 80
M = round(Int, N / 4)
rng = Random.default_rng()
Random.seed!(rng, 1234)

## Generate artificial data.
x1s = rand(rng, Float32, M) * 4.5f0;
x2s = rand(rng, Float32, M) * 4.5f0;
xt1s = Array([[x1s[i] + 0.5f0; x2s[i] + 0.5f0] for i in 1:M])
x1s = rand(rng, Float32, M) * 4.5f0;
x2s = rand(rng, Float32, M) * 4.5f0;
append!(xt1s, Array([[x1s[i] - 5.0f0; x2s[i] - 5.0f0] for i in 1:M]))

x1s = rand(rng, Float32, M) * 4.5f0;
x2s = rand(rng, Float32, M) * 4.5f0;
xt0s = Array([[x1s[i] + 0.5f0; x2s[i] - 5.0f0] for i in 1:M])
x1s = rand(rng, Float32, M) * 4.5f0;
x2s = rand(rng, Float32, M) * 4.5f0;
append!(xt0s, Array([[x1s[i] - 5.0f0; x2s[i] + 0.5f0] for i in 1:M]))

## Store all the data for later.
xs = [xt1s; xt0s]
ts = [ones(2 * M); zeros(2 * M)]

## Plot data points.
function plot_data()
    x1 = first.(xt1s)
    y1 = last.(xt1s)
    x2 = first.(xt0s)
    y2 = last.(xt0s)

    plt = Plots.scatter(x1, y1; color="red", clim=(0, 1))
    Plots.scatter!(plt, x2, y2; color="blue", clim=(0, 1))

    return plt
end

plot_data()

# ## Building the Neural Network

# The next step is to define a feedforward neural network where we express our parameters as distributions, and not single points as with traditional neural networks. For this we will use `Dense` to define liner layers and compose them via `Chain`, both are neural network primitives from `Lux`. The network `nn` we will creat will have two hidden layers with `tanh` activations and one output layer with `sigmoid` activation, as shown below.

# The `nn` is an instance that acts as a function and can take data, parameters and current state as inputs and output predictions. We will define distributions on the neural network parameters.

## Construct a neural network using Lux
nn = Chain(Dense(2, 3, tanh), Dense(3, 2, tanh), Dense(2, 1, sigmoid))

## Initialize the model weights and state
ps, st = Lux.setup(rng, nn)

Lux.parameterlength(nn) # number of paraemters in NN

# The probabilistic model specification below creates a parameters variable, which has IID normal variables. The parameters represents all parameters of our neural net (weights and biases).

## Create a regularization term and a Gaussian prior variance term.
alpha = 0.09
sig = sqrt(1.0 / alpha)

# Construct named tuple from a sampled parameter vector. We could also use ComponentArrays here and simply broadcast to avoid doing this. But let's do it this way to avoid dependencies.
function vector_to_parameters(ps_new::AbstractVector, ps::NamedTuple)
    @assert length(ps_new) == Lux.parameterlength(ps)
    i = 1
    function get_ps(x)
        z = reshape(view(ps_new, i:(i + length(x) - 1)), size(x))
        i += length(x)
        return z
    end
    return fmap(get_ps, ps)
end

## Specify the probabilistic model.
@model function bayes_nn(xs, ts)
    global st

    ## Sample the parameters
    nparameters = Lux.parameterlength(nn)
    parameters ~ MvNormal(zeros(nparameters), sig .* ones(nparameters))

    ## Forward NN to make predictions
    preds, st = nn(xs, vector_to_parameters(parameters, ps), st)

    ## Observe each prediction.
    for i in 1:length(ts)
        ts[i] ~ Bernoulli(preds[i])
    end
end

# Inference can now be performed by calling sample. We use the HMC sampler here.

## Perform inference.
N = 5000
ch = sample(bayes_nn(hcat(xs...), ts), HMC(0.05, 4), N)

# Now we extract the parameter samples from the sampled chain as theta (this is of size `5000 x 20` where `5000` is the number of iterations and `20` is the number of parameters). We'll use these primarily to determine how good our model's classifier is.

## Extract all weight and bias parameters.
theta = MCMCChains.group(ch, :parameters).value;

# ## Prediction Visualization

## A helper to run the nn through data `x` using parameters `theta`
nn_forward(x, theta) = nn(x, vector_to_parameters(theta, ps), st)[1]

## Plot the data we have.
plot_data()

## Find the index that provided the highest log posterior in the chain.
_, i = findmax(ch[:lp])

## Extract the max row value from i.
i = i.I[1]

## Plot the posterior distribution with a contour plot
x1_range = collect(range(-6; stop=6, length=25))
x2_range = collect(range(-6; stop=6, length=25))
Z = [nn_forward([x1, x2], theta[i, :])[1] for x1 in x1_range, x2 in x2_range]
contour!(x1_range, x2_range, Z)

# The contour plot above shows that the MAP method is not too bad at classifying our data. Now we can visualize our predictions.

# $p(\tilde{x} | X, \alpha) = \int_{\theta} p(\tilde{x} | \theta) p(\theta | X, \alpha) \approx \sum_{\theta \sim p(\theta | X, \alpha)}f_{\theta}(\tilde{x})$

# The `nn_predict` function takes the average predicted value from a network parameterized by weights drawn from the MCMC chain.

## Return the average predicted value across multiple weights.
function nn_predict(x, theta, num)
    return mean([nn_forward(x, view(theta, i, :))[1] for i in 1:10:num])
end

# Next, we use the `nn_predict` function to predict the value at a sample of points where the x1 and x2 coordinates range between -6 and 6. As we can see below, we still have a satisfactory fit to our data, and more importantly, we can also see where the neural network is uncertain about its predictions much easier---those regions between cluster boundaries.

# Plot the average prediction.
plot_data()

n_end = 1500
x1_range = collect(range(-6; stop=6, length=25))
x2_range = collect(range(-6; stop=6, length=25))
Z = [nn_predict([x1, x2], theta, n_end)[1] for x1 in x1_range, x2 in x2_range]
contour!(x1_range, x2_range, Z)

# Suppose we are interested in how the predictive power of our Bayesian neural network evolved between samples. In that case, the following graph displays an animation of the contour plot generated from the network weights in samples 1 to 1,000.

## Number of iterations to plot.
n_end = 1000

anim = @gif for i in 1:n_end
    plot_data()
    Z = [nn_forward([x1, x2], theta[i, :])[1] for x1 in x1_range, x2 in x2_range]
    contour!(x1_range, x2_range, Z; title="Iteration $i", clim=(0, 1))
end every 5

ext: jl


text:
module PLaplacianTests

using Test
using Gridap
import Gridap: ∇
using LinearAlgebra

domain = (0,1,0,1)
cells = (4,4)
model = CartesianDiscreteModel(domain,cells)

u(x) = x[1] + x[2]
∇u(x) = VectorValue(1,1)
f(x) = 0
∇(::typeof(u)) = ∇u

const p = 3
flux(∇u) = norm(∇u)^(p-2) * ∇u
dflux(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du

order = 3
T = Float64

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),dirichlet_tags="boundary")
U = TrialFESpace(V,u)

Ω = Triangulation(model)

degree = 2*order
dΩ = Measure(Ω,degree)

res(u,v) = ∫( ∇(v)⋅(flux∘∇(u)) )*dΩ
jac(u,du,v) = ∫( ∇(v)⋅(dflux∘(∇(du),∇(u))) )*dΩ

op = FEOperator(res,jac,U,V)

nls = NLSolver(show_trace=false, method=:newton)
solver = FESolver(nls)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)

e = u - uh

l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

# now with autodiff

op = FEOperator(res,U,V)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)

e = u - uh

l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

end # module


text:
abstract type AbstractNavMessageFile end

struct GPSNavMessageFile <: AbstractNavMessageFile
    prn::Int
    toc::DateTime
    clock_bias::Float64
    clock_drift::Float64
    clock_drift_rate::Float64
    IODE::Float64
    C_rs::Float64
    Δn::Float64
    M_0::Float64
    C_uc::Float64
    e::Float64
    Cus::Float64
    sqrt_A::Float64
    t_oe::Float64
    C_ic::Float64
    ω₀::Float64
    C_is::Float64
    i₀::Float64
    C_rc::Float64
    ω::Float64
    Ω_dot::Float64
    IDOT::Float64
    l2codes::Float64
    gps_week::Float64
    l2pflag::Float64
    ura::Float64
    sv_health::Float64
    T_GD::Float64
    IODC::Float64
    ttm::Float64
    fit_interval::Float64
end


text:
module DiffEqPy

using Reexport
@reexport using DifferentialEquations
using PyCall
pysolve = PyCall.pyfunctionret(solve,Any,Vararg{PyCall.PyAny})

const inspect = PyNULL()
function __init__()
    copy!(inspect, pyimport("inspect"))
end

export pysolve

function DiffEqBase.numargs(f::PyCall.PyObject)
    haskey(f,:py_func) ? _f = f[:py_func] : _f = f
    if PyCall.pyversion < v"3.0.0"
        return length(first(inspect[:getargspec](_f)))
    else
        return length(first(inspect[:getfullargspec](_f)))
    end
end

end # module


text:
"""
    struct ExtractEmpty end

Concrete type to dispatch on for extraction of empty samples.
"""
struct ExtractEmpty end

"""
    extractempty

A singleton of type [`ExtractEmpty`](@ref) is used to signal
downstream extractors that they should extract an empty sample.
"""
const extractempty = ExtractEmpty()

const MissingOrNothing = Union{Missing, Nothing}
const HierarchicType = Union{AbstractDict, AbstractVector, StringOrNumber, MissingOrNothing, ExtractEmpty}

_make_array_node(x, v, store_input) = store_input ? ArrayNode(x, v) : ArrayNode(x)
_make_bag_node(x, bags, v, store_input) = store_input ? BagNode(x, bags, v) : BagNode(x, bags)

"""
returns empty bag of 0 observations
"""
make_empty_bag(x, v) = BagNode(x, Mill.AlignedBags(Vector{UnitRange{Int64}}()))

include("auxiliary.jl")
include("extractarray.jl")
include("extractdict.jl")
include("extractcategorical.jl")
include("extractscalar.jl")
include("extractstring.jl")
include("extractvector.jl")
include("extract_keyasfield.jl")
include("multirepresentation.jl")


text:
@testset "oracle" begin
    data, labels = OneClassActiveLearning.load_data(TEST_DATA_FILE)
    init_strategy = SimpleCombinedStrategy(FixedGammaStrategy(GaussianKernel(2.0)), FixedCStrategy(0.5))
    expected_feedback = Set(OneClassActiveLearning.LABEL_ENCODING.label)

    @testset "initialize" begin
        @test_throws ErrorException OneClassActiveLearning.initialize_oracle(OneClassActiveLearning, data, labels)
        @test_throws ArgumentError OneClassActiveLearning.initialize_oracle(QuerySynthesisOCCOracle, data, labels)
        @test_throws ArgumentError OneClassActiveLearning.initialize_oracle(QuerySynthesisGMMOracle, data, labels)
    end

    @testset "PoolOracle" begin
        oracle = OneClassActiveLearning.initialize_oracle(PoolOracle, data, labels)
        @test isa(oracle, PoolOracle)
        @test ask_oracle(oracle, Array(1:3)) == labels[1:3]
    end

    @testset "NoisyOracle" begin
        oracle = OneClassActiveLearning.initialize_oracle(NoisyOracle, data, labels, Dict{Symbol, Any}(
            :p => 0.5,
            :subtype => :PoolOracle,
            :subparams => Dict{Symbol, Any}()
        ))
        query_ids = Array(1:3)
        feedback = ask_oracle(oracle, query_ids)
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
        oracle = OneClassActiveLearning.initialize_oracle(NoisyOracle, data, labels, Dict{Symbol, Any}(
            :p => 1,
            :subtype => :PoolOracle
        ))
        feedback = ask_oracle(oracle, query_ids)
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
        @test all(feedback .!== labels[query_ids])
    end

    @testset "QuerySynthesisFunctionOracle" begin
        oracle = QuerySynthesisFunctionOracle(x -> fill(:inlier, size(x, 2)))
        @test ask_oracle(oracle, ones(2, 1)) == [:inlier]
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end

    @testset "QuerySynthesisKNNOracle" begin
        @test_throws ArgumentError OneClassActiveLearning.initialize_oracle(QuerySynthesisKNNOracle, data, labels, Dict{Symbol, Any}(:k => 2))
        oracle = OneClassActiveLearning.initialize_oracle(QuerySynthesisKNNOracle, data, labels)
        @test ask_oracle(oracle, data[:, 1:1]) == [labels[1]]
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end

    @testset "QuerySynthesisGMMOracle" begin
        gmm = rand(GMM, 1, 2)
        oracle = QuerySynthesisGMMOracle(gmm, 0.1)
        feedback = ask_oracle(oracle, rand(2, 1))
        @test length(feedback) == 1
        @test issubset(Set(feedback), expected_feedback)
        f = open(TEST_OUTPUT_FILE, "w")
        serialize(f, gmm)
        close(f)
        oracle_param = Dict{Symbol, Any}(:file => TEST_OUTPUT_FILE)
        @test_throws ErrorException OneClassActiveLearning.initialize_oracle(QuerySynthesisGMMOracle, data, labels, oracle_param)
        f = open(TEST_OUTPUT_FILE, "w")
        serialize(f, oracle)
        close(f)
        deserialized_oracle = OneClassActiveLearning.initialize_oracle(QuerySynthesisGMMOracle, data, labels, oracle_param)
        feedback = ask_oracle(deserialized_oracle, rand(2, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end

    @testset "QuerySynthesisOCCOracle" begin
        oracle = OneClassActiveLearning.initialize_oracle(QuerySynthesisOCCOracle, data, labels, Dict{Symbol, Any}(
            :classifier_type => SVDD.RandomOCClassifier,
            :init_strategy => init_strategy,
            :solver => TEST_SOLVER
        ))
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end

    @testset "QuerySynthesisSVMOracle" begin
        oracle = OneClassActiveLearning.initialize_oracle(QuerySynthesisSVMOracle, data, labels, Dict{Symbol, Any}(
            :init_strategy => init_strategy,
        ))
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 1))
        @test length(feedback) == 1
        @test issubset(Set(feedback), expected_feedback)
        oracle = OneClassActiveLearning.initialize_oracle(QuerySynthesisSVMOracle, data, labels, Dict{Symbol, Any}(
            :gamma_search_range_oracle => [0.5, 1.0],
        ))
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end

    @testset "QuerySynthesisCVWrapperOracle" begin
        @test_throws ArgumentError OneClassActiveLearning.initialize_oracle(QuerySynthesisCVWrapperOracle, data, labels)
        @test_throws ArgumentError OneClassActiveLearning.initialize_oracle(QuerySynthesisCVWrapperOracle, data, labels, Dict{Symbol, Any}(
            :subtype => QuerySynthesisCVWrapperOracle,
        ))
        oracle = OneClassActiveLearning.Oracles.initialize_oracle(QuerySynthesisCVWrapperOracle, data, labels, Dict{Symbol, Any}(
            :subtype => QuerySynthesisSVMOracle,
            :gamma_search_range_oracle => [0.1, 1],
            :num_folds => 2,
        ))
        feedback = ask_oracle(oracle, rand(TEST_DATA_NUM_DIMENSIONS, 3))
        @test length(feedback) == 3
        @test issubset(Set(feedback), expected_feedback)
    end
end


text:
@doc raw"""
    JS_BWG_HNE{K}(α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ)
    JS_BWG_HNE{K}(x)

9 parameter model for vertical, northward and eastward displacement of a partical in waves with JONSWAP marginal spectra and bimodal wrapped Gaussian spreading.

# Type parameter
The type parameter `K` is a non-negative integer representing the ammount of aliasing to be done.
This is useful as some devices have high-pass filters which limit the ammount of aliasing that will be present in a recorded series.
In particular, a `JS_BWG_HNE{K}` model for a series recorded every `Δ` seconds would be appropriate for a series sampled every `Δ` seconds and filtered with a high-pass filter at `(2K-1)π/Δ`.

# Arguments
- `α`: The scale parameter.
- `ωₚ`: The peak angular frequency.
- `γ`: The peak enhancement factor.
- `r`: The tail decay index.
- `ϕₘ`: The mean direction.
- `β`: The limiting peak separation.
- `ν`: The peak separation shape.
- `σₗ`: The limiting angular width.
- `σᵣ`: The angular width shape.

# Vector constructor arguments
- `x`: Vector of parameters `[α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ]`.

# Background
The model is a transformation of a model for the frequency-direction spectra of a wave process.
Such a frequency-direction spectra, denoted ``f(ω,ϕ)`` is typically decomposed by writing:
```math
f(ω,ϕ) = f(ω)D(ω,ϕ)
```
where ``f(ω)`` is the marginal spectral density function and ``D(ω,ϕ)`` is the spreading function.
This model uses a JONSWAP marginal sdf and a bimodal wrapped Gaussian spreading function.
The JONSWAP spectral density function is
```math
f(ω) = αω^{-r}\exp\left \{-\frac{r}{4}\left(\frac{|ω|}{ωₚ}\right)^{-4}\right \}γ^{δ(|ω|)}
```
where
```math
δ(ω) = \exp\left\{-\tfrac{1}{2 (0.07+0.02\cdot\mathbb{1}_{ωₚ>|ω|})^2}\left (\tfrac{|ω|}{ωₚ}-1\right )^2\right\}.
```
The bimodal wrapped Gaussian spreading function is
```math
D(ω,ϕ) = \frac{1}{2σ(ω)\sqrt{2π}}\sum\limits_{k=-∞}^{∞} \sum\limits_{i=1}^{2} \exp\left\{-\frac{1}{2}\left(\frac{ϕ-ϕ_{mi}(ω)-2\pi k}{σ(ω)}\right)^2\right\}
```
where
```math
\begin{align*}
    ϕ_{m1}(ω) &= ϕ_m + β \exp\{-ν \min(ω_p/|ω|,1)\}/2 \\
    ϕ_{m2}(ω) &= ϕ_m - β \exp\{-ν \min(ω_p/|ω|,1)\}/2 \\
    σ(ω) &= σ_l - \frac{σ_r}{3}\left\{ 4\left(\frac{ω_p}{|ω|}\right)^2 - \left(\frac{ω_p}{|ω|}\right)^8 \right\}
\end{align*}
```

"""
struct JS_BWG_HNE{K} <: UnknownAcvTimeSeriesModel{3,Float64}
    α::Float64
    ωₚ::Float64
    γ::Float64
    r::Float64
    ϕₘ::Float64
    β::Float64
    ν::Float64
    σₗ::Float64
    σᵣ::Float64

    r_over4::Float64
    ωₚ²::Float64
    ωₚ³::Float64
    ωₚ⁴::Float64
    ωₚ⁶::Float64
    logγ::Float64
    σᵣ_over3::Float64
    cosϕₘ::Float64
    sinϕₘ::Float64
    cos2ϕₘ::Float64
    sin2ϕₘ::Float64
    function JS_BWG_HNE{K}(α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ) where {K}
        α > 0 || throw(ArgumentError("JS_BWG_HNE requires α > 0"))
        ωₚ > 0 || throw(ArgumentError("JS_BWG_HNE requires ωₚ > 0"))
        γ >= 1 || throw(ArgumentError("JS_BWG_HNE requires γ > 1"))
        r > 1 || throw(ArgumentError("JS_BWG_HNE requires r > 1"))
        β >= 0 || throw(ArgumentError("JS_BWG_HNE requires β > 0"))
        ν >= 0 || throw(ArgumentError("JS_BWG_HNE requires ν > 0"))
        σₗ >= 0 || throw(ArgumentError("JS_BWG_HNE requires σₗ > 0"))
        σᵣ >= 0 || throw(ArgumentError("JS_BWG_HNE requires σᵣ > 0"))
        new{K}(α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ,
        r/4,ωₚ^2,ωₚ^3,ωₚ^4,ωₚ^6,log(γ),
        σᵣ/3,cos(ϕₘ),sin(ϕₘ),cos(2ϕₘ),sin(2ϕₘ)
        )
    end
    function JS_BWG_HNE{K}(x::AbstractVector{Float64}) where {K}
        @boundscheck checkparameterlength(x,JS_BWG_HNE{K})
        @inbounds JS_BWG_HNE{K}(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])
    end
end

# functions to throw informative error if type parameter not provided
JS_BWG_HNE(x::AbstractVector{Float64}) = JS_BWG_HNE(ones(9)...)
JS_BWG_HNE(α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ) = error("JS_BWG_HNE process requires the ammount of aliasing specified as a type parameter. Use JS_BWG_HNE{K}() where K ∈ N₀.")

WhittleLikelihoodInference.npars(::Type{JS_BWG_HNE{K}}) where {K} = 9
WhittleLikelihoodInference.nalias(::JS_BWG_HNE{K}) where {K} = K

WhittleLikelihoodInference.lowerbounds(::Type{JS_BWG_HNE{K}}) where {K} = [0,0,1,1,-Inf,0,0,0,0]
WhittleLikelihoodInference.upperbounds(::Type{JS_BWG_HNE{K}}) where {K} = [Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]

@propagate_inbounds @fastmath function WhittleLikelihoodInference.add_sdf!(out, model::JS_BWG_HNE{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:6)
    @inbounds begin
        s_om = sign(ω)
        ω = abs(ω)

        if ω > 1e-10
            α,ωₚ,γ,r,β,ν,σₗ = model.α,model.ωₚ,model.γ,model.r,model.β,model.ν,model.σₗ

            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            ω_over_ωₚ = ω / ωₚ
            ωₚ_over_ω = 1/ω_over_ωₚ
            δ = exp(-1 / (2σ1²) * (ω_over_ωₚ - 1)^2)
            ω⁻⁴ = ω^(-4)
            ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴
            sdf = (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2

            PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ⁴_over_ω⁴)^2)
            σ² = σ^2

            cosPS = cos(PeakSep)
            exp2part = (2 * exp(2.0 * σ²))
            cos2part = (model.cos2ϕₘ*cosPS) / exp2part
            imexppart = sdf * 1im * s_om / exp(0.5 * σ²) * cos(PeakSep/2)

            out[1] += sdf
            out[2] += imexppart * model.cosϕₘ
            out[3] += imexppart * model.sinϕₘ
            out[4] += sdf * (0.5 + cos2part)
            out[5] += sdf * model.sin2ϕₘ * cosPS / exp2part
            out[6] += sdf * (0.5 - cos2part)
        end # add zero otherwise
    end
    return nothing
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.grad_add_sdf!(out, model::JS_BWG_HNE{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:6,1:9)
    @inbounds begin
        s_om = sign(ω)
        ω = abs(ω)

        if ω > 1e-10
            α,ωₚ,γ,r,β,ν,σₗ,σᵣ = model.α,model.ωₚ,model.γ,model.r,model.β,model.ν,model.σₗ,model.σᵣ

            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            ω_over_ωₚ = ω / ωₚ
            ωₚ_over_ω = 1/ω_over_ωₚ
            δ = exp(-1 / (2σ1²) * (ω_over_ωₚ - 1)^2)
            ω⁻⁴ = ω^(-4)
            ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴
            sdf = (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2

            ∂S∂α = sdf / α
            ∂S∂ωₚ = sdf * (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*model.ωₚ³*ω⁻⁴)
            ∂S∂γ = sdf * δ / γ
            ∂S∂r = sdf * (-log(ω)-ωₚ⁴_over_ω⁴/4)

            PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ⁴_over_ω⁴)^2)
            σ² = σ^2

            cosϕₘ = model.cosϕₘ
            sinϕₘ = model.sinϕₘ
            cos2ϕₘ = model.cos2ϕₘ
            sin2ϕₘ = model.sin2ϕₘ
            cosPShalf = cos(PeakSep/2)
            cosPS = cos(PeakSep)
            sinPShalf = sin(PeakSep/2)
            sinPS = sin(PeakSep)

            exp2part = (exp(2.0 * σ²))
            cos2part = (cos2ϕₘ*cosPS) / exp2part
            exphalf = exp(0.5 * σ²)
            imexppart = 1im * s_om / exphalf * cosPShalf

            # hzz = 1
            hxx = (0.5 + cos2part/2)
            hyy = (0.5 - cos2part/2)
            hxz = imexppart * cosϕₘ
            hyz = imexppart * sinϕₘ
            hyx = sin2ϕₘ * cosPS / exp2part / 2

            #### h_derivs
            ## ϕₘ
            hyxpart∂ϕₘ = cosPS / exp2part
            hxx∂ϕₘ = -sin2ϕₘ * hyxpart∂ϕₘ
            hyy∂ϕₘ = -hxx∂ϕₘ
            hzpart∂ϕₘ = -1im * s_om * cosPShalf / exphalf
            hxz∂ϕₘ = sinϕₘ * hzpart∂ϕₘ
            hyz∂ϕₘ = -cosϕₘ * hzpart∂ϕₘ
            hyx∂ϕₘ = cos2ϕₘ * hyxpart∂ϕₘ

            ## β, ν
            PeakSep∂β = PeakSep / β
            PeakSep∂ν = -((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0) * β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            hx∂βνtemp = sinPS / 2 / exp2part
            hz∂βνtemp = -0.5im * sinPShalf / exphalf * s_om

            hx∂βt = hx∂βνtemp * PeakSep∂β
            hx∂νt = hx∂βνtemp * PeakSep∂ν
            hz∂βt = hz∂βνtemp * PeakSep∂β
            hz∂νt = hz∂βνtemp * PeakSep∂ν

            hxx∂β = -cos2ϕₘ * hx∂βt
            hyy∂β = -hxx∂β
            hxz∂β = cosϕₘ * hz∂βt
            hyz∂β = sinϕₘ * hz∂βt
            hyx∂β = -sin2ϕₘ * hx∂βt

            hxx∂ν = -cos2ϕₘ * hx∂νt
            hyy∂ν = -hxx∂ν
            hxz∂ν = cosϕₘ * hz∂νt
            hyz∂ν = sinϕₘ * hz∂νt
            hyx∂ν = -sin2ϕₘ * hx∂νt

            ## σₗ, σᵣ
            # ∂σ∂σₗ = 1
            ∂σ∂σᵣ = (σ-σₗ)/σᵣ
            hx∂σtemp = cosPS * 2 / exp2part * σ
            hz∂σtemp = -1im * cosPShalf / exphalf * σ * s_om

            hx∂σᵣt = hx∂σtemp * ∂σ∂σᵣ
            hz∂σᵣt = hz∂σtemp * ∂σ∂σᵣ

            hxx∂σₗ = -cos2ϕₘ * hx∂σtemp
            hyy∂σₗ = -hxx∂σₗ
            hxz∂σₗ = cosϕₘ * hz∂σtemp
            hyz∂σₗ = sinϕₘ * hz∂σtemp
            hyx∂σₗ = -sin2ϕₘ * hx∂σtemp

            hxx∂σᵣ = -cos2ϕₘ * hx∂σᵣt
            hyy∂σᵣ = -hxx∂σᵣ
            hxz∂σᵣ = cosϕₘ * hz∂σᵣt
            hyz∂σᵣ = sinϕₘ * hz∂σᵣt
            hyx∂σᵣ = -sin2ϕₘ * hx∂σᵣt

            ## ωₚ
            PeakSep∂ωₚ = -(ω > ωₚ) * ν/ω * β * exp(-ν*ωₚ/ω)
            ∂σ∂ωₚ = (-8σᵣ/3*(ωₚ/(ω^2)-ωₚ^7/(ω^8)))

            hx∂ωₚt = (sinPS / 2 * PeakSep∂ωₚ + 2*cosPS*σ*∂σ∂ωₚ)/ exp2part
            hz∂ωₚt = -1im * (sinPShalf / 2 * PeakSep∂ωₚ + cosPShalf*σ*∂σ∂ωₚ)/ exphalf * s_om

            hxx∂ωₚ = -cos2ϕₘ * hx∂ωₚt
            hyy∂ωₚ = -hxx∂ωₚ
            hxz∂ωₚ = cosϕₘ * hz∂ωₚt
            hyz∂ωₚ = sinϕₘ * hz∂ωₚt
            hyx∂ωₚ = -sin2ϕₘ * hx∂ωₚt

            ## α
            out[1, 1] += ∂S∂α
            out[2, 1] += ∂S∂α * hxz
            out[3, 1] += ∂S∂α * hyz
            out[5, 1] += ∂S∂α * hyx
            out[4, 1] += ∂S∂α * hxx
            out[6, 1] += ∂S∂α * hyy

            ## ωₚ
            out[1, 2] += ∂S∂ωₚ
            out[2, 2] += ∂S∂ωₚ * hxz + sdf * hxz∂ωₚ
            out[3, 2] += ∂S∂ωₚ * hyz + sdf * hyz∂ωₚ
            out[4, 2] += ∂S∂ωₚ * hxx + sdf * hxx∂ωₚ
            out[5, 2] += ∂S∂ωₚ * hyx + sdf * hyx∂ωₚ
            out[6, 2] += ∂S∂ωₚ * hyy + sdf * hyy∂ωₚ

            ## γ
            out[1, 3] += ∂S∂γ
            out[2, 3] += ∂S∂γ * hxz
            out[3, 3] += ∂S∂γ * hyz
            out[4, 3] += ∂S∂γ * hxx
            out[5, 3] += ∂S∂γ * hyx
            out[6, 3] += ∂S∂γ * hyy

            ## r
            out[1, 4] += ∂S∂r
            out[2, 4] += ∂S∂r * hxz
            out[3, 4] += ∂S∂r * hyz
            out[4, 4] += ∂S∂r * hxx
            out[5, 4] += ∂S∂r * hyx
            out[6, 4] += ∂S∂r * hyy

            ## ϕₘ
            # out[1, 5] += 0
            out[2, 5] += sdf * hxz∂ϕₘ
            out[3, 5] += sdf * hyz∂ϕₘ
            out[4, 5] += sdf * hxx∂ϕₘ
            out[5, 5] += sdf * hyx∂ϕₘ
            out[6, 5] += sdf * hyy∂ϕₘ

            ## β
            # out[1, 6] += 0
            out[2, 6] += sdf * hxz∂β
            out[3, 6] += sdf * hyz∂β
            out[4, 6] += sdf * hxx∂β
            out[5, 6] += sdf * hyx∂β
            out[6, 6] += sdf * hyy∂β

            ## ν
            # out[1, 7] += 0
            out[2, 7] += sdf * hxz∂ν
            out[3, 7] += sdf * hyz∂ν
            out[4, 7] += sdf * hxx∂ν
            out[5, 7] += sdf * hyx∂ν
            out[6, 7] += sdf * hyy∂ν

            ## σₗ
            # out[1, 8] += 0
            out[2, 8] += sdf * hxz∂σₗ
            out[3, 8] += sdf * hyz∂σₗ
            out[4, 8] += sdf * hxx∂σₗ
            out[5, 8] += sdf * hyx∂σₗ
            out[6, 8] += sdf * hyy∂σₗ

            ## σᵣ
            # out[1, 9] += 0
            out[2, 9] += sdf * hxz∂σᵣ
            out[3, 9] += sdf * hyz∂σᵣ
            out[4, 9] += sdf * hxx∂σᵣ
            out[5, 9] += sdf * hyx∂σᵣ
            out[6, 9] += sdf * hyy∂σᵣ

        end # add zero otherwise
    end

    return nothing
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.hess_add_sdf!(out, model::JS_BWG_HNE{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:6,1:45)
    @inbounds begin
        s_om = sign(ω)
        ω = abs(ω)

        if ω > 1e-10
            α,ωₚ,γ,r,β,ν,σₗ,σᵣ = model.α,model.ωₚ,model.γ,model.r,model.β,model.ν,model.σₗ,model.σᵣ

            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            ω_over_ωₚ = ω / ωₚ
            ωₚ_over_ω = 1/ω_over_ωₚ
            δ = exp(-1 / (2σ1²) * (ω_over_ωₚ - 1)^2)
            ω⁻⁴ = ω^(-4)
            ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴
            sdf = (α * ω^(-r) * exp(-(model.r_over4) * ωₚ⁴_over_ω⁴) * γ^δ) / 2

            ∂S∂α = sdf / α
            ∂S∂ωₚ = sdf * (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*model.ωₚ³*ω⁻⁴)
            ∂S∂γ = sdf * δ / γ
            ∂S∂r = sdf * (-log(ω)-ωₚ⁴_over_ω⁴/4)

            PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ⁴_over_ω⁴)^2)
            σ² = σ^2

            cosϕₘ = model.cosϕₘ
            sinϕₘ = model.sinϕₘ
            cos2ϕₘ = model.cos2ϕₘ
            sin2ϕₘ = model.sin2ϕₘ
            cosPShalf = cos(PeakSep/2)
            cosPS = cos(PeakSep)
            sinPShalf = sin(PeakSep/2)
            sinPS = sin(PeakSep)

            exp2part = (exp(2.0 * σ²))
            cos2part = (cos2ϕₘ*cosPS) / exp2part
            exphalf = exp(0.5 * σ²)
            imexppart = 1im * s_om / exphalf * cosPShalf

            # hzz = 1
            hxx = (0.5 + cos2part/2)
            hyy = (0.5 - cos2part/2)
            hxz = imexppart * cosϕₘ
            hyz = imexppart * sinϕₘ
            hyx = sin2ϕₘ * cosPS / exp2part / 2

            #### h_derivs
            ## ϕₘ
            hyxpart∂ϕₘ = cosPS / exp2part
            hxx∂ϕₘ = -sin2ϕₘ * hyxpart∂ϕₘ
            hyy∂ϕₘ = -hxx∂ϕₘ
            hzpart∂ϕₘ = -1im * s_om * cosPShalf / exphalf
            hxz∂ϕₘ = sinϕₘ * hzpart∂ϕₘ
            hyz∂ϕₘ = -cosϕₘ * hzpart∂ϕₘ
            hyx∂ϕₘ = cos2ϕₘ * hyxpart∂ϕₘ

            ## β, ν
            PeakSep∂β = PeakSep / β
            PeakSep∂ν = -((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0) * β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            hx∂βνtemp = sinPS / 2 / exp2part
            hz∂βνtemp = -0.5im * sinPShalf / exphalf * s_om

            hx∂βt = hx∂βνtemp * PeakSep∂β
            hx∂νt = hx∂βνtemp * PeakSep∂ν
            hz∂βt = hz∂βνtemp * PeakSep∂β
            hz∂νt = hz∂βνtemp * PeakSep∂ν

            hxx∂β = -cos2ϕₘ * hx∂βt
            hyy∂β = -hxx∂β
            hxz∂β = cosϕₘ * hz∂βt
            hyz∂β = sinϕₘ * hz∂βt
            hyx∂β = -sin2ϕₘ * hx∂βt

            hxx∂ν = -cos2ϕₘ * hx∂νt
            hyy∂ν = -hxx∂ν
            hxz∂ν = cosϕₘ * hz∂νt
            hyz∂ν = sinϕₘ * hz∂νt
            hyx∂ν = -sin2ϕₘ * hx∂νt

            ## σₗ, σᵣ
            # ∂σ∂σₗ = 1
            ∂σ∂σᵣ = (σ-σₗ)/σᵣ
            hx∂σtemp = cosPS * 2 / exp2part * σ
            hz∂σtemp = -1im * cosPShalf / exphalf * σ * s_om

            hx∂σᵣt = hx∂σtemp * ∂σ∂σᵣ
            hz∂σᵣt = hz∂σtemp * ∂σ∂σᵣ

            hxx∂σₗ = -cos2ϕₘ * hx∂σtemp
            hyy∂σₗ = -hxx∂σₗ
            hxz∂σₗ = cosϕₘ * hz∂σtemp
            hyz∂σₗ = sinϕₘ * hz∂σtemp
            hyx∂σₗ = -sin2ϕₘ * hx∂σtemp

            hxx∂σᵣ = -cos2ϕₘ * hx∂σᵣt
            hyy∂σᵣ = -hxx∂σᵣ
            hxz∂σᵣ = cosϕₘ * hz∂σᵣt
            hyz∂σᵣ = sinϕₘ * hz∂σᵣt
            hyx∂σᵣ = -sin2ϕₘ * hx∂σᵣt

            ## ωₚ
            PeakSep∂ωₚ = -(ω > ωₚ) * ν/ω * β * exp(-ν*ωₚ/ω)
            ∂σ∂ωₚ = (-8σᵣ/3*(ωₚ/(ω^2)-ωₚ^7/(ω^8)))

            hx∂ωₚt = (sinPS / 2 * PeakSep∂ωₚ + 2*cosPS*σ*∂σ∂ωₚ)/ exp2part
            hz∂ωₚt = -1im * (sinPShalf / 2 * PeakSep∂ωₚ + cosPShalf*σ*∂σ∂ωₚ)/ exphalf * s_om

            hxx∂ωₚ = -cos2ϕₘ * hx∂ωₚt
            hyy∂ωₚ = -hxx∂ωₚ
            hxz∂ωₚ = cosϕₘ * hz∂ωₚt
            hyz∂ωₚ = sinϕₘ * hz∂ωₚt
            hyx∂ωₚ = -sin2ϕₘ * hx∂ωₚt

            ## Peak seperation derivatives
            PeakSep∂β∂ν = PeakSep∂ν/β
            PeakSep∂β∂ωₚ = PeakSep∂ωₚ/β
            PeakSep∂β2 = 0
            PeakSep∂ν∂ωₚ = -(ω > ωₚ) * (ωₚ_over_ω * PeakSep∂ωₚ + PeakSep/ω)
            PeakSep∂ν2 = (((ωₚ_over_ω)^2 - 1) * (ω > ωₚ) + 1) * PeakSep
            PeakSep∂ωₚ2 = (ν/ω) ^ 2 * PeakSep * (ω > ωₚ)
            ## sigma derivatives
            ## most of these are 0, except
            ∂σ∂σᵣ∂ωₚ = ∂σ∂ωₚ / σᵣ
            ∂σ∂ωₚ2 = (-8σᵣ/3*(1/(ω^2)-7*(model.ωₚ⁶)/(ω^8)))

            ## g1
            g1 = -0.5 / exp2part * sinPS
            g1_β = g1 * PeakSep∂β
            g1_ν = g1 * PeakSep∂ν
            g1_ωₚ = g1 * PeakSep∂ωₚ

            minus_exp2part_over2 = -1/exp2part/2
            g1_β2 = minus_exp2part_over2 * PeakSep∂β^2 * cosPS
            g1_βν = minus_exp2part_over2 * (PeakSep∂β∂ν * sinPS + PeakSep∂β * PeakSep∂ν * cosPS)
            g1_βσₗ = 2/exp2part * sinPS * PeakSep∂β * σ
            g1_βσᵣ = g1_βσₗ * ∂σ∂σᵣ

            g1_ν2 = minus_exp2part_over2 * (PeakSep∂ν2 * sinPS + PeakSep∂ν^2 * cosPS)
            g1_νσₗ = 2/exp2part * sinPS * PeakSep∂ν * σ
            g1_νσᵣ = g1_νσₗ * ∂σ∂σᵣ

            g1_ωₚ2 = 1/exp2part * (-sinPS/2 * PeakSep∂ωₚ2 + PeakSep∂ωₚ * (2σ * ∂σ∂ωₚ * sinPS - 0.5 * PeakSep∂ωₚ * cosPS))

            g1_ωₚβ = minus_exp2part_over2 * (PeakSep∂β∂ωₚ * sinPS + PeakSep∂β * PeakSep∂ωₚ * cosPS)
            g1_ωₚν = minus_exp2part_over2 * (PeakSep∂ν∂ωₚ * sinPS + PeakSep∂ωₚ * PeakSep∂ν * cosPS)
            g1_ωₚσₗ = 2/exp2part * sinPS * PeakSep∂ωₚ * σ
            g1_ωₚσᵣ = g1_ωₚσₗ * ∂σ∂σᵣ

            ## g2
            g2_σₗ = -2 / exp2part * cosPS * σ
            g2_σᵣ = g2_σₗ * ∂σ∂σᵣ
            g2_ωₚ = g2_σₗ * ∂σ∂ωₚ

            one_minus_4σ² = (1 - 4σ²)
            g2_σₗ2 = -2 * cosPS / exp2part * one_minus_4σ²
            g2_σₗσᵣ = g2_σₗ2 * ∂σ∂σᵣ
            g2_σᵣ2 = g2_σₗσᵣ * ∂σ∂σᵣ
            g2_ωₚ2 = 2/exp2part * (PeakSep∂ωₚ * sinPS * ∂σ∂ωₚ * σ - cosPS * (∂σ∂ωₚ2 * σ + ∂σ∂ωₚ^2 * one_minus_4σ²))

            g2_ωₚβ = 2/exp2part * sinPS * PeakSep∂β * ∂σ∂ωₚ * σ
            g2_ωₚν = 2/exp2part * sinPS * PeakSep∂ν * ∂σ∂ωₚ * σ
            g2_ωₚσₗ = -2 * cosPS / exp2part * (∂σ∂ωₚ * one_minus_4σ²)
            g2_ωₚσᵣ = -2 * cosPS / exp2part * (σ * ∂σ∂σᵣ∂ωₚ + ∂σ∂σᵣ * ∂σ∂ωₚ * one_minus_4σ²)

            ## g3
            g3 = s_om * 0.5im / exphalf * sinPShalf
            g3_β = g3 * PeakSep∂β
            g3_ν = g3 * PeakSep∂ν
            g3_ωₚ = g3 * PeakSep∂ωₚ

            g3_β2 = s_om * 0.5im / exphalf * (PeakSep∂β2 * sinPShalf + PeakSep∂β^2 * cosPShalf / 2)
            g3_βν = s_om * 0.5im / exphalf * (PeakSep∂β∂ν * sinPShalf + PeakSep∂β * PeakSep∂ν * cosPShalf / 2)
            g3_ν2 = s_om * 0.5im / exphalf * (PeakSep∂ν2 * sinPShalf + PeakSep∂ν^2 * cosPShalf / 2)
            g3_ωₚ2 = s_om * 0.5im / exphalf * ((PeakSep∂ωₚ2 - σ * ∂σ∂ωₚ * PeakSep∂ωₚ) * sinPShalf + PeakSep∂ωₚ^2*cosPShalf/2)
            g3_βσₗ = -s_om * 0.5im / exphalf * sinPShalf * PeakSep∂β * σ
            g3_βσᵣ = g3_βσₗ * ∂σ∂σᵣ
            g3_νσₗ = -s_om * 0.5im / exphalf * sinPShalf * PeakSep∂ν * σ
            g3_νσᵣ = g3_νσₗ * ∂σ∂σᵣ

            g3_ωₚβ = s_om * 1im/exphalf/2 * (PeakSep∂β∂ωₚ * sinPShalf + PeakSep∂β * PeakSep∂ωₚ * cosPShalf/2)
            g3_ωₚν = s_om * 1im/exphalf/2 * (PeakSep∂ν∂ωₚ * sinPShalf + PeakSep∂ωₚ * PeakSep∂ν * cosPShalf/2)
            g3_ωₚσₗ = -s_om * 1im/exphalf/2 * sinPShalf * PeakSep∂ωₚ * σ
            g3_ωₚσᵣ = g3_ωₚσₗ * ∂σ∂σᵣ
            ## g4
            g4 = s_om * 1im / exphalf * cosPShalf * σ
            g4_σₗ = g4
            g4_σᵣ = g4 * ∂σ∂σᵣ
            g4_ωₚ = g4 * ∂σ∂ωₚ

            g4_σₗ2 = s_om * 1im * cosPShalf / exphalf * (1-σ²)
            g4_σₗσᵣ = g4_σₗ2 * ∂σ∂σᵣ
            g4_σᵣ2 = g4_σₗσᵣ * ∂σ∂σᵣ
            g4_ωₚ2 = -s_om * 1im / exphalf * (PeakSep∂ωₚ * ∂σ∂ωₚ * sinPShalf * σ/2 - cosPShalf * (∂σ∂ωₚ2*σ + ∂σ∂ωₚ^2 * (1 - σ²)) )

            g4_ωₚβ = -s_om * 0.5im/exphalf * sinPShalf * PeakSep∂β * ∂σ∂ωₚ * σ
            g4_ωₚν = -s_om * 0.5im/exphalf * sinPShalf * PeakSep∂ν * ∂σ∂ωₚ * σ
            g4_ωₚσₗ = s_om * 1im * cosPShalf / exphalf * (∂σ∂ωₚ * (1 - σ²))
            g4_ωₚσᵣ = s_om * 1im * cosPShalf / exphalf * (σ * ∂σ∂σᵣ∂ωₚ + ∂σ∂σᵣ * ∂σ∂ωₚ * (1 - σ²))

            ## finished required g functions

            #### hxx

            ## ϕₘ (6)
            hxx∂ϕₘbit = -2sin2ϕₘ
            hxx∂ϕₘ2 = -2cos2ϕₘ * hyxpart∂ϕₘ
            hxx∂ϕₘ∂β = hxx∂ϕₘbit * g1_β
            hxx∂ϕₘ∂ν = hxx∂ϕₘbit * g1_ν
            hxx∂ϕₘ∂σₗ = hxx∂ϕₘbit * g2_σₗ
            hxx∂ϕₘ∂σᵣ = hxx∂ϕₘbit * g2_σᵣ
            hxx∂ϕₘ∂ωₚ = hxx∂ϕₘbit * (g1_ωₚ + g2_ωₚ)
            ## ωₚ (5)
            hxx∂ωₚ2 = cos2ϕₘ * (g1_ωₚ2 + g2_ωₚ2)
            hxx∂ωₚ∂β = cos2ϕₘ * (g1_ωₚβ + g2_ωₚβ)
            hxx∂ωₚ∂ν = cos2ϕₘ * (g1_ωₚν + g2_ωₚν)
            hxx∂ωₚ∂σₗ = cos2ϕₘ * (g1_ωₚσₗ + g2_ωₚσₗ)
            hxx∂ωₚ∂σᵣ = cos2ϕₘ * (g1_ωₚσᵣ + g2_ωₚσᵣ)

            ## β (4)
            hxx∂β2 = cos2ϕₘ * g1_β2
            hxx∂β∂ν = cos2ϕₘ * g1_βν
            hxx∂β∂σₗ = cos2ϕₘ * g1_βσₗ
            hxx∂β∂σᵣ = cos2ϕₘ * g1_βσᵣ
            ## ν (3)
            hxx∂ν2 = cos2ϕₘ * g1_ν2
            hxx∂ν∂σₗ = cos2ϕₘ * g1_νσₗ
            hxx∂ν∂σᵣ = cos2ϕₘ * g1_νσᵣ
            ## σₗ (2)
            hxx∂σₗ2 = cos2ϕₘ * g2_σₗ2
            hxx∂σₗ∂σᵣ = cos2ϕₘ * g2_σₗσᵣ
            ## σᵣ (1)
            hxx∂σᵣ2 = cos2ϕₘ * g2_σᵣ2

            #### hyx

            ## ϕₘ (6)
            hyx∂ϕₘ2 = -2sin2ϕₘ * hyxpart∂ϕₘ
            hyx∂ϕₘ∂β = 2cos2ϕₘ * g1_β
            hyx∂ϕₘ∂ν = 2cos2ϕₘ * g1_ν
            hyx∂ϕₘ∂σₗ = 2cos2ϕₘ * g2_σₗ
            hyx∂ϕₘ∂σᵣ = 2cos2ϕₘ * g2_σᵣ
            hyx∂ϕₘ∂ωₚ = 2cos2ϕₘ * (g1_ωₚ + g2_ωₚ)
            ## ωₚ (5)
            hyx∂ωₚ2 = sin2ϕₘ * (g1_ωₚ2 + g2_ωₚ2)
            hyx∂ωₚ∂β = sin2ϕₘ * (g1_ωₚβ + g2_ωₚβ)
            hyx∂ωₚ∂ν = sin2ϕₘ * (g1_ωₚν + g2_ωₚν)
            hyx∂ωₚ∂σₗ = sin2ϕₘ * (g1_ωₚσₗ + g2_ωₚσₗ)
            hyx∂ωₚ∂σᵣ = sin2ϕₘ * (g1_ωₚσᵣ + g2_ωₚσᵣ)

            ## β (4)
            hyx∂β2 = sin2ϕₘ * g1_β2
            hyx∂β∂ν = sin2ϕₘ * g1_βν
            hyx∂β∂σₗ = sin2ϕₘ * g1_βσₗ
            hyx∂β∂σᵣ = sin2ϕₘ * g1_βσᵣ
            ## ν (3)
            hyx∂ν2 = sin2ϕₘ * g1_ν2
            hyx∂ν∂σₗ = sin2ϕₘ * g1_νσₗ
            hyx∂ν∂σᵣ = sin2ϕₘ * g1_νσᵣ
            ## σₗ (2)
            hyx∂σₗ2 = sin2ϕₘ * g2_σₗ2
            hyx∂σₗ∂σᵣ = sin2ϕₘ * g2_σₗσᵣ
            ## σᵣ (1)
            hyx∂σᵣ2 = sin2ϕₘ * g2_σᵣ2

            #### hxz

            ## ϕₘ (6)
            hxz∂ϕₘ2 = cosϕₘ * hzpart∂ϕₘ # minus already added here
            hxz∂ϕₘbit = sinϕₘ # double negative
            hxz∂ϕₘ∂β = hxz∂ϕₘbit * g3_β
            hxz∂ϕₘ∂ν = hxz∂ϕₘbit * g3_ν
            hxz∂ϕₘ∂σₗ = hxz∂ϕₘbit * g4_σₗ
            hxz∂ϕₘ∂σᵣ = hxz∂ϕₘbit * g4_σᵣ
            hxz∂ϕₘ∂ωₚ = hxz∂ϕₘbit * (g3_ωₚ + g4_ωₚ)
            ## ωₚ (5)
            hxz∂ωₚ2 = -cosϕₘ * (g3_ωₚ2 + g4_ωₚ2)
            hxz∂ωₚ∂β = -cosϕₘ * (g3_ωₚβ + g4_ωₚβ)
            hxz∂ωₚ∂ν = -cosϕₘ * (g3_ωₚν + g4_ωₚν)
            hxz∂ωₚ∂σₗ = -cosϕₘ * (g3_ωₚσₗ + g4_ωₚσₗ)
            hxz∂ωₚ∂σᵣ = -cosϕₘ * (g3_ωₚσᵣ + g4_ωₚσᵣ)

            ## β (4)
            hxz∂β2 = -cosϕₘ * g3_β2
            hxz∂β∂ν = -cosϕₘ * g3_βν
            hxz∂β∂σₗ = -cosϕₘ * g3_βσₗ
            hxz∂β∂σᵣ = -cosϕₘ * g3_βσᵣ
            ## ν (3)
            hxz∂ν2 = -cosϕₘ * g3_ν2
            hxz∂ν∂σₗ = -cosϕₘ * g3_νσₗ
            hxz∂ν∂σᵣ = -cosϕₘ * g3_νσᵣ
            ## σₗ (2)
            hxz∂σₗ2 = -cosϕₘ * g4_σₗ2
            hxz∂σₗ∂σᵣ = -cosϕₘ * g4_σₗσᵣ
            ## σᵣ (1)
            hxz∂σᵣ2 = -cosϕₘ * g4_σᵣ2

            #### hyz

            ## ϕₘ (6)
            hyz∂ϕₘ2 = sinϕₘ * hzpart∂ϕₘ
            hyz∂ϕₘbit = -cosϕₘ  # add in minus for conjugate
            hyz∂ϕₘ∂β = hyz∂ϕₘbit * g3_β
            hyz∂ϕₘ∂ν = hyz∂ϕₘbit * g3_ν
            hyz∂ϕₘ∂σₗ = hyz∂ϕₘbit * g4_σₗ
            hyz∂ϕₘ∂σᵣ = hyz∂ϕₘbit * g4_σᵣ
            hyz∂ϕₘ∂ωₚ = hyz∂ϕₘbit * (g3_ωₚ + g4_ωₚ)
            ## ωₚ (5)
            hyz∂ωₚ2 = -sinϕₘ * (g3_ωₚ2 + g4_ωₚ2)
            hyz∂ωₚ∂β = -sinϕₘ * (g3_ωₚβ + g4_ωₚβ)
            hyz∂ωₚ∂ν = -sinϕₘ * (g3_ωₚν + g4_ωₚν)
            hyz∂ωₚ∂σₗ = -sinϕₘ * (g3_ωₚσₗ + g4_ωₚσₗ)
            hyz∂ωₚ∂σᵣ = -sinϕₘ * (g3_ωₚσᵣ + g4_ωₚσᵣ)

            ## β (4)
            hyz∂β2 = -sinϕₘ * g3_β2
            hyz∂β∂ν = -sinϕₘ * g3_βν
            hyz∂β∂σₗ = -sinϕₘ * g3_βσₗ
            hyz∂β∂σᵣ = -sinϕₘ * g3_βσᵣ
            ## ν (3)
            hyz∂ν2 = -sinϕₘ * g3_ν2
            hyz∂ν∂σₗ = -sinϕₘ * g3_νσₗ
            hyz∂ν∂σᵣ = -sinϕₘ * g3_νσᵣ
            ## σₗ (2)
            hyz∂σₗ2 = -sinϕₘ * g4_σₗ2
            hyz∂σₗ∂σᵣ = -sinϕₘ * g4_σₗσᵣ
            ## σᵣ (1)
            hyz∂σᵣ2 = -sinϕₘ * g4_σᵣ2

            #### out

            ## α2
            # out[1, 1] += 0
            # out[2, 1] += 0
            # out[3, 1] += 0
            # out[4, 1] += 0
            # out[5, 1] += 0
            # out[6, 1] += 0

            ## α ωₚ
            ∂S∂αωₚ = ∂S∂ωₚ / α
            out[1, 2] += ∂S∂αωₚ
            out[2, 2] += ∂S∂αωₚ * hxz + ∂S∂α * hxz∂ωₚ
            out[3, 2] += ∂S∂αωₚ * hyz + ∂S∂α * hyz∂ωₚ
            out[4, 2] += ∂S∂αωₚ * hxx + ∂S∂α * hxx∂ωₚ
            out[5, 2] += ∂S∂αωₚ * hyx + ∂S∂α * hyx∂ωₚ
            out[6, 2] += ∂S∂αωₚ * hyy + ∂S∂α * hyy∂ωₚ

            ## α γ
            ∂S∂αγ = ∂S∂γ / α
            out[1, 3] += ∂S∂αγ
            out[2, 3] += ∂S∂αγ * hxz
            out[3, 3] += ∂S∂αγ * hyz
            out[4, 3] += ∂S∂αγ * hxx
            out[5, 3] += ∂S∂αγ * hyx
            out[6, 3] += ∂S∂αγ * hyy

            ## α r
            ∂S∂αr = ∂S∂r / α
            out[1, 4] += ∂S∂αr
            out[2, 4] += ∂S∂αr * hxz
            out[3, 4] += ∂S∂αr * hyz
            out[4, 4] += ∂S∂αr * hxx
            out[5, 4] += ∂S∂αr * hyx
            out[6, 4] += ∂S∂αr * hyy

            ## α ϕₘ
            # out[1, 5] += 0
            out[2, 5] += ∂S∂α * hxz∂ϕₘ
            out[3, 5] += ∂S∂α * hyz∂ϕₘ
            out[4, 5] += ∂S∂α * hxx∂ϕₘ
            out[5, 5] += ∂S∂α * hyx∂ϕₘ
            out[6, 5] += ∂S∂α * hyy∂ϕₘ

            ## α β
            # out[1, 6] += 0
            out[2, 6] += ∂S∂α * hxz∂β
            out[3, 6] += ∂S∂α * hyz∂β
            out[4, 6] += ∂S∂α * hxx∂β
            out[5, 6] += ∂S∂α * hyx∂β
            out[6, 6] += ∂S∂α * hyy∂β

            ## α ν
            # out[1, 7] += 0
            out[2, 7] += ∂S∂α * hxz∂ν
            out[3, 7] += ∂S∂α * hyz∂ν
            out[4, 7] += ∂S∂α * hxx∂ν
            out[5, 7] += ∂S∂α * hyx∂ν
            out[6, 7] += ∂S∂α * hyy∂ν

            ## α σₗ
            # out[1, 8] += 0
            out[2, 8] += ∂S∂α * hxz∂σₗ
            out[3, 8] += ∂S∂α * hyz∂σₗ
            out[4, 8] += ∂S∂α * hxx∂σₗ
            out[5, 8] += ∂S∂α * hyx∂σₗ
            out[6, 8] += ∂S∂α * hyy∂σₗ

            ## α σᵣ
            # out[1, 9] += 0
            out[2, 9] += ∂S∂α * hxz∂σᵣ
            out[3, 9] += ∂S∂α * hyz∂σᵣ
            out[4, 9] += ∂S∂α * hxx∂σᵣ
            out[5, 9] += ∂S∂α * hyx∂σᵣ
            out[6, 9] += ∂S∂α * hyy∂σᵣ

            ## ωₚ2
            ∂S∂ωₚ2part = δ * log(γ) * ω / σ1²
            ∂S∂ωₚUsepart = (∂S∂ωₚ2part * (ω-ωₚ) / (ωₚ^3)) - (r * (ωₚ^3) / (ω^4))
            ∂S∂ωₚ2 = ∂S∂ωₚ * ∂S∂ωₚUsepart + sdf * (∂S∂ωₚ2part *((-3ω + 2ωₚ)/ωₚ^4 + ω * (ω-ωₚ)^2/ωₚ^6/σ1²) - 3r*ωₚ^2/ω^4)
            out[1, 10] += ∂S∂ωₚ2
            out[2, 10] += ∂S∂ωₚ2 * hxz + 2∂S∂ωₚ * hxz∂ωₚ + sdf*hxz∂ωₚ2
            out[3, 10] += ∂S∂ωₚ2 * hyz + 2∂S∂ωₚ * hyz∂ωₚ + sdf*hyz∂ωₚ2
            out[4, 10] += ∂S∂ωₚ2 * hxx + 2∂S∂ωₚ * hxx∂ωₚ + sdf*hxx∂ωₚ2
            out[5, 10] += ∂S∂ωₚ2 * hyx + 2∂S∂ωₚ * hyx∂ωₚ + sdf*hyx∂ωₚ2
            out[6, 10] += ∂S∂ωₚ2 * hyy + 2∂S∂ωₚ * hyy∂ωₚ - sdf*hxx∂ωₚ2

            ## ωₚγ
            ∂S∂ωₚγ = ∂S∂γ * ∂S∂ωₚUsepart + sdf * δ / γ * ω / σ1²*(ω-ωₚ)/ωₚ^3
            out[1, 11] += ∂S∂ωₚγ
            out[2, 11] += ∂S∂ωₚγ * hxz + ∂S∂γ * hxz∂ωₚ
            out[3, 11] += ∂S∂ωₚγ * hyz + ∂S∂γ * hyz∂ωₚ
            out[4, 11] += ∂S∂ωₚγ * hxx + ∂S∂γ * hxx∂ωₚ
            out[5, 11] += ∂S∂ωₚγ * hyx + ∂S∂γ * hyx∂ωₚ
            out[6, 11] += ∂S∂ωₚγ * hyy + ∂S∂γ * hyy∂ωₚ

            ## ωₚr
            ∂S∂ωₚr = ∂S∂r * ∂S∂ωₚUsepart - sdf * ωₚ^3/ω^4
            out[1, 12] += ∂S∂ωₚr
            out[2, 12] += ∂S∂ωₚr * hxz + ∂S∂r * hxz∂ωₚ
            out[3, 12] += ∂S∂ωₚr * hyz + ∂S∂r * hyz∂ωₚ
            out[4, 12] += ∂S∂ωₚr * hxx + ∂S∂r * hxx∂ωₚ
            out[5, 12] += ∂S∂ωₚr * hyx + ∂S∂r * hyx∂ωₚ
            out[6, 12] += ∂S∂ωₚr * hyy + ∂S∂r * hyy∂ωₚ

            ## ωₚϕₘ
            # out[1, 13] += 0
            out[2, 13] += ∂S∂ωₚ * hxz∂ϕₘ + sdf * hxz∂ϕₘ∂ωₚ
            out[3, 13] += ∂S∂ωₚ * hyz∂ϕₘ + sdf * hyz∂ϕₘ∂ωₚ
            out[4, 13] += ∂S∂ωₚ * hxx∂ϕₘ + sdf * hxx∂ϕₘ∂ωₚ
            out[5, 13] += ∂S∂ωₚ * hyx∂ϕₘ + sdf * hyx∂ϕₘ∂ωₚ
            out[6, 13] += ∂S∂ωₚ * hyy∂ϕₘ - sdf * hxx∂ϕₘ∂ωₚ

            ## ωₚβ
            # out[1, 14] += 0
            out[2, 14] += ∂S∂ωₚ * hxz∂β + sdf * hxz∂ωₚ∂β
            out[3, 14] += ∂S∂ωₚ * hyz∂β + sdf * hyz∂ωₚ∂β
            out[4, 14] += ∂S∂ωₚ * hxx∂β + sdf * hxx∂ωₚ∂β
            out[5, 14] += ∂S∂ωₚ * hyx∂β + sdf * hyx∂ωₚ∂β
            out[6, 14] += ∂S∂ωₚ * hyy∂β - sdf * hxx∂ωₚ∂β

            ## ωₚν
            # out[1, 15] += 0
            out[2, 15] += ∂S∂ωₚ * hxz∂ν + sdf * hxz∂ωₚ∂ν
            out[3, 15] += ∂S∂ωₚ * hyz∂ν + sdf * hyz∂ωₚ∂ν
            out[4, 15] += ∂S∂ωₚ * hxx∂ν + sdf * hxx∂ωₚ∂ν
            out[5, 15] += ∂S∂ωₚ * hyx∂ν + sdf * hyx∂ωₚ∂ν
            out[6, 15] += ∂S∂ωₚ * hyy∂ν - sdf * hxx∂ωₚ∂ν

            ## ωₚσₗ
            # out[1, 16] += 0
            out[2, 16] += ∂S∂ωₚ * hxz∂σₗ + sdf * hxz∂ωₚ∂σₗ
            out[3, 16] += ∂S∂ωₚ * hyz∂σₗ + sdf * hyz∂ωₚ∂σₗ
            out[4, 16] += ∂S∂ωₚ * hxx∂σₗ + sdf * hxx∂ωₚ∂σₗ
            out[5, 16] += ∂S∂ωₚ * hyx∂σₗ + sdf * hyx∂ωₚ∂σₗ
            out[6, 16] += ∂S∂ωₚ * hyy∂σₗ - sdf * hxx∂ωₚ∂σₗ

            ## ωₚσᵣ
            # out[1, 17] += 0
            out[2, 17] += ∂S∂ωₚ * hxz∂σᵣ + sdf * hxz∂ωₚ∂σᵣ
            out[3, 17] += ∂S∂ωₚ * hyz∂σᵣ + sdf * hyz∂ωₚ∂σᵣ
            out[4, 17] += ∂S∂ωₚ * hxx∂σᵣ + sdf * hxx∂ωₚ∂σᵣ
            out[5, 17] += ∂S∂ωₚ * hyx∂σᵣ + sdf * hyx∂ωₚ∂σᵣ
            out[6, 17] += ∂S∂ωₚ * hyy∂σᵣ - sdf * hxx∂ωₚ∂σᵣ

            ## γ2
            ∂S∂γ2 = δ/γ * (∂S∂γ - sdf/γ)
            out[1, 18] += ∂S∂γ2
            out[2, 18] += ∂S∂γ2 * hxz
            out[3, 18] += ∂S∂γ2 * hyz
            out[4, 18] += ∂S∂γ2 * hxx
            out[5, 18] += ∂S∂γ2 * hyx
            out[6, 18] += ∂S∂γ2 * hyy

            ## γr
            ∂S∂γr = δ/γ * ∂S∂r
            out[1, 19] += ∂S∂γr
            out[2, 19] += ∂S∂γr * hxz
            out[3, 19] += ∂S∂γr * hyz
            out[4, 19] += ∂S∂γr * hxx
            out[5, 19] += ∂S∂γr * hyx
            out[6, 19] += ∂S∂γr * hyy

            ## γϕₘ
            # out[1, 20] += 0
            out[2, 20] += ∂S∂γ * hxz∂ϕₘ
            out[3, 20] += ∂S∂γ * hyz∂ϕₘ
            out[4, 20] += ∂S∂γ * hxx∂ϕₘ
            out[5, 20] += ∂S∂γ * hyx∂ϕₘ
            out[6, 20] += ∂S∂γ * hyy∂ϕₘ

            ## γβ
            # out[1, 21] += 0
            out[2, 21] += ∂S∂γ * hxz∂β
            out[3, 21] += ∂S∂γ * hyz∂β
            out[4, 21] += ∂S∂γ * hxx∂β
            out[5, 21] += ∂S∂γ * hyx∂β
            out[6, 21] += ∂S∂γ * hyy∂β

            ## γν
            # out[1, 22] += 0
            out[2, 22] += ∂S∂γ * hxz∂ν
            out[3, 22] += ∂S∂γ * hyz∂ν
            out[4, 22] += ∂S∂γ * hxx∂ν
            out[5, 22] += ∂S∂γ * hyx∂ν
            out[6, 22] += ∂S∂γ * hyy∂ν

            ## γσₗ
            # out[1, 23] += 0
            out[2, 23] += ∂S∂γ * hxz∂σₗ
            out[3, 23] += ∂S∂γ * hyz∂σₗ
            out[4, 23] += ∂S∂γ * hxx∂σₗ
            out[5, 23] += ∂S∂γ * hyx∂σₗ
            out[6, 23] += ∂S∂γ * hyy∂σₗ

            ## γσᵣ
            # out[1, 24] += 0
            out[2, 24] += ∂S∂γ * hxz∂σᵣ
            out[3, 24] += ∂S∂γ * hyz∂σᵣ
            out[4, 24] += ∂S∂γ * hxx∂σᵣ
            out[5, 24] += ∂S∂γ * hyx∂σᵣ
            out[6, 24] += ∂S∂γ * hyy∂σᵣ

            ##  r2
            ∂S∂r2 = (-log(ω) - 0.25*(ωₚ/ω)^4) * ∂S∂r
            out[1, 25] += ∂S∂r2
            out[2, 25] += ∂S∂r2 * hxz
            out[3, 25] += ∂S∂r2 * hyz
            out[4, 25] += ∂S∂r2 * hxx
            out[5, 25] += ∂S∂r2 * hyx
            out[6, 25] += ∂S∂r2 * hyy

            ## rϕₘ
            # out[1, 26] += 0
            out[2, 26] += ∂S∂r * hxz∂ϕₘ
            out[3, 26] += ∂S∂r * hyz∂ϕₘ
            out[4, 26] += ∂S∂r * hxx∂ϕₘ
            out[5, 26] += ∂S∂r * hyx∂ϕₘ
            out[6, 26] += ∂S∂r * hyy∂ϕₘ

            ## rβ
            # out[1, 27] += 0
            out[2, 27] += ∂S∂r * hxz∂β
            out[3, 27] += ∂S∂r * hyz∂β
            out[4, 27] += ∂S∂r * hxx∂β
            out[5, 27] += ∂S∂r * hyx∂β
            out[6, 27] += ∂S∂r * hyy∂β

            ## rν
            # out[1, 28] += 0
            out[2, 28] += ∂S∂r * hxz∂ν
            out[3, 28] += ∂S∂r * hyz∂ν
            out[4, 28] += ∂S∂r * hxx∂ν
            out[5, 28] += ∂S∂r * hyx∂ν
            out[6, 28] += ∂S∂r * hyy∂ν

            ## rσₗ
            # out[1, 29] += 0
            out[2, 29] += ∂S∂r * hxz∂σₗ
            out[3, 29] += ∂S∂r * hyz∂σₗ
            out[4, 29] += ∂S∂r * hxx∂σₗ
            out[5, 29] += ∂S∂r * hyx∂σₗ
            out[6, 29] += ∂S∂r * hyy∂σₗ

            ## rσᵣ
            # out[1, 30] += 0
            out[2, 30] += ∂S∂r * hxz∂σᵣ
            out[3, 30] += ∂S∂r * hyz∂σᵣ
            out[4, 30] += ∂S∂r * hxx∂σᵣ
            out[5, 30] += ∂S∂r * hyx∂σᵣ
            out[6, 30] += ∂S∂r * hyy∂σᵣ

            ## ϕₘ2
            # out[1, 31] += 0
            out[2, 31] += sdf * hxz∂ϕₘ2
            out[3, 31] += sdf * hyz∂ϕₘ2
            out[4, 31] += sdf * hxx∂ϕₘ2
            out[5, 31] += sdf * hyx∂ϕₘ2
            out[6, 31] += sdf * -hxx∂ϕₘ2

            ## ϕₘβ
            # out[1, 32] += 0
            out[2, 32] += sdf * hxz∂ϕₘ∂β
            out[3, 32] += sdf * hyz∂ϕₘ∂β
            out[4, 32] += sdf * hxx∂ϕₘ∂β
            out[5, 32] += sdf * hyx∂ϕₘ∂β
            out[6, 32] += sdf * -hxx∂ϕₘ∂β

            ## ϕₘν
            # out[1, 33] += 0
            out[2, 33] += sdf * hxz∂ϕₘ∂ν
            out[3, 33] += sdf * hyz∂ϕₘ∂ν
            out[4, 33] += sdf * hxx∂ϕₘ∂ν
            out[5, 33] += sdf * hyx∂ϕₘ∂ν
            out[6, 33] += sdf * -hxx∂ϕₘ∂ν

            ## ϕₘσₗ
            # out[1, 34] += 0
            out[2, 34] += sdf * hxz∂ϕₘ∂σₗ
            out[3, 34] += sdf * hyz∂ϕₘ∂σₗ
            out[4, 34] += sdf * hxx∂ϕₘ∂σₗ
            out[5, 34] += sdf * hyx∂ϕₘ∂σₗ
            out[6, 34] += sdf * -hxx∂ϕₘ∂σₗ

            ## ϕₘσᵣ
            # out[1, 35] += 0
            out[2, 35] += sdf * hxz∂ϕₘ∂σᵣ
            out[3, 35] += sdf * hyz∂ϕₘ∂σᵣ
            out[4, 35] += sdf * hxx∂ϕₘ∂σᵣ
            out[5, 35] += sdf * hyx∂ϕₘ∂σᵣ
            out[6, 35] += sdf * -hxx∂ϕₘ∂σᵣ

            ## β2
            # out[1, 36] += 0
            out[2, 36] += sdf * hxz∂β2
            out[3, 36] += sdf * hyz∂β2
            out[4, 36] += sdf * hxx∂β2
            out[5, 36] += sdf * hyx∂β2
            out[6, 36] += sdf * -hxx∂β2

            ## βν
            # out[1, 37] += 0
            out[2, 37] += sdf * hxz∂β∂ν
            out[3, 37] += sdf * hyz∂β∂ν
            out[4, 37] += sdf * hxx∂β∂ν
            out[5, 37] += sdf * hyx∂β∂ν
            out[6, 37] += sdf * -hxx∂β∂ν

            ## βσₗ
            # out[1, 38] += 0
            out[2, 38] += sdf * hxz∂β∂σₗ
            out[3, 38] += sdf * hyz∂β∂σₗ
            out[4, 38] += sdf * hxx∂β∂σₗ
            out[5, 38] += sdf * hyx∂β∂σₗ
            out[6, 38] += sdf * -hxx∂β∂σₗ

            ## βσᵣ
            # out[1, 39] += 0
            out[2, 39] += sdf * hxz∂β∂σᵣ
            out[3, 39] += sdf * hyz∂β∂σᵣ
            out[4, 39] += sdf * hxx∂β∂σᵣ
            out[5, 39] += sdf * hyx∂β∂σᵣ
            out[6, 39] += sdf * -hxx∂β∂σᵣ

            ## ν2
            # out[1, 40] += 0
            out[2, 40] += sdf * hxz∂ν2
            out[3, 40] += sdf * hyz∂ν2
            out[4, 40] += sdf * hxx∂ν2
            out[5, 40] += sdf * hyx∂ν2
            out[6, 40] += sdf * -hxx∂ν2

            ## νσₗ
            # out[1, 41] += 0
            out[2, 41] += sdf * hxz∂ν∂σₗ
            out[3, 41] += sdf * hyz∂ν∂σₗ
            out[4, 41] += sdf * hxx∂ν∂σₗ
            out[5, 41] += sdf * hyx∂ν∂σₗ
            out[6, 41] += sdf * -hxx∂ν∂σₗ

            ## νσᵣ
            # out[1, 42] += 0
            out[2, 42] += sdf * hxz∂ν∂σᵣ
            out[3, 42] += sdf * hyz∂ν∂σᵣ
            out[4, 42] += sdf * hxx∂ν∂σᵣ
            out[5, 42] += sdf * hyx∂ν∂σᵣ
            out[6, 42] += sdf * -hxx∂ν∂σᵣ

            ## σₗ2
            # out[1, 43] += 0
            out[2, 43] += sdf * hxz∂σₗ2
            out[3, 43] += sdf * hyz∂σₗ2
            out[4, 43] += sdf * hxx∂σₗ2
            out[5, 43] += sdf * hyx∂σₗ2
            out[6, 43] += sdf * -hxx∂σₗ2

            ## σₗσᵣ
            # out[1, 44] += 0
            out[2, 44] += sdf * hxz∂σₗ∂σᵣ
            out[3, 44] += sdf * hyz∂σₗ∂σᵣ
            out[4, 44] += sdf * hxx∂σₗ∂σᵣ
            out[5, 44] += sdf * hyx∂σₗ∂σᵣ
            out[6, 44] += sdf * -hxx∂σₗ∂σᵣ

            ## σᵣ2
            # out[1, 45] += 0
            out[2, 45] += sdf * hxz∂σᵣ2
            out[3, 45] += sdf * hyz∂σᵣ2
            out[4, 45] += sdf * hxx∂σᵣ2
            out[5, 45] += sdf * hyx∂σᵣ2
            out[6, 45] += sdf * -hxx∂σᵣ2

        end # otherwise add 0
    end

    return nothing
end

function WhittleLikelihoodInference.coherancy(model::JS_BWG_HNE{K}, ω::Real) where {K}
    s_om = sign(ω)
    ω = abs(ω)
    α,ωₚ,γ,r,ϕₘ,β,ν,σₗ,σᵣ = model.α,model.ωₚ,model.γ,model.r,model.ϕₘ,model.β,model.ν,model.σₗ,model.σᵣ

    ω_over_ωₚ = ω / ωₚ
    ωₚ_over_ω = 1/ω_over_ωₚ
    ω⁻⁴ = ω^(-4)
    ωₚ⁴_over_ω⁴ = model.ωₚ⁴ * ω⁻⁴

    PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
    σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ⁴_over_ω⁴)^2)
    σ² = σ^2

    cosPS = cos(PeakSep)
    exp2part = (2 * exp(2.0 * σ²))
    cos2part = (model.cos2ϕₘ*cosPS) / exp2part
    imexppart = 1im * s_om / exp(0.5 * σ²) * cos(PeakSep/2)
    out = zeros(ComplexF64, 6)

    out[1] = 1
    out[2] = imexppart * model.cosϕₘ / sqrt(0.5 + cos2part)
    out[3] = imexppart * model.sinϕₘ / sqrt(0.5 - cos2part)
    out[4] = 1
    out[5] = model.sin2ϕₘ * cosPS / exp2part / sqrt(0.5 + cos2part) / sqrt(0.5 - cos2part)
    out[6] = 1

    return WhittleLikelihoodInference.SHermitianCompact{3, ComplexF64, 6}(out)
end


text:
using LinearAlgebra: norm
using Printf: @printf

"""
    QuadSS(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate a function over the range (-∞, ∞) using the
*sinh-sinh quadrature*. It utilizes the change of variables to transform the
integrand into a form well-suited to the trapezoidal rule.

`QuadSS` tries to calculate integral values `maxlevel` times at a maximum;
the step size of a trapezoid is started from `h0` and is halved in each
following repetition for finer accuracy. The repetition is terminated when the
difference from the previous estimation gets smaller than a certain threshold.
The threshold is determined by the runtime parameters, see below.

The type `T` represents the accuracy of interval. The integrand should accept
values `x<:T` as its parameter.

---

    I, E = (q::QuadSS)(f::Function;
                       atol::Real=zero(T),
                       rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                       where {T<:AbstractFloat}

Numerically integrate `f(x)` over the interval (-∞, ∞) and return the integral
value `I` and an estimated error `E`. The `E` is not exactly equal to the
difference from the true value. However, one can expect that the integral value
`I` is converged if `E <= max(atol, rtol*norm(I))` is true. Otherwise, the
obtained `I` would be unreliable; the number of repetitions exceeds the
`maxlevel` before converged.

The integrand `f` can also return any value other than a scalar, as far as
`+`, `-`, multiplication by real values, and `LinearAlgebra.norm`, are
implemented. For example, `Vector` or `Array` of numbers are acceptable
although, unfortunately, it may not be very performant.

# Examples
```jldoctest
julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qss = QuadSS(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qss(f);

julia> I ≈ 2π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = qss(x -> [1/(1 + x^2), 2/(1 + x^2)]);

julia> I ≈ [π, 2π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true
```
"""
struct QuadSS{T<:AbstractFloat,N}
    h0::T
    origin::Tuple{T,T}
    table0::Vector{Tuple{T,T}}
    tables::NTuple{N,Vector{Tuple{T,T}}}
end
function QuadSS(T::Type{<:AbstractFloat}; maxlevel::Integer=12, h0::Real=one(T))
    @assert maxlevel > 0
    @assert h0 > 0
    origin = samplepoint(QuadSS, zero(T))
    table0 = generate_table(QuadSS, h0, 1)
    tables = Vector{Tuple{T,T}}[]
    for level in 1:maxlevel
        h = h0/2^level
        table = generate_table(QuadSS, h, 2)
        push!(tables, table)
    end
    return QuadSS{T,maxlevel}(T(h0), origin, table0, Tuple(tables))
end

function (q::QuadSS{T,N})(f::Function; atol::Real=zero(T),
                          rtol::Real=atol>0 ? zero(T) : sqrt(eps(T))) where {T<:AbstractFloat,N}
    f⁺ = f
    f⁻ = u -> f(-u)
    eval⁺(t) = f⁺(t[1])*t[2]
    eval⁻(t) = f⁻(t[1])*t[2]
    x0, w0 = q.origin
    Σ = f(x0)*w0
    istart⁺ = startindex(f⁺, q.table0, 1)
    istart⁻ = startindex(f⁻, q.table0, 1)
    Σ += mapsum(eval⁺, q.table0, istart⁺)
    Σ += mapsum(eval⁻, q.table0, istart⁻)
    h0 = q.h0
    I = h0*Σ
    E = zero(eltype(I))
    for level in 1:N
        table = q.tables[level]
        istart⁺ = startindex(f⁺, table, 2*istart⁺ - 1)
        istart⁻ = startindex(f⁻, table, 2*istart⁻ - 1)
        Σ += mapsum(eval⁺, table, istart⁺)
        Σ += mapsum(eval⁻, table, istart⁻)
        h = h0/2^level
        prevI = I
        I = h*Σ
        E = estimate_error(T, prevI, I)
        tol = max(norm(I)*rtol, atol)
        !(E > tol) && break
    end
    return I, E
end

function Base.show(io::IO, ::MIME"text/plain", q::QuadSS{T,N}) where {T<:AbstractFloat,N}
    @printf(io, "DoubleExponentialFormulas.QuadSS{%s}: maxlevel=%d, h0=%.3e",
            string(T), N, q.h0)
end

function generate_table(::Type{QuadSS}, h::T, step::Int) where {T<:AbstractFloat}
    table = Tuple{T,T}[]
    k = 1
    while true
        t = k*h
        xk, wk = samplepoint(QuadSS, t)
        xk ≥ floatmax(T) && break
        wk ≥ floatmax(T) && break
        push!(table, (xk, wk))
        k += step
    end
    reverse!(table)
    return table
end

function samplepoint(::Type{QuadSS}, t)
    sinht = sinh(t)
    ϕ = sinh(sinht*π/2)
    ϕ′ = (cosh(t)*π/2)*cosh(sinht*π/2)
    return ϕ, ϕ′
end


text:
"""

$(SIGNATURES)

Request a NIPA table from the BEA data API and return an object of type [`BeaNipaTable`](@ref).

Arguments
---------
* `b` -- a [`Bea`](@ref) connection
* `TableID` -- the integer Table ID for the desired NIPA table
* `frequency` -- "A" for annual, "Q" for quarerly
* `startyear` -- first year of data requested, in YYYY format
* `endyear` -- last year of data requested, in YYYY format

"""
function get_nipa_table(b::Bea, TableID::Int, frequency::AbstractString, startyear::Int, endyear::Int)
    url = b.url
    key = b.key
    dataset = b.dataset
    bea_method = "GetData"
    years = collect(startyear:1:endyear)

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "DatasetName" => dataset,
                     "TableID" => TableID,
                     "Frequency" => frequency,
                     "Year" => years,
                     "ResultFormat" => "JSON")

    response = get(url; query = querydict)
    response_json = Requests.json(response)

    # Extract the vector of Dicts containing the data
    data_dicts = response_json["BEAAPI"]["Results"]["Data"]
    # Retrieve data from each Dict
    all_data = [parse_data_dict(d) for d in data_dicts]
    values = [tup[1] for tup in all_data]
    dates = [tup[2] for tup in all_data]
    linenums = [tup[3] for tup in all_data]
    # Create data frame of values, reshape, and create names
    dflong = DataFrame(date = dates, line = linenums, value = values)
    df = unstack(dflong, :date, :line, :value) # Convert to "wide" format
    newnames = [Symbol(string("line", name)) for name in names(df)[2:end]]
    names!(df, [:date; newnames])

    # Create OrderedDict of line descriptions
    linedesc = [tup[4] for tup in all_data]
    linekeys = OrderedDict{AbstractString, AbstractString}(zip(linenums, linedesc))

    # Extract metadata for the table
    tablenotes = response_json["BEAAPI"]["Results"]["Notes"]
    tableinfo = split(tablenotes[1]["NoteText"])
    tablenum = rstrip(tableinfo[2], '.')
    tabledesc = join(tableinfo[3:end], " ")
    if length(tablenotes) > 1
        noteref = [note["NoteRef"] for note in tablenotes[2:end]]
        notetext = [note["NoteText"] for note in tablenotes[2:end]]
        notes = OrderedDict{AbstractString, AbstractString}(zip(noteref, notetext))
    else
        notes = ""
    end

    frequency == "Q"? freq = "Quarterly": freq = "Annual"
    out = BeaNipaTable(tablenum, TableID, tabledesc, linekeys, notes, freq, startyear, endyear, df)
    return out
end

"""
$(SIGNATURES)

Extract information for a single observation and return as a tuple.  (Internal method for [`get_nipa_table`](@ref).)
"""
function parse_data_dict(dict::Dict)
    # This function
    linenum = dict["LineNumber"]
    linedesc = dict["LineDescription"]
    # Add footnote indicator to line description
    if dict["NoteRef"] != dict["TableID"]
        fn = split(dict["NoteRef"], ',')
        linedesc = string(linedesc, " (", fn[2], ")")
    end
    # Remove commas from data values
    dataval = float(replace(dict["DataValue"], ",", ""))

    # Change date from string to Date() type
    timeperiod = dict["TimePeriod"]
    year = parse(Int, timeperiod[1:4])
    # Quarterly data
    if ismatch(r"Q", timeperiod)
        quarter = parse(Int, timeperiod[end])
        date = Date(year, (quarter*3 - 2), 1)
    # Annual data
    else
        date = Date(year, 12, 31)
    end

    return (dataval, date, linenum, linedesc)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

include("traits/variable.jl")
include("traits/domain.jl")
include("traits/data.jl")
include("traits/optional.jl")

# type alias for convenience
const DomainOrData = Union{Domain,Data}


text:
# Hard class group of Carlo
# Qx, x  = PolynomialRing(FlintQQ, "x");
# f = x^4-x^3-378*x^2+2081*x-67;
# K, a = number_field(f);
# C = cyclotomic_extension(K, 7)
# K = absolute_field(C);

f = x^24-2*x^23-2266*x^22+15510*x^21+2122618*x^20-25022546*x^19-994971857*x^18+17780224498*x^17+204666064586*x^16-6330369896072*x^15+6277040017522*x^14+1024499591458410*x^13-9491657191588213*x^12-25302258130671850*x^11+1000337781397362862*x^10-7709908154554062116*x^9+28992452711792779249*x^8-52916196926862805102*x^7+29422192713421978873*x^6+13659727327936090732*x^5+15671251179023563628*x^4+15074383085029309092*x^3+24075625195314205057*x^2-10144837465570789446*x+64250789525986279259


text:
@testset "operators" begin

@testset "squash" begin

    @test squash(ValGraph(126)) isa ValGraph{Int8}
    @test squash(ValDiGraph(127)) isa ValDiGraph{UInt8}
    @test squash(ValOutDiGraph{UInt8}(255)) isa ValOutDiGraph{Int16} # we try to avoid nv(g) == typemax(eltype(g))

end

@testset "reverse" begin

    g1 = DummyValGraph(ValGraph(2, edgeval_types=(a=String,)))
    add_edge!(g1, 1, 2, val="xyz")
    @test reverse(g1) === g1

    @test is_directed(typeof(SimpleValueGraphs.Reverse(g1))) == false

    g2 = DummyValGraph(ValDiGraph(5, vertexval_types=(a=Int,), vertexval_init=v -> (v,),
                                  edgeval_types=(b=String,), graphvals=(c=100,)))
    add_edge!(g2, 1, 2, val="abc")
    add_edge!(g2, 2, 1, val="def")
    add_edge!(g2, 3, 4, val="ghi")
    add_edge!(g2, 4, 4, val="jkl")

    rg = reverse(g2)

    @test reverse(rg) == g2

    @test is_directed(typeof(g2)) == true

    @test rg isa SimpleValueGraphs.Reverse
    @test rg.graph === g2

    @test nv(rg) == nv(g2) == 5
    add_vertex!(rg, val=6)
    @test nv(rg) == nv(g2) == 6

    @test get_vertexval(rg, 6) == get_vertexval(g2, 6) == 6
    set_vertexval!(rg, 6, :, (600, ))
    @test get_vertexval(rg, 6, :) == get_vertexval(g2, 6, :) == (a=600,)
    set_vertexval!(rg, 6, 1, 600)
    @test get_vertexval(rg, 6, 1) == get_vertexval(g2, 6, 1) == 600
    set_vertexval!(rg, 6, :a, 600)
    @test get_vertexval(rg, 6, :a) == get_vertexval(g2, 6, :a) == 600

    @test get_edgeval(rg, 1, 2) == get_edgeval(g2, 2, 1) == "def"
    set_edgeval!(rg, 1, 2, :, ("DEF",))
    @test get_edgeval(rg, 1, 2, :) == get_edgeval(g2, 2, 1, :) == (b="DEF",)
    set_edgeval!(rg, 2, 1, 1, "ABC")
    @test get_edgeval(rg, 2, 1, 1) == get_edgeval(g2, 1, 2, 1) == "ABC"
    set_edgeval!(rg, 4, 4, :b, "JKL")
    @test get_edgeval(rg, 4, 4, :b) == get_edgeval(g2, 4, 4, :b) == "JKL"

    @test ne(rg) == ne(g2) == 4
    add_edge!(rg, 1, 6, val="mno")
    @test get_edgeval(rg, 1, 6) == get_edgeval(g2, 6, 1) == "mno"
    @test has_edge(rg, 1, 6)
    @test has_edge(g2, 6, 1)
    @test !has_edge(rg, 6, 1)
    @test !has_edge(g2, 1, 6)
    @test ne(rg) == ne(g2) == 5

    @test rem_edge!(rg, 1, 6) == true
    @test !has_edge(rg, 1, 6)
    @test !has_edge(g2, 6, 1)
    @test ne(rg) == ne(g2) == 4
    @test rem_edge!(rg, 1, 6) == false
    @test ne(rg) == ne(g2) == 4

    @test get_graphval(rg) == get_graphval(g2) == 100
    set_graphval!(rg, 1, 200)
    @test get_graphval(rg, 1) == get_graphval(g2, 1) == 200
    set_graphval!(rg, :c, 300)
    @test get_graphval(rg, :c) == get_graphval(g2, :c) == 300
    set_graphval!(rg, :, (300,))
    @test get_graphval(rg, :) == get_graphval(g2, :) == (c=300,)

    @test outneighbors(rg, 1) == inneighbors(g2, 1)
    @test outneighbors(rg, 2) == inneighbors(g2, 2)
    @test outneighbors(rg, 3) == inneighbors(g2, 3)
    @test outneighbors(rg, 4) == inneighbors(g2, 4)

    @test inneighbors(rg, 1) == outneighbors(g2, 1)
    @test inneighbors(rg, 2) == outneighbors(g2, 2)
    @test inneighbors(rg, 3) == outneighbors(g2, 3)
    @test inneighbors(rg, 4) == outneighbors(g2, 4)

    @test outedgevals(rg, 1, 1) == inedgevals(g2, 1, 1)
    @test outedgevals(rg, 2, :b) == inedgevals(g2, 2, :b)
    @test outedgevals(rg, 3, :) == inedgevals(g2, 3, :)

    @test inedgevals(rg, 2, 1) == outedgevals(g2, 2, 1)
    @test inedgevals(rg, 3, :b) == outedgevals(g2, 3, :b)
    @test inedgevals(rg, 4, :) == outedgevals(g2, 4, :)

    # TODO Tests currently disabled as we do not have zero implemented on ValDiGraph
    # @test zero(typeof(rg)) isa typeof(rg)
    # @test nv(zero(typeof(rg))) == 0
end

end


text:
module SMLMSim

using SMLMData
using Distributions

include("typedefs.jl")
include("molecules.jl")
include("cameras.jl")
include("patterns.jl")
include("sim.jl")
include("interface.jl")

end


text:
__precompile__(true)

"""
Main module for `CarrierCapture.jl` -- A set of codes to compute carrier capture and recombination rates in semiconducting compounds.

Two structs are exported from this module for public use:

- [`potential`](@ref): Potential.
- [`conf_coord`](@ref): Configuration coordinate.
- [`Plotter`](@ref): Auxiliary submodule for plotting.
- [`Brooglie`](@ref): 1D Shrödinger equation solver.

"""
module CarrierCapture
    using DocStringExtensions
    # using Reexport
    using LsqFit, Polynomials, Dierckx, Interpolations, Roots

    # using Plots, LaTeXStrings
    using DataFrames
    # @reexport using Plots, LaTeXStrings
    # @reexport using DataFrames

    include("Brooglie.jl")

    include("Potential.jl")
    include("CaptureRate.jl")
    include("Plotter.jl")
    include("paramScan.jl")
end # module


text:
#=
Idmrg environments are only to be used internally.
They have to be updated manually, without any kind of checks
=#

struct IDMRGEnv{H<:MPOHamiltonian,V<:MPSTensor}
    opp :: H
    lw :: PeriodicArray{V,2}
    rw :: PeriodicArray{V,2}
end

function IDMRGEnv(state::InfiniteMPS,env::MPOHamInfEnv)
    state === env.dependency || recalculate!(env,state);
    IDMRGEnv(env.opp,deepcopy(env.lw),deepcopy(env.rw));
end

leftenv(envs::IDMRGEnv,pos::Int,state=nothing) = envs.lw[pos,:];
function setleftenv!(envs::IDMRGEnv,pos,lw)
    envs.lw[pos,:] = lw[:]
end

rightenv(envs::IDMRGEnv,pos::Int,state=nothing) = envs.rw[pos,:];
function setrightenv!(envs::IDMRGEnv,pos,rw)
    envs.rw[pos,:] = rw[:]
end


text:
export proxLLR!, normLLR

"""
    proxLLR!(x::Vector{T}, λ::Float64=1e-6; kargs...) where T

proximal map for LLR regularization using singular-value-thresholding

# Arguments
* `x::Vector{T}`                - Vector to apply proximal map to
* `λ::Float64`                  - regularization parameter
* `shape::Tuple{Int}=[]`        - dimensions of the image
* `blockSize::Tuple{Int}=[2;2]` - size of patches to perform singluar value thresholding on
* `randshift::Bool=true`        - randomly shifts the patches to ensure translation invariance
"""
function proxLLR!(x::Vector{T}, λ::Float64=1e-6; shape::NTuple=[],
   blockSize::Array{Int64,1}=[2; 2], randshift::Bool=true, kargs...) where T
  # xᵖʳᵒˣ = zeros(T,size(x))
  N = prod(shape)
  # K = floor(Int,length(x)/N)
  x[:] = vec( svt(x[:], shape, λ; blockSize=blockSize, randshift=randshift, kargs...) )

  return x
end

function svt(x::Vector{T}, shape::Tuple, λ::Float64=1e-6;
   blockSize::Array{Int64,1}=[2; 2], randshift::Bool=true, kargs...) where T

  x = reshape( x, tuple( shape...,floor(Int64, length(x)/prod(shape)) ) )

  Wy = blockSize[1]
  Wz = blockSize[2]

  if randshift
    # Random.seed!(1234)
    shift_idx = [rand(1:Wy) rand(1:Wz) 0]
    x = circshift(x, shift_idx)
  end

  ny, nz, K = size(x)

  # reshape into patches
  L = floor(Int,ny*nz/Wy/Wz) # number of patches, assumes that image dimensions are divisble by the blocksizes

  xᴸᴸᴿ = zeros(T,Wy*Wz,L,K)
  for i=1:K
    xᴸᴸᴿ[:,:,i] = im2colDistinct(x[:,:,i], (Wy,Wz))
  end
  xᴸᴸᴿ = permutedims(xᴸᴸᴿ,[1 3 2])

  # threshold singular values
  for i = 1:L
    if xᴸᴸᴿ[:,:,i] == zeros(T, Wy*Wz,K)
      continue
    end
    SVDec = svd(xᴸᴸᴿ[:,:,i])
    proxL1!(SVDec.S,λ)
    xᴸᴸᴿ[:,:,i] = SVDec.U*Matrix(Diagonal(SVDec.S))*SVDec.Vt
  end

  # reshape into image
  xᵗʰʳᵉˢʰ = zeros(T,size(x))
  for i = 1:K
    xᵗʰʳᵉˢʰ[:,:,i] = col2imDistinct( xᴸᴸᴿ[:,i,:], (Wy,Wz), (ny,nz) )
  end

  if randshift
    xᵗʰʳᵉˢʰ = circshift(xᵗʰʳᵉˢʰ, -1*shift_idx)
  end

  if !isempty(shape)
    xᵗʰʳᵉˢʰ = reshape( xᵗʰʳᵉˢʰ, prod(shape),floor( Int, length(xᵗʰʳᵉˢʰ)/prod(shape) ) )
  end

  return xᵗʰʳᵉˢʰ
end

"""
    normLLR(x::Vector{T}, λ::Float64; kargs...) where T

returns the value of the LLR-regularization term.
Arguments are the same is in `proxLLR!`
"""
function normLLR(x::Vector{T}, λ::Float64; shape::NTuple=[], L=1, blockSize::Array{Int64,1}=[2; 2], randshift::Bool=true, kargs...) where T

  N = prod(shape)
  K = floor(Int,length(x)/(N*L))
  normᴸᴸᴿ = 0.
  for i = 1:L
    normᴸᴸᴿ +=  blockNuclearNorm(x[(i-1)*N*K+1:i*N*K], shape; blockSize=blockSize, randshift=randshift, kargs...)
  end

  return λ*normᴸᴸᴿ
end

function blockNuclearNorm(x::Vector{T}, shape::Tuple; blockSize::Array{Int64,1}=[2; 2],
      randshift::Bool=true, kargs...) where T
    x = reshape( x, tuple( shape...,floor(Int64, length(x)/prod(shape)) ) )

    Wy = blockSize[1]
    Wz = blockSize[2]

    if randshift
      srand(1234)
      shift_idx = [rand(1:Wy) rand(1:Wz) 0]
      x = circshift(x, shift_idx)
    end

    ny, nz, K = size(x)

    # reshape into patches
    L = floor(Int,ny*nz/Wy/Wz) # number of patches, assumes that image dimensions are divisble by the blocksizes

    xᴸᴸᴿ = zeros(T,Wy*Wz,L,K)
    for i=1:K
      xᴸᴸᴿ[:,:,i] = im2colDistinct(x[:,:,i], (Wy,Wz))
    end
    xᴸᴸᴿ = permutedims(xᴸᴸᴿ,[1 3 2])

    # L1-norm of singular values
    normᴸᴸᴿ = 0.
    for i = 1:L
      SVDec = svd(xᴸᴸᴿ[:,:,i])
      normᴸᴸᴿ += norm(SVDec.S,1)
    end

    return normᴸᴸᴿ
end


text:
mutable struct SuperObjCollision
    superObj::Array{Object3D,1}
    function SuperObjCollision()
        new(Array{Object3D,1}())
    end
end

mutable struct SuperObjMass
    superObj::Array{Object3D,1}
    function SuperObjMass()
        new(Array{Object3D,1}())
    end
end

mutable struct SuperObjForce
    superObj::Array{Object3D,1}
    function SuperObjForce()
        new(Array{Object3D,1}())
    end
end

mutable struct SuperObjVisu
    superObj::Array{Object3D,1}
    function SuperObjVisu()
        new(Array{Object3D,1}())
    end
end

mutable struct SuperObjsRow
    superObjCollision::SuperObjCollision
    superObjMass::SuperObjMass
    superObjForce::SuperObjForce
    superObjVisu::SuperObjVisu
    noCPair::Array{Int64,1}
    function SuperObjsRow()
        new(SuperObjCollision(), SuperObjMass(), SuperObjForce(), SuperObjVisu(), Array{Int64,1}())
    end
end


text:
function FieldAltitudeDetector(arg0::RealFieldElement, arg1::BodyShape)
    return FieldAltitudeDetector((RealFieldElement, BodyShape), arg0, arg1)
end

function FieldAltitudeDetector(arg0::RealFieldElement, arg1::RealFieldElement, arg2::BodyShape)
    return FieldAltitudeDetector((RealFieldElement, RealFieldElement, BodyShape), arg0, arg1, arg2)
end

function FieldAltitudeDetector(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::BodyShape)
    return FieldAltitudeDetector((RealFieldElement, RealFieldElement, RealFieldElement, BodyShape), arg0, arg1, arg2, arg3)
end

function g(obj::FieldAltitudeDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_altitude(obj::FieldAltitudeDetector)
    return jcall(obj, "getAltitude", RealFieldElement, ())
end

function get_body_shape(obj::FieldAltitudeDetector)
    return jcall(obj, "getBodyShape", BodyShape, ())
end


text:
#=
Set
=#

"""
Space complexity: O(n)
"""
struct SimpleSet{T} <: AbstractSet{T}
    d::SimpleDict{T,Nothing}

    function SimpleSet{T}(; capacity::Integer = 8) where {T}
        new{T}(SimpleDict{T,Nothing}(capacity = capacity))
    end
end

#=
Constructors
=#

function SimpleSet{T}(xs::AbstractArray{T}) where T
    ss = SimpleSet{T}()
    for x in xs
        push!(ss, x)
    end
    return ss
end

SimpleSet(xs::AbstractArray{T}) where T = SimpleSet{T}(xs)

#=
Core methods
=#

"""
Time complexity: O(1)
"""
function Base.push!(ss::SimpleSet{T}, item::T) where T
    ss.d[item] = nothing
    return ss
end

"""
Time complexity: O(1)
"""
Base.in(item::T, ss::SimpleSet{T}) where T = haskey(ss.d, item)

"""
Time complexity: O(1)
"""
function Base.delete!(ss::SimpleSet{T}, item::T) where T
    delete!(ss.d, item)
    return ss
end

#=
Interface
=#

Base.length(ss::SimpleSet) = ss.d.count

Base.iterate(ss::SimpleSet, i::Integer = 1) = iterate(ss.d, i)


text:
include(joinpath(pwd(), "src/controller/newton_structure_solver/methods.jl"))

sim = get_simulation("flamingo", "flat_2D_lc", "flat")
model = sim.model
env = sim.env

ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, flat_2D_lc,
    joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))

copy_traj = deepcopy(ref_traj)

ip_opts = eval(interior_point_options(:interior_point))(
			κ_init = 1.0e-4,
			κ_tol = 2.0 * 1.0e-4,
			r_tol = 1.0e-8,
			diff_sol = true,
			solver = :empty_solver)

im_traj = ImplicitTraj(ref_traj, sim,
	ip_type = :interior_point,
	κ = 1.0e-4,
	mode = :configuration,
	opts=ip_opts)

nq = model.dim.q
n = 2 * nq
m = model.dim.u
T = 15

nz = n * (T - 1) + m * (T - 1)
nd = n * (T - 1)

u_idx = [collect((t - 1) * (m + n) .+ (1:m)) for t = 1:T-1]
x_idx = [collect((t - 1) * (m + n) + m .+ (1:n)) for t = 1:T-1]
n_idx = [(t - 1) * n .+ (1:n) for t = 1:T-1]

# implicit_dynamics!(im_traj, sim, copy_traj, κ = 1.0e-4)
# im_traj.dq2
# [norm(im_traj.dq2[t], Inf) for t = 1:T-1]
# im_traj.ip[1].z

implicit_dynamics!(im_traj, model, env, ref_traj.u, ref_traj.q[1:T-1], ref_traj.q[2:T], T, copy_traj, κ = 1.0e-4)
[norm(im_traj.dq2[t] - ref_traj.q[t+2], Inf) for t = 1:T-1]

# A = [[zeros(nq, nq) I; im_traj.δq0[t] im_traj.δq1[t]] for t = 1:T-1]
# B = [[zeros(nq, m); im_traj.δu1[t]] for t = 1:T-1]
Q = [Diagonal(0.1 * ones(n)) for t = 1:T]
V = [Diagonal(0.001 * ones(nq)) for t = 1:T]
R = [Diagonal(0.01 * ones(m)) for t = 1:T-1]

x_ref = [[ref_traj.q[t]; ref_traj.q[t+1]] for t = 1:T]
u_ref = [ref_traj.u[t] for t = 1:T-1]

x_init = [ref_traj.q[1]; ref_traj.q[2]]
u = [ref_traj.u[t] for t = 1:T-1]
x = deepcopy(x_ref[2:T])
z = [vcat([[u[t]; x[t]] for t = 1:T-1]...); zeros(nd)]

function constraints(z)
	c = zeros(eltype(z), nd)
	x_traj = [z[x_idx[t]] for t = 1:T-1]
	u_traj = [z[u_idx[t]] for t = 1:T-1]
	qa = [x_init[1:nq], [x_traj[t][1:nq] for t = 1:T-1]...]
	qb = [x_init[nq .+ (1:nq)], [x_traj[t][nq .+ (1:nq)] for t = 1:T-1]...]

	# linear_contact_implicit_dynamics!(lci_traj, u_traj, qa, qb, T-1)
	implicit_dynamics!(im_traj, model, env, u_traj, qa, qb, T, copy_traj, κ = 1.0e-4)

	for t = 1:T-1
		x2 = z[x_idx[t]]
		x1 = (t == 1 ? x_init : z[x_idx[t-1]])
		u1 = z[u_idx[t]]
		c[n_idx[t]] = x2 - [x1[nq .+ (1:nq)]; im_traj.dq2[t]]#A[t] * x1 - B[t] * u1
		# c[n_idx[t]] = x2 - A[t] * x1 - B[t] * u1
		# c[n_idx[t]] = x2 - [zeros(nq, nq) I; im_traj.δq0[t] im_traj.δq1[t]] * x1 - [zeros(nq, m); im_traj.δu1[t]] * u1
	end

	return c
end

constraints(z)
norm(constraints(z), Inf)

function constraints_jacobian(z)
	C = zeros(eltype(z), nd, nz)

	x_traj = [z[x_idx[t]] for t = 1:T-1]
	u_traj = [z[u_idx[t]] for t = 1:T-1]
	qa = [x_init[1:nq], [x_traj[t][1:nq] for t = 1:T-1]...]
	qb = [x_init[nq .+ (1:nq)], [x_traj[t][nq .+ (1:nq)] for t = 1:T-1]...]

	# linear_contact_implicit_dynamics!(lci_traj, u_traj, qa, qb, T-1)
	implicit_dynamics!(im_traj, model, env, u_traj, qa, qb, T, copy_traj, κ = 1.0e-4)

	for t = 1:T-1
		C[n_idx[t], u_idx[t]] = -[zeros(nq, m); im_traj.δu1[t]]
		C[n_idx[t], x_idx[t]] = Diagonal(ones(n))
		t == 1 && continue
		C[n_idx[t], x_idx[t-1]] = -[zeros(nq, nq) I; im_traj.δq0[t] im_traj.δq1[t]]
	end

	return C
end

constraints_jacobian(z)

function objective_gradient(z)
	j = zeros(eltype(z), nz)
	for t = 1:T-1
		u1 = z[u_idx[t]]
		x1 = z[x_idx[t]]
		j[u_idx[t]] = R[t] * (u1 - u_ref[t])
		j[x_idx[t]] = Q[t+1] * (x1 - x_ref[t+1])
		j[x_idx[t][1:nq]] -= V[t+1] * (x1[nq .+ (1:nq)] - x1[1:nq])
		j[x_idx[t][nq .+ (1:nq)]] += V[t+1] * (x1[nq .+ (1:nq)] - x1[1:nq])
	end
	return j
end

objective_gradient(z)
norm(objective_gradient(z))

function objective_hessian(z)
	J = zeros(eltype(z), nz, nz)
	for t = 1:T-1
		J[u_idx[t], u_idx[t]] = R[t]
		J[x_idx[t], x_idx[t]] = Q[t+1]
		J[x_idx[t][1:nq], x_idx[t][nq .+ (1:nq)]] -= V[t+1]
		J[x_idx[t][nq .+ (1:nq)], x_idx[t][1:nq]] -= V[t+1]
	end
	return J
end

objective_hessian(z)

ρ = 0.0
function hessian(z)
	J = objective_hessian(z)
	C = constraints_jacobian(z)

	return [J + ρ * I C'; C -ρ * I]
end

function get_lagrangian(z)
	y = z[nz .+ (1:nd)]
	return objective_gradient(z) + constraints_jacobian(z)' * y
end

function gradient(z)
	y = z[nz .+ (1:nd)]
	return [objective_gradient(z) + constraints_jacobian(z)' * y; constraints(z)]
end

# norm(get_lagrangian(z), 1)

res_norm = norm(gradient(z), 1)

Δ = hessian(z) \ gradient(z)
α = 1.0
z_cand = z - α * Δ

# norm(get_lagrangian(z_cand), 1)
# norm(constraints(z_cand), 1)
res_cand_norm = norm(gradient(z_cand), 1)

z = z_cand


text:
function expr_exec(context, grn_obj)
  ccall((:grn_expr_exec, libgroonga),
        Void,
        (Ptr{Void}, Ptr{Void}, Int),
        context,
        grn_obj,
        0)
end


text:
# NOTE: You pass in weights using a row vector, e.g. [500 300 100], then I convert it to a 1-dim column vector

function make_dense_nn_regression(count_weights::Array{Int64,1},size_output::Int;regu_coef = 0.0)
    count_weights = vec(count_weights)
    count_hidden_layers = size(count_weights, 1)
    count_inner_product_layers = count_hidden_layers + 1

    # Make model name string
    model_name = "dense"
    for i = 1:count_hidden_layers
        model_name = model_name * "_" * string(count_weights[i])
    end

    # Returns a function that returns a function to create model
    nn_factory = function()
        common_layers = Array(InnerProductLayer, count_inner_product_layers)

        # Create the input, hidden, and output layers
        for i = 1:count_inner_product_layers
            # Specify the from and to Symbols
            if i == 1
                from_Symbol = :data
            else
                from_Symbol = Symbol( "ip$(i-1)" )
            end
            to_Symbol = Symbol( "ip$(i)" )
            # Create the layer
            if i != count_inner_product_layers
                common_layers[i] = InnerProductLayer(
                                        name = string(to_Symbol),
                                        output_dim = count_weights[i],
                                        neuron = Neurons.ReLU(),
                                        bottoms = [from_Symbol],
                                        tops = [to_Symbol],
                weight_regu = L2Regu(regu_coef))
            else
                common_layers[i] = InnerProductLayer(
                                            name = string(to_Symbol),
                                            output_dim = size_output,
                                            bottoms = [from_Symbol],
                                            tops = [to_Symbol],
                                            weight_regu = L2Regu(regu_coef))
            end
        end

        last_dense = Symbol("ip$(count_inner_product_layers)")

        storage_layer = MemoryOutputLayer( name = "storage", bottoms = [last_dense])

        loss_layer = SquareLossLayer(name = "loss",
                                    bottoms = [Symbol( "ip$(count_inner_product_layers)" ),:ground_truth],
				                    #normalize = :no ) # <= I don't think normalize is necessary...
                                    )

        return (common_layers, storage_layer, loss_layer)
    end
    return (nn_factory, model_name)
end


text:
dtexec = [0.0];
constraint_time = nothing;


text:
## This file is outdated

export timeevolution, ETDRK4, BDF2, BDF4, BDF22

function RK(L,y,h)
    k1=L(y)
    k2=L(y+.5h*k1)
    k3=L(y+.5h*k2)
    k4=L(y+h*k3)

    y+h*(k1+2k2+2k3+k4)/6.
end

function BDF2(B,A::Operator,g::Function,bcs,u0,h,m,glp,tol=1000eps())
    SBDF2 = [B;I-2.0/3.0*h*A]

    u1=u0
    u2=chop(RK(g,u1,h),tol)
    u2,u1  = chop(SBDF2\[bcs;1/3.0*(4u2-u1)],tol),u2
    push!(glp,u2)
    yield()

    for k=1:m
        u2,u1 = chop(RK(g,u2,h),tol),u2
        u2,u1  = chop(SBDF2\[bcs;1/3.0*(4u2-u1)],tol),u2
        push!(glp,u2)
        yield()
    end

    u2
end

function BDF2(L::Operator,N::Function,u0,h,m,glp,tol=1000eps())
    SBDF2 = I-2.0/3.0*h*L

    u1=u0
    u2=chop(RK(N,u1,h),tol)
    u2,u1  = chop(SBDF2\(1/3.0*(4u2-u1)),tol),u2
    push!(glp,u2)
    yield()

    for k=1:m
        u2,u1 = chop(RK(N,u2,h),tol),u2
        u2,u1  = chop(SBDF2\(1/3.0*(4u2-u1)),tol),u2
        push!(glp,u2)
        yield()
    end

    u2
end

function ETDRK4(L::Operator,N::Function,u,t,h,m,glp,tol=10eps())
    z = L*h
    ez,ez2,h2ez2m1 = cache(exp(z)),cache(exp(z/2)),cache(h/2*expm1(z/2))
    hezα,h2ezβ,hezγ = cache(h*expα(z)),cache(2h*expβ(z)),cache(h*expγ(z))
    push!(glp,u)
    yield()

    for k=1:m
        ez2u = ez2*u
        Nut = chop(N(u,t),tol)

        a = ez2u + h2ez2m1*Nut
        Nath2 = chop(N(a,t+h/2),tol)
        b = ez2u + h2ez2m1*Nath2
        Nbth2 = chop(N(b,t+h/2),tol)
        c = ez2*a + h2ez2m1*(2Nbth2-Nut)

        u = chop(ez*u + hezα*Nut + h2ezβ*(Nath2+Nbth2) + hezγ*N(c,t+h),tol)
        push!(glp,u)
        yield()
        t+=h
    end

    u,t
end

## Multivariate

function BDF4(B::Vector,op::Operator,bcs::Vector,uin::MultivariateFun,h::Real,m::Integer,glp)
    nt=size(uin,2)
    d=domain(uin)
    SBE   = discretize([B;I-h*op],d,nt)            # backward euler for first 2 time steps
    SBDF2 = discretize([B;I-2.0/3.0*h*op],d,nt)    # BDF formula for subsequent itme steps
    SBDF3 = discretize([B;I-6.0/11.0*h*op],d,nt)    # BDF formula for subsequent itme steps
    SBDF4 = discretize([B;I-12.0/25.0*h*op],d,nt)    # BDF formula for subsequent itme steps

    u1=uin
    u2=SBE\[bcs;u1]
    push!(glp,pad(u2,80,80))
    yield()
    u3=SBDF2\[bcs;1/3.0*(4u2-u1)]
    push!(glp,pad(u3,80,80))#updates window
    yield()
    u4=SBDF3\[bcs;1/11.0*(18u3-9u2+2u1)]
    push!(glp,pad(u4,80,80))#updates window
    yield()

    for k=1:m
        u4,u3,u2,u1  = SBDF4\[bcs;1/25.0*(48u4-36u3+16u2-3u1)],u4,u3,u2

        push!(glp,pad(u4,80,80))#updates window
        yield()
    end
    u4
end

BDF4(B::Vector,op::Operator,uin::MultivariateFun,h::Real,m::Integer,glp)=BDF4(B,op,zeros(length(B)),uin,h,m,glp)

function BDF22(B::Vector,op::Operator,bcs::Vector,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,m::Integer,glp)
    nt=size(uin[1],2)
    SBE  = discretize([B;I-h^2*op],domain(uin[1]),nt)            # backward euler for first 2 time steps
    SBDF = discretize([B;I-4.0/9.0*h^2*op],domain(uin[1]),nt)    # BDF formula for subsequent itme steps

    u1,u2=uin
    u3 =SBE\[bcs;2u2-u1]
    push!(glp,pad(u3,80,80))
    yield()
    u4 =SBE\[bcs;2u3-u2]
    push!(glp,pad(u4,80,80))
    yield()

    for k=1:m
        u4,u3,u2,u1  = SBDF\[bcs;1/9.0*(24u4-22u3+8u2-u1)],u4,u3,u2

        push!(glp,pad(u4,80,80)) #updates window
        yield()
    end
    u4
end

function BDF22(B::Vector,op::Operator,g::Function,bcs::Vector,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,m::Integer,glp)
    nt=size(uin[1],2)
    SBE  = discretize([B;I-h^2*op],domain(uin[1]),nt)            # backward euler for first 2 time steps
    SBDF = discretize([B;I-4.0/9.0*h^2*op],domain(uin[1]),nt)    # BDF formula for subsequent itme steps

    u1,u2=uin
    u3 =SBE\[bcs;2u2-u1]
    push!(glp,pad(u3,80,80))
    yield()
    u4,u3,u2,u1=u3,u2,u1,u1
    for k=1:m
        u4,u3,u2,u1 =2u4 - u3 +h^2*g(u4),u4,u3,u2
        u4,u3,u2,u1  = SBDF\[bcs;1/9.0*(24u4-22u3+8u2-u1)],u4,u3,u2
        push!(glp,pad(u4,80,80))
        yield()
    end
    u4
end

BDF22(B::Vector,op::Operator,uin::MultivariateFun,h::Real,m::Integer,glp)=BDF22(B,op,zeros(length(B)),(uin,uin),h,m,glp)
BDF22(B::Vector,op::Operator,g::Function,uin::MultivariateFun,h::Real,m::Integer,glp)=BDF22(B,op,g,zeros(length(B)),(uin,uin),h,m,glp)

## GLPlot routines

#u_t = op*u
function timeevolution(B::Vector,op,bcs::Vector,uin::MultivariateFun,h::Real,m::Integer,glp)
    require("GLPlot")
    setplotter("GLPlot")
    nt=size(uin,2)
    d=domain(uin)
    SBE   = discretize([B;I-h*op],d,nt)            # backward euler for first 2 time steps
    SBDF2 = discretize([B;I-2.0/3.0*h*op],d,nt)    # BDF formula for subsequent itme steps
    SBDF3 = discretize([B;I-6.0/11.0*h*op],d,nt)    # BDF formula for subsequent itme steps
    SBDF4 = discretize([B;I-12.0/25.0*h*op],d,nt)    # BDF formula for subsequent itme steps

    u1=chop(uin,1000eps())
    u2=chop(SBE\[bcs;u1],1000eps())
    plot(pad(u2,80,80),glp...)#updates window
    u3=chop(SBDF2\[bcs;1/3.0*(4u2-u1)],1000eps())
    plot(pad(u3,80,80),glp...)#updates window
    u4=chop(SBDF3\[bcs;1/11.0*(18u3-9u2+2u1)],1000eps())
    plot(pad(u4,80,80),glp...)#updates window

    for k=1:m
        u4,u3,u2,u1  = chop(SBDF4\[bcs;1/25.0*(48u4-36u3+16u2-3u1)],1000eps()),u4,u3,u2

        plot(pad(u4,80,80),glp...)#updates window
    end
    u4
end

function timeevolution(B::Vector,op,bcs::Vector,uin::MultivariateFun,h::Real,m=5000)
    require("GLPlot")
    setplotter("GLPlot")
    timeevolution(B,op,bcs,uin,h,m,plot(pad(uin,80,80)))
end

timeevolution(B::Vector,op,uin::MultivariateFun,h::Real,dat...)=timeevolution(B,op,zeros(length(B)),uin,h,dat...)

timeevolution(B::Vector,op,bcs::Vector,uin::Fun,dat...)=timeevolution(B,op,bcs,ProductFun(uin),dat...)
timeevolution(B::Vector,op,uin::Fun,dat...)=timeevolution(B,op,ProductFun(uin),dat...)

timeevolution(B::Operator,dat...)=timeevolution([B],dat...)

#u_tt = op*u
function timeevolution2(B::Vector,op,bcs::Vector,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,m::Integer,glp)
    require("GLPlot")
    setplotter("GLPlot")
    nt=size(uin[1],2)
    SBE  = discretize([B;I-h^2*op],domain(uin[1]),nt)            # backward euler for first 2 time steps
    SBDF = discretize([B;I-4.0/9.0*h^2*op],domain(uin[1]),nt)    # BDF formula for subsequent itme steps

    u1,u2=chop(uin[1],1000eps()),chop(uin[2],1000eps())
    u3 =chop(SBE\[bcs;2u2-u1],1000eps())
    u4 =chop(SBE\[bcs;2u3-u2],1000eps())

    for k=1:m
        u4,u3,u2,u1  = chop(SBDF\[bcs;1/9.0*(24u4-22u3+8u2-u1)],1000eps()),u4,u3,u2

        plot(pad(u4,80,80),glp...)#updates window
    end
    u4
end

#u_tt = op*u
function timeevolution2(B::Vector,op,g::Function,bcs::Vector,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,m::Integer,glp)
    require("GLPlot")
    setplotter("GLPlot")
    nt=size(uin[1],2)
    SBE  = discretize([B;I-h^2*op],domain(uin[1]),nt)            # backward euler for first 2 time steps
    SBDF = discretize([B;I-4.0/9.0*h^2*op],domain(uin[1]),nt)    # BDF formula for subsequent itme steps

    u1,u2=chop(uin[1],1000eps()),chop(uin[2],1000eps())
    u3 =chop(SBE\[bcs;2u2-u1],1000eps())
    u4 =chop(2u3 - u2 +h^2*g(u3),1000eps())
    u4,u3,u2,u1  = chop(SBDF\[bcs;1/9.0*(24u4-22u3+8u2-u1)],1000eps()),u4,u3,u2
    plot(pad(u4,80,80),glp...)#updates window

    for k=1:m
        u4,u3,u2,u1  = chop(2u4 - u3 +h^2*g(u4),1000eps()),u4,u3,u2
        u4,u3,u2,u1  = chop(SBDF\[bcs;1/9.0*(24u4-22u3+8u2-u1)],1000eps()),u4,u3,u2
        plot(pad(u4,80,80),glp...)#updates window
    end
    u4
end

function timeevolution2(B::Vector,op,uin::Tuple{MultivariateFun,MultivariateFun},bcs::Vector,h::Real,m=5000)
    require("GLPlot")
    setplotter("GLPlot")
    timeevolution2(B,op,bcs,uin,h,m,plot(pad(uin[end],80,80)))
end

function timeevolution2(B::Vector,op,g::Function,uin::Tuple{MultivariateFun,MultivariateFun},bcs::Vector,h::Real,m=5000)
    require("GLPlot")
    setplotter("GLPlot")
    timeevolution2(B,op,g,bcs,uin,h,m,plot(pad(uin[end],80,80)))
end

timeevolution2(B::Vector,op,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,dat...)=timeevolution2(B,op,uin,zeros(length(B)),h,dat...)
timeevolution2(B::Vector,op,uin::MultivariateFun,dat...)=timeevolution2(B,op,(uin,uin),dat...)
timeevolution2(B::Operator,dat...)=timeevolution2([B],dat...)

timeevolution2(B::Vector,op,g::Function,uin::Tuple{MultivariateFun,MultivariateFun},h::Real,dat...)=timeevolution2(B,op,g,uin,zeros(length(B)),h,dat...)
timeevolution2(B::Vector,op,g::Function,uin::MultivariateFun,dat...)=timeevolution2(B,op,g,(uin,uin),dat...)

timeevolution2(B::Vector,op,g::Function,uin::Tuple{Fun,Fun},dat...)=timeevolution2(B,op,g,(ProductFun(uin[1]),ProductFun(uin[2])),dat...)
timeevolution2(B::Vector,op,g::Function,uin::Fun,dat...)=timeevolution2(B,op,g,ProductFun(uin),dat...)
timeevolution2(B::Vector,op,uin::Fun,dat...)=timeevolution2(B,op,ProductFun(uin),dat...)

timeevolution(o::Integer,dat...)=o==2?timeevolution2(dat...):timeevolution(dat...)


text:
export dataLayer

type dataLayer<:Layer
  bottom ::Array{Layer}
  top    ::Array{Layer}
  topRange ::Vector

  rₒ     ::Array{Float64}  #data
  rᵢ     ::Array{Float64} #objective

  lock   ::Bool
end
dataLayer(inDim::Int,outDim::Int)=
  dataLayer(Layer[],Layer[],Vector[],rand(outDim,batchsize),rand(inDim,batchsize),false)