text:
"""
    struct BuffersFaceIntegral

Stores all temporary buffers needed during face integrals.
Avoids costly re-allocation.
"""
struct BuffersFaceIntegral
    boundary_state::Array{Float64,1}
    boundary_flux::Array{Float64,2}
    numericalflux::Array{Float64,1}

    function BuffersFaceIntegral(ndofs)
        boundary_state = zeros(ndofs)
        boundary_flux = zeros(2,ndofs)
        numericalflux = zeros(ndofs)
        new(boundary_state, boundary_flux, numericalflux)
    end
end

"""
    local_lax_friedrichs(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)

Computes the local Lax-Friedrichs (or Rusanov) numerical flux for degrees of freedom
`dofs`, degrees of freedoms of the neighbor `dofsneigh`, flux `flux`, flux of neighbor
`fluxneigh`, cellsize `dx`.
All quantities are assumed to be represented by a basis on the reference line.
The face is parametrized by a index `normalidx`, where 1 stands for a face in x-direction
and 2 for a face in y-direction.
The sign of the outer normal of the face is given by `normalsign`.

The numerical flux is stored in `numericalflux`.
Method also returns the maximal eigenvalue.
"""
function local_lax_friedrichs(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)
    maxeigenval_left = max_eigenval(eq, dofs, normalidx)
    maxeigenval_right = max_eigenval(eq, dofsneigh, normalidx)
    maxeigenval = max(maxeigenval_left, maxeigenval_right)

    @views @inbounds numericalflux .= 0.5 .* (flux[normalidx,:] .+ fluxneigh[normalidx,:]) .+ 0.5 .* maxeigenval .* normalsign .* (dofs .- dofsneigh)
end

"""
    evaluate_face_integral(eq,globals, dofs, dofsneigh, flux, fluxneigh, buffers, cell, face, celldu)

Computes face integrals for equation `eq`, cell `cell` on face `face`.
Buffers are passed in `buffers`.
Result is stored in `celldu`.
"""
function evaluate_face_integral(eq,globals, dofs, dofsneigh, flux, fluxneigh, buffers, cell, face, celldu)
    normalidx = globals.normalidxs[face]
    normalsign = globals.normalsigns[face]

    # Compute Riemann solver (in normal)
    local_lax_friedrichs(eq, dofs, dofsneigh, flux, fluxneigh, cell.size[1], normalidx, normalsign, buffers.numericalflux)

    # Hack: Set numerical flux for materials to zero
    @views buffers.numericalflux[get_nvars(eq) + 1 : end] .= 0

    celldu .-= normalsign * area(cell) * buffers.numericalflux
end


text:
using CoinGeckoAPI
using Test

r = ping();

r = get_price("bitcoin", "usd");

@testset "CoinGeckoAPI.jl" begin
    # Write your tests here.

    @test 1 ==  1;

    @test ping()["gecko_says"] == "(V3) To the Moon!"

    @test get_price("bitcoin", "usd")["bitcoin"][:usd] |> typeof == Int64 ; # check if price is a number (whatever that is)

end

ext: jl


text:
using NumericalMethodsforEngineers, IterativeSolvers

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

f = cholesky(a)
y = f.L \ b
c = f.U \ y

# p_2_11 version

x0 = [1., 1., 1.]
tol = 1.0e-5
limit = 1000

function steepest_descent(a::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64};
  tol=1.0e-10, limit=1000, trace=false)

  x = copy(x0)
  trace && println("\nFirst few iterations:")
  r = b - a * x
  iters = 0
  while true
    iters += 1
    u = a * r
    α = r⋅r / dot(r, u)
    xnew = x + α * r
    r -= α * u
    if trace && iters <= 5
      println("After $iters iterations: x = $x")
    end
    (checkit(xnew, x, tol) || (iters == limit)) && break
    x = xnew
  end
  converged = iters <= limit
  iters, converged, x
end

iters, converged, x = steepest_descent(a, b, x0; tol=tol, limit=limit, trace=true)
println("\nSolution Vector: $x")
println("Converged: $(iters < limit)")
println("Number of iters: $(iters)")

@test round.(x; digits=3) == round.(c; digits=3)

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

tol = 1.0e-10
iters = 20
x = cg(a, b; maxiter=iters);

println("\nSolution Vector: $x")

@test round.(x; digits=6) == round.(c; digits=6)


text:
using Documenter, DocumenterMarkdown, DocumenterCitations, Hecke, Nemo, AbstractAlgebra, Pkg

@show Hecke.pkgdir

# Overwrite printing to make the header not full of redundant nonsense
# Turns
#   Hecke.Order - Method
# into
#   Order - Method

# To remove the '-'
# Documenter.Utilities.print_signature(io::IO, signature)        = print(io, signature)

# To remove the "Method", "Type", "Module" use the following
# Documenter.Utilities.doccat(b::Base.Docs.Binding, ::Type)  = ""
# doccat(::Type)     = ""
# doccat(::Module)   = ""

# Remove the module prefix
Base.print(io::IO, b::Base.Docs.Binding) = print(io, b.var)

DocMeta.setdocmeta!(Hecke, :DocTestSetup, :(using Hecke); recursive = true)

bib = CitationBibliography(joinpath(Hecke.pkgdir, "docs", "src", "Hecke.bib"), sorting = :nyt)

if Hecke.html_build[]
  makedocs(
      bib,
      format = Documenter.HTML(prettyurls = !local_build),
      doctest= true,
      modules = [Hecke],
      sitename = "Hecke documentation",
      pages = [
      "index.md",
      "Number fields" => [ "number_fields/intro.md",
                           "number_fields/fields.md",
                           "number_fields/elements.md",
                           "number_fields/internal.md"],
      "Function fields" => [ "function_fields/intro.md",
			     "function_fields/basics.md",
			     "function_fields/elements.md",
			     "function_fields/internal.md",
			     "function_fields/degree_localization.md"],
      "Orders" => [ "orders/introduction.md",
                    "orders/orders.md",
                    "orders/elements.md",
                    "orders/ideals.md",
                    "orders/frac_ideals.md"
                  ],
#      "Maximal Orders" => [ "MaximalOrders/Introduction.md",
#                            "MaximalOrders/Creation.md",
#                            "MaximalOrders/Elements.md",
#                            "MaximalOrders/Ideals.md"
#                          ],
      "abelian/introduction.md",
      "class_fields/intro.md",
      "sparse/intro.md",
      "FacElem.md"
      ]

  )
else
  makedocs(
      bib,
      doctest= true,
      modules = [Hecke],
      format = Markdown(),
  )

  docsdir = joinpath(@__DIR__, "build/")

  function _super_cool_example(f, overwrite = true)
    inside = false
    new_file = ""
    collapsing = false
    open(f) do file
      for ln in eachline(file);
        if startswith(ln, "<a id='Example")
          continue
        end
        if startswith(ln, "#### Example")
          if startswith(ln, "#### Example +")
            collapsing = true
          end
          ln = ""
          inside = true
        end
        if inside
          if startswith(ln, "```julia-repl")
            if collapsing
              line = "??? note \"Example\"\n    ```julia"
              collapsing = false
            else
              line = "!!! note \"Example\"\n    ```julia"
            end
          else
            line = "    " * ln
          end
        else
          line = ln;
        end

        if startswith(ln, "```") && !occursin("julia-repl", ln)
          inside = false
        end
        new_file = new_file * "\n" * line
      end
    end
    rm(f)
    open(f, "w") do file
      write(file, new_file)
    end
  end

  for (root, dirs, files) in walkdir(docsdir)
    for file in files
      filename = joinpath(root, file) # path to files
      #run(`sed -i 's/.*dash; \*Method.*/---/g' $filename`)
      #run(`sed -i 's/.*dash; \*Type.*/---/g' $filename`)
      #run(`sed -i 's/.*dash; \*Function.*/---/g' $filename`)
      #run(`sed -i '/>source<\/a>/d' $filename`)
      #run(`sed -i '/>\#<\/a>/d' $filename`)
      _super_cool_example(filename)
    end
  end
end

deploydocs(
  repo = "github.com/thofma/Hecke.jl.git",
  deps = Deps.pip("pymdown-extensions", "pygments", "mkdocs", "python-markdown-math", "mkdocs-material", "mkdocs-cinder"),
  target = "site",
  make = () -> run(`mkdocs build`),
)

#makedocs(
#    modules = [Hecke, Nemo, AbstractAlgebra],
#    clean   = true,
#    format = :html,
#    sitename = "Hecke",
#    doctest = !false,
#)
#
## Hack around to get syntax highlighting working
##cd(joinpath(dirname(pathof(Hecke)), "..", "docs"))
##
##cp("application-f78e5cb881.palette.css", "build/application-f78e5cb881.palette.css", force = true)
##cp("application-e2807e330f.css", "build/application-e2807e330f.css", force = true)
#
#deploydocs(
#     deps = nothing,
##    deps = Deps.pip()#="pygments",
##                    "mkdocs==0.16.3",
##                    "python-markdown-math",
##                    "mkdocs-cinder")=#,
#    repo = "github.com/thofma/Hecke.jl.git",
#    target = "build",
#    make = nothing,
#    osname = "linux",
#    julia = "1.0",
#)
#
## Try out the following deps = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4" ,"python-markdown-math", "pygments", "pymdown-extensions")


text:
@doc raw"""
`ODE`: Ordinary Differential Equation

Defines an initial value problem
```math
\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,
```
with vector field ``v``, initial condition ``q_{0}`` and the solution
``q`` taking values in ``\mathbb{R}^{d}``.

### Parameters

* `vType <: Function`: type of `v`
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `v`: function computing the vector field
* `invariants`: functions for the computation of invariants, either a `NamedTuple` containing the equation's invariants or `NullInvariants`
* `parameters`: type constraints for parameters, either a `NamedTuple` containing the equation's parameters or `NullParameters`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions, either a `AbstractArray` or `NullPeriodicity`

The function `v` providing the vector field must have the interface
```julia
    function v(t, q, v, params)
        v[1] = ...
        v[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, `v` is the
vector which holds the result of evaluating the vector field ``v`` on `t` and
`q`, and params are additional parameters.

### Constructors

```julia
ODE(v, invariants, parameters, periodicity)
ODE(v; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity())
```

"""
struct ODE{vType <: Callable,
           invType <: OptionalInvariants,
           parType <: OptionalParameters,
           perType <: OptionalPeriodicity} <: AbstractEquationODE{invType,parType,perType}

    v::vType

    invariants::invType
    parameters::parType
    periodicity::perType

    function ODE(v, invariants, parameters, periodicity)
        @assert !isempty(methods(v))
        # @assert hasmethod(v, (Real, AbstractArray, AbstractArray, OptionalParameters))

        new{typeof(v), typeof(invariants), typeof(parameters), typeof(periodicity)}(
                v, invariants, parameters, periodicity)
    end
end

ODE(v; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = ODE(v, invariants, parameters, periodicity)

GeometricBase.invariants(equation::ODE) = equation.invariants
GeometricBase.parameters(equation::ODE) = equation.parameters
GeometricBase.periodicity(equation::ODE) = equation.periodicity

hasvectorfield(::ODE) = true

function check_initial_conditions(::ODE, ics::NamedTuple)
    haskey(ics, :q) || return false
    return true
end

function check_methods(equ::ODE, tspan, ics, params)
    applicable(equ.v, tspan[begin], ics.q, zero(ics.q), params) || return false
    return true
end

function datatype(equ::ODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return eltype(ics.q)
end

function arrtype(equ::ODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    typeof(ics.q)
end

_get_v(equ::ODE, params) = (t,q,v) -> equ.v(t, q, v, params)
_get_v̄(equ::ODE, params) = _get_v(equ, params)

_functions(equ::ODE) = (v = equ.v,)
_functions(equ::ODE, params::OptionalParameters) = (v = _get_v(equ, params),)


text:
using Distributions
using Optim

function readAndProcCoinData(filename)
 data = readdlm(filename)
 ndata = length(data)
 heads = zeros(ndata)
 for i=1:ndata
   heads[i] = length(matchall(r"H", data[i]))
 end
 return int64(heads)
end

function sigmoid(x)
 return 1.0./(1.0.+exp(-x))
end

function loglik(x)
 dists = Any[]
 theta = x[1:4]
 cs = x[5:end]
 ntheta = length(theta)
 N = 50
 for i=1:ntheta
   push!(dists, Binomial(N, sigmoid(theta[i])))
 end

 cs = reshape(cs, nprocdata, 4)
 cs = sigmoid(cs)
 cs = cs ./ sum(cs,2)

 ll = 0
 for k=1:length(procdata)
 for i=1:ntheta
   ll += logpdf(dists[i], procdata[k])*cs[k,i]
 end
 end
 ll *= -1
end

function getMeSomeMLE(procdata, theta)
 optimize(loglik, theta, show_trace=true, method = :bfgs, ftol = 1e-12, iterations=5000)

end

# read in and process the data
procdata = readAndProcCoinData("coindata.txt")
nprocdata = length(procdata)

# optimize and return the maximum likelihood estimates
theta = rand(4)
cs = rand(nprocdata*4)
x = [theta; cs]
results = getMeSomeMLE(procdata, x)
x_opt = results.minimum
theta_opt = x_opt[1:4]
cs_opt = x_opt[5:end]
println("Maximum Likelihood Estimates: ", sigmoid(theta_opt))

max_stars_repo_name: haroldsoh/puzzles

max_issues_repo_name: haroldsoh/puzzles

max_forks_repo_name: haroldsoh/puzzles


text:
"""-----------------------------------------------------------------------------
  This routine is the example in ImageFeatures.jl modified to use the HOFASM
  method and diplay the matchings between the points. This method is primarily
  included to provide a template for a visual interface for future
  experimentation.

-----------------------------------------------------------------------------"""
function align_photos(img_1,img_2)

    img1 = load(img_1)
    img2 = load(img_2)

    #grayscale for fast corners
    img1 = Gray.(img1)
    img2 = Gray.(img2)

    #------#  Get image features

    orb_params = ORB(num_keypoints = 100)
    #brief_params = BRIEF(size = 256, window = 10, seed = 123)

    keypoints_1 = Keypoints(fastcorners(img1, 12, 0.4))
    keypoints_2 = Keypoints(fastcorners(img2, 12, 0.4))

    desc_1, ret_keypoints_1 = create_descriptor(img1, orb_params)
    desc_2, ret_keypoints_2 = create_descriptor(img2, orb_params)
    #------#:from imageFeatures.jl demo
    n = length(desc_1)
    m = length(desc_2)

    source_triangles = brute_force_triangles(desc_1)
    target_triangles = brute_force_triangles(desc_2)

    #build tensors
    index_tensor_indices, index_tensor_vals , bases_tensor_indices, bases_tensor_vals =
    build_index_and_bases_tensors(source_triangles, target_triangles, 5.0)

    marg_ten_pairs, kron_time =
        @timed Make_HOFASM_tensor_pairs(index_tensor_indices,bases_tensor_indices,
                                        bases_tensor_vals,n,m)

    x, iteration_time = @timed HOFASM_iterations(marg_ten_pairs,n,m)

    matching = build_assignment(x,use_colsort=true)
    kp_matchings = [[ret_keypoints_1[j],ret_keypoints_2[i]] for (i,j) in matching]

    grid = draw_matches(img1, img2, kp_matchings)

    imshow(grid)
    return grid
end

"""-----------------------------------------------------------------------------
  This is a routine works to run the HOFASM alignment method on a pair of source
  and target embeddings. This method is primarily included to facilitate
  additional experimentation.

-----------------------------------------------------------------------------"""
function align_embeddings(source::Array{Float64,2},target::Array{Float64,2};method="HOFASM")

    n = size(source,1)
    m = size(target,1)

    source_triangles = brute_force_triangles(source)
    target_triangles = brute_force_triangles(target)

    if method == "HOFASM"
        #build tensors
        index_tensor_indices, index_tensor_vals , bases_tensor_indices, bases_tensor_vals =
        build_index_and_bases_tensors(source_triangles, target_triangles, 5.0)

        marg_ten_pairs, kron_time =
            @timed Make_HOFASM_tensor_pairs(index_tensor_indices,bases_tensor_indices,
                                            bases_tensor_vals,n,m)

        x, iteration_time = @timed HOFASM_iterations(marg_ten_pairs,n,m)

    elseif method == "HOM"
        indices, vals = produce_HOM_tensor(source_triangles,target_triangles,m,tol=1e-12)
        marginalized_tensor = sym_mode1_marginalization(indices,vals,n*m)
        x =  Array(reshape(HOM_graduated_assignment(marginalized_tensor,m),n,m)')
    else
        raise()
    end

    return x, build_assignment(x,use_colsort=true)
end

"""-----------------------------------------------------------------------------
  Serialized timing experiments for the HOFASM and HOM methods. Experiments vary
  over different sizes of graphs with a fixed noise level (sigma = 1).

  Inputs
  ------
  * trials - (Int):
     How many trials to run. Extra trial is added and skipped to account
     for compile time.
  * seed - (Bool):
     Whether or not to seed the experiments. Seeds are generated seeding Random.jl
     with 0, and generating as many seeds as there are experiments.
     seeds.
  * method - (String):
     Whether to run 'HOFASM' or 'HOM' for the experiment. sub_method is passed to
     the synthetic_HOFASM to determine which type of marginalization to use.

  Outputs
  -------
  * Experiments - (Dict{Int,Array{1,Floats}}):
    A dictionary linking each of the experiment's n value linked to the runtimes
    of the trials run.
-----------------------------------------------------------------------------"""
function timing_experiments(trials = 100,seed=false;method="HOFASM",sub_method="new")
 #   kron_formation_times = []
    run_times = []
    n_values = [10,20,30,40,60,80,100]
    Experiments = Dict()
    sigma = 1.0

    if seed
        Random.seed!(0)
        seeds = rand(UInt64, length(n_values), trials+1)
    end

    n_exp_j = 1
    for n in n_values

        kron_results = []
        runtime_results = []
#        kron_time = 0.0
        run_time = 0.0
        for i in 1:trials+1
            if method == "HOFASM"
                _ , kt, rt, _ = synthetic_HOFASM(n,sigma,method=sub_method,seed = seeds[n_exp_j,i])
            elseif method == "HOM"
                kt = 0.0
                _, rt, _ = synthetic_HOM(n,sigma,seed = seeds[n_exp_j,i])
            end

            if i == 0
                continue #skip due to compile time
            end
            push!(runtime_results,rt+kt)

        end
        Experiments[n] = runtime_results
             n_exp_j += 1

        println("finished problem size $n")
    end

    return Experiments
end

"""-----------------------------------------------------------------------------
  This function runs timing experiments using the Distributed.jl package.
  Experiments vary over different sizes of graphs with a fixed noise
  level (sigma = 1).

Inputs
------
* nprocs - (Int):
   How many processes have been allocated.
* trials - (Int):
   How many trials to run. Extra trial is added and skipped to account
   for compile time.
* seed - (Bool):
   Whether or not to seed the experiments. Seeds are generated seeding Random.jl
   with 0, and generating as many seeds as there are experiments. Lists of seeds
   for each processer are generated, so the number of processes may affect final
   seeds.
* method - (String):
   Whether to run 'HOFASM' or 'HOM' for the experiment. sub_method is passed to
   the synthetic_HOFASM to determine which type of marginalization to use.

Outputs
-------
* Experiments - (Dict{Int,Array{1,Floats}}):
  A dictionary linking each of the experiment's n value linked to the runtimes
  of the trials run.

-----------------------------------------------------------------------------"""
function distributed_timing_experiments(nprocs::Int,trials = 100,seed=false;method="HOFASM",sub_method="new")

    @assert average_over % nprocs == 0 #must divide for simplicitly

    @everywhere include_string(Main,$(read("HOFASM.jl",String)),"HOFASM.jl")
 #   kron_formation_times = []

    n_values = [10,20,30,40,60,80,100]

    trials_per_process = Int(trials/nprocs)

    if seed
        Random.seed!(0)
        seeds = rand(UInt64, length(n_values),nprocs,trials_per_process + 1)
    end

    Experiments = Dict()
    sigma = 1.0

    #
    #  trials computed on each process
    #
    function sub_tests(trials,n,sigma,method,sub_method;seeds::Array{UInt,1}=nothing)

        if seeds !== nothing
            @assert length(seeds) == trials+1
        end
        results = []
        for i in 1:trials + 1

            if method == "HOFASM"
                _ , kt, rt, _  = synthetic_HOFASM(n,sigma,method=sub_method;seed = seeds[i])
                #_ , kt, rt, _ = synthetic_HOFASM(n,sigma,method=sub_method)
            elseif method == "HOM"
                kt = 0.0
                _, rt, _  = synthetic_HOM(n,sigma;seed = seeds[i])
            end

            if i == 1
                continue #skip run with compile time
            else
                push!(results,kt+rt)
            end

        end
        return results
    end

    n_exp_j = 1
    for n in n_values

        runtime_results = []
        futures = []

        for p in 2:nprocs +1  #leave process 1 as the parent process
            f = @spawnat p sub_tests(trials_per_process,n,sigma,method,sub_method;seeds = seeds[n_exp_j,p-1,:])
            push!(futures,f)
        end

        runtime_results = []
        for future in futures
            r = fetch(future)
            runtime_results = [runtime_results; r] #concatonate on
        end

        Experiments[n] = runtime_results
        n_exp_j += 1
        println("finished problem size $n")
    end

#    kron_formation_times,
    return Experiments
end

"""-----------------------------------------------------------------------------
  Serialized accuracy experiments for the HOFASM and HOM methods. Experiments
  vary over different noise levels (sigma) with a fixed noise level (n = 30).

Inputs
------
* trials - (Int):
   How many trials to average over.
* seed - (Bool):
   Whether or not to seed the experiments. Seeds are generated seeding Random.jl
   with 0, and generating as many seeds as there are experiments. Lists of seeds
   for each processer are generated, so the number of processes may affect final
   seeds.
* method - (String):
   Whether to run 'HOFASM' or 'HOM' for the experiment. sub_method is passed to
   the synthetic_HOFASM to determine which type of marginalization to use.

Outputs
-------
* Experiments - (Dict{Int,Array{1,Floats}}):
  A dictionary linking each of the experiment's sigma value linked to the
  runtimes of the trials run.

-----------------------------------------------------------------------------"""
function accuracy_experiments(trials::Int = 2,seed=false;method="HOFASM",sub_method="new")

    sigmas = [.025,.05,.0725,.1,.125,.15,.175,.2]

    Experiments = Dict()

    n = 30

    if seed
        Random.seed!(0)
        seeds = rand(UInt64,length(sigmas), trials)
    else
        exp_seed = nothing
    end

    sigma_exp_j = 1
    for sigma in sigmas

        colsort_accuracies = []

        for i in 1:trials
            if seed
                exp_seed = seeds[sigma_exp_j,i]
            end

            if method == "HOFASM"
                X , _, _, p  = synthetic_HOFASM(n,sigma,method=sub_method;seed=exp_seed)
            elseif method == "HOM"
                X,_, p  = synthetic_HOM(n,sigmal;seed=seeds[sigma_exp_j,i])
            end

            assignment = build_assignment(X,use_colsort=true)

            colsort_accuracy = sum([1.0 for (i,j) in assignment if p[i] == j])/n

            push!(colsort_accuracies,colsort_accuracy)
        end

        Experiments[sigma] = colsort_accuracies
        sigma_exp_j += 1
    end

    return Experiments
end

"""-----------------------------------------------------------------------------
  This function runs timing experiments using the Distributed.jl package.
  Experiments vary over different noise levels (sigma) with a fixed noise level
  (n = 30).

Inputs
------
* trials - (Int):
   How many trials to run. Extra trial is added and skipped to account
   for compile time.
* seed - (Bool):
   Whether or not to seed the experiments. Seeds are generated seeding Random.jl
   with 0, and generating as many seeds as there are experiments. Lists of seeds
   for each processer are generated, so the number of processes may affect final
   seeds.
* method - (String):
   Whether to run 'HOFASM' or 'HOM' for the experiment. sub_method is passed to
   the synthetic_HOFASM to determine which type of marginalization to use.

Outputs
-------
* Experiments - (Dict{Int,Array{1,Floats}}):
  A dictionary linking each of the experiment's n value linked to the runtimes
  of the trials run.

-----------------------------------------------------------------------------"""
function distributed_accuracy_experiments(average_over::Int = 2,seed::Bool=false;method="HOFASM",sub_method="new")

    @everywhere include_string(Main,$(read("HOFASM.jl",String)),"HOFASM.jl")
    sigmas = [.025,.05,.0725,.1,.125,.15,.175,.2]

    if seed
        Random.seed!(0)
        seeds = rand(UInt64, length(sigmas), average_over)
    else
        seed = nothing
    end

    Experiments = Dict()

    n = 30

    sigma_exp_j = 1
    for sigma in sigmas

        futures = []
#        kron_time = 0.0

        for i in 1:average_over

            if seed
                exp_seed = seeds[sigma_exp_j,i]
            end

            if method == "HOFASM"
                f = @spawn synthetic_HOFASM(n,sigma,method=sub_method;seed=exp_seed)
            elseif method == "HOM"
                f = @spawn synthetic_HOM(n,sigma;seed=seeds[sigma_exp_j,i])
            end
            push!(futures,f)
            #transpose is need to maintain row major format
        end

        colsort_accuracies = []
        for future in futures
            if method == "HOFASM"
                X , _, _, p  = fetch(future)
            elseif method == "HOM"
                X, _, p  = fetch(future)
            end

            assignment = build_assignment(X,use_colsort=true)
            colsort_accuracy = sum([1.0 for (i,j) in assignment if p[i] == j])/n

            push!(colsort_accuracies,colsort_accuracy)
        end

        Experiments[sigma] = colsort_accuracies
        sigma_exp_j += 1
    end

    return Experiments
end

function outlier_experiments(average_over::Int = 2;method="HOFASM")

    sigma = .1
    max_outliers = 20

    colsort_accuracies = []
    colsort_variances = []

    n = 20

    for outliers = 0:2:max_outliers

#        kron_time = 0.0
        test_argsort_accuracies = []
        test_colsort_accuracies = []

        for i in 1:average_over
            if method =="HOM"
                X , _, p = synthetic_HOM(n,sigma,outliers)
            else
                X , _, _, p = synthetic_HOFASM(n,sigma,outliers)
            end

            assignment = build_assignment(X,use_colsort=true)
            colsort_accuracy = sum([1.0 for (i,j) in assignment if p[i] == j])/(n+outliers)

            push!(test_colsort_accuracies,colsort_accuracy)
        end

        push!(colsort_variances,var(test_colsort_accuracies))
        push!(colsort_accuracies,sum(test_colsort_accuracies)/average_over)
    end

    return colsort_variances, colsort_accuracies
end

function scaling_experiments(average_over::Int = 2; method ="HOFASM")

    sigma = .05
    outliers = 5
    scalings = range(1,stop=3,length=10)

    colsort_accuracies = []
    colsort_variances = []

    n = 20

    for scale in scalings

        test_colsort_accuracies = []

        for i in 1:average_over

            if method =="HOM"
                X , _, p = synthetic_HOM(n,sigma,outliers,scale)
            else
                X , _, _, p = synthetic_HOFASM(n,sigma,outliers,scale)
            end

            assignment = build_assignment(X,use_colsort=true)
            colsort_accuracy = sum([1.0 for (i,j) in assignment if p[i] == j])/(n+outliers)

            push!(test_colsort_accuracies,colsort_accuracy)

        end

        push!(colsort_variances,var(test_colsort_accuracies))
        push!(colsort_accuracies,sum(test_colsort_accuracies)/average_over)
    end

    return colsort_variances, colsort_accuracies
end

#=       ---------------Currently not functioning-----------------

# displays the matching results of the source and target points in blue and red respectively
# outliers have a reduced alpha value
function display_synthetic_result(source,target,n,matching)

#    fig = figure("pyplot_scatterplot",figsize=(10,10))
#    ax = PyPlot.axes()

    PyPlot.scatter(source[1:n,1],source[1:n,2],c = "blue")
    PyPlot.scatter(source[n:end,1],source[n:end,2],c = "blue",alpha=0.6)

    PyPlot.scatter(target[1:n,1],target[1:n,2],c = "green")
    PyPlot.scatter(target[n:end,1],target[n:end,2],c = "green",alpha=0.6)

    for (i,j) in matching
        if i == j
            PyPlot.plot(source[i,:], target[j,:],c="black")
        else
            PyPlot.plot(source[i,:], target[j,:],c="red")
        end
    end

    PyPlot.show()
end

function run_and_display_synthetic_exp(n,sigma,outliers,scale)

    X, _,_,source_points, target_points = synthetic_problem_new(n,sigma,outliers,scale)
    display_synthetic_result(source_points,target_points,n,build_assignment(X,use_colsort=true))

end

=#

# Helper Functions from
#    https://www.juliabloggers.com/image-stitching-part-2/

# this function takes the two images and concatenates them horizontally.
# to horizontally concatenate, both images need to be made the same
# vertical size
function pad_display(img1, img2)
    img1h = size(img1, 1)
    img2h = size(img2, 1)
    mx = max(img1h, img2h);
    hcat(vcat(img1, zeros(RGB{Float64},
                max(0, mx - img1h), size(img1, 2))),
        vcat(img2, zeros(RGB{Float64},
                max(0, mx - img2h), size(img2, 2))))
end

function draw_matches(img1, img2, matches)
    # instead of having grid = [img1 img2], we'll use the new
    # pad_display() function
    grid = pad_display(parent(img1), parent(img2));
    offset = CartesianIndex(0, size(img1, 2));
    for m in matches
        draw!(grid, LineSegment(m[1], m[2] + offset))
    end
    grid
end


text:
module Evolution

import Base: shuffle!, getindex
using SimpleStructs

include("common.jl")
include("types.jl")
include("genome.jl")
include("phenotype.jl")
include("genomes/arraygenome.jl")
include("phenotypes/minimizer.jl")

export

    Random,
    RandomNormal,
    Elitism,

    pmap!,

    AbstractGenome,
    IndexableGenome,
    randswap!,

    AbstractPhenotype,
    genome,
    fitness!,
    fitness,
    best,
    worst,

    crossover,
    crossover!,
    mutate!,

    ArrayGenome,
    MinimizerPhenotype

end # module


text:
macro run_exps(db)
    #addprocs(1)
    #n = procs()[end]
    #db = fetch(@spawnat n $(quote
    _results = gensym()
    all_implementations = gensym()
    functions = gensym()
    parameters = gensym()
    results = gensym()
    current_dir = gensym()
    x = gensym()
    esc(quote
        using JuliaDB
        using Query

        # Extract the implementations
        $all_implementations = Int[]
        $x = @from i in db.runs_table begin
            @from j in db.experiments_table
            @where j.Execute && i.Experiment_ID == j.ID && !i.Executed && !i.Executing
            @select i.Implementation_ID
            @collect
        end
        $all_implementations = unique(x)

        # Extract the functions, parameters and results
        $functions = Dict{Int, Symbol}()
        $parameters = Dict{Int, Vector{Tuple{Symbol, Int}}}()
        $results = Dict{Int, Vector{Int}}()
        for implementation in $all_implementations
            $parameters[implementation] = Tuple{Symbol,Int}[]
            $results[implementation] = Int[]

            i = db.implementations_table[implementation]
            $functions[implementation] = i.Function
            if i.Repository != "" && i.Commit != "" && i.Commit != ComputExp.HEAD(i.Repository)
                ComputExp.STASH(i.Repository)
                ComputExp.CHECKOUT(i.Repository, i.Commit)
            end
            if i.Repository != ""
                eval(parse("using $(i.Repository)"))
            end
            include(i.Script)

            params = @from j in db.implementation_parameters_table begin
                @where j.Implementation_ID == implementation
                @select {j.Parameter_symbol, j.ID}
                @collect
            end
            for j in params
                push!($parameters[implementation], (j.Parameter_symbol, j.ID))
            end

            # Retrieves the results of the implementation, order is assumed the same order returned by the function
            $results[implementation] = @from j in db.implementation_results_table begin
                @where j.Implementation_ID == implementation
                @select j.ID
                @collect
            end
        end

        $current_dir = pwd()

        # Run the experiments
        for exper in db.experiments_table
            if exper.Execute
                # Get the run ids of the experiment
                exper_runs = @from i in db.runs_table begin
                    @where i.Experiment_ID == exper.ID
                    @select i.ID
                    @collect
                end

                # Change the directory
                cd(exper.Directory)

                # Run the run
                for rn in exper_runs
                    r = db.runs_table[rn]
                    if !r.Executed && !r.Executing
                        column(db.runs_table, :Executing)[rn] = true
                        problem_features = Dict(@from i in db.problem_features_table begin
                            @where i.Problem_ID == r.Problem_ID
                            @select {i.Feature_symbol, i.Feature_value}
                            @collect
                        end)

                        run_parameters = Dict(@from i in db.run_implementation_parameter_table begin
                            @from j in db.implementation_parameters_table
                            @where i.Run_ID == r.ID && i.Implementation_parameter_ID == j.ID
                            @select {j.Parameter_symbol, i.Value}
                            @collect
                        end)

                        func = db.implementations_table[r.Implementation_ID].Function
                        $_results = eval(func)(; problem = problem_features, parameters = run_parameters)

                        column(db.runs_table, :Executed)[rn] = true
                        column(db.runs_table, :Executing)[rn] = false

                        for (rid, rv) in zip($results[r.Implementation_ID], $_results)
                            resultid = ComputExp.newid(db.run_implementation_result_table)
                            ComputExp.add_row!(db.run_implementation_result_table, @NT(ID = resultid, Implementation_result_ID = rid, Run_ID = r.ID, Value = rv))
                        end
                    end
                end
                #column(db.experiments_table, :Execute)[exper.ID] = false
            end
        end
        cd($current_dir)
    #end))
    #rmprocs(n)
    end)
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbtransp(srf::NURBS2D{I,F}
             )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}

Transpose a NURBS surface, by swapping U and V directions.

Utility function that transposes a NURBS surface, by swapping U and V
directions. NURBS curves cannot be transposed.

# Arguments:
- `srf`: NURBS surface, see [`nrbmak`](@ref).

# Output:
- `tsrf`: NURBS surface with U and V diretions transposed.

# Examples:
```julia
julia> tsrf = nrbtransp(srf)
```
"""
function nrbtransp(srf::NURBS2D{I,F}
                 )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}
return nrbmak(permutedims(srf.coefs,[1,3,2]),reverse(srf.knots));
end # nrbtransp


text:
a = 5
b = 2a^2 + 30a + 9

using PyPlot
x = range(0,stop=5,length=101)
y = cos.(2x .+ 5)
plot(x, y, linewidth=2.0, linestyle="--")
title("a nice cosinus")
xlabel("x axis")
ylabel("y axis")


text:
using ReservoirComputing, CellularAutomata

const input = ones(2, 10)
const output = zeros(2, 10)
const g = 6
const rule = 90

reca = RECA(input, DCA(rule);
    generations = g,
    input_encoding = RandomMapping(6, 10))

output_layer = train(reca, output, StandardRidge(0.001))
prediction = reca(Predictive(input), output_layer)
final_pred = convert(AbstractArray{Int}, prediction .> 0.5)
@test final_pred == output

rm1 = RandomMapping(6, 10)
rm2 = RandomMapping(6, expansion_size = 10)
rm3 = RandomMapping(permutations=6, expansion_size=10)
@test rm1 == rm2 == rm3

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
# Julia wrapper for header: /usr/include/scip/pub_history.h
# Automatically generated using Clang.jl wrap_c

function SCIPhistoryGetVSIDS(history, dir)
    ccall((:SCIPhistoryGetVSIDS, libscip), Cdouble, (Ptr{SCIP_HISTORY}, SCIP_BRANCHDIR), history, dir)
end

function SCIPhistoryGetCutoffSum(history, dir)
    ccall((:SCIPhistoryGetCutoffSum, libscip), Cdouble, (Ptr{SCIP_HISTORY}, SCIP_BRANCHDIR), history, dir)
end

function SCIPvaluehistoryGetNValues(valuehistory)
    ccall((:SCIPvaluehistoryGetNValues, libscip), Cint, (Ptr{SCIP_VALUEHISTORY},), valuehistory)
end

function SCIPvaluehistoryGetHistories(valuehistory)
    ccall((:SCIPvaluehistoryGetHistories, libscip), Ptr{Ptr{SCIP_HISTORY}}, (Ptr{SCIP_VALUEHISTORY},), valuehistory)
end

function SCIPvaluehistoryGetValues(valuehistory)
    ccall((:SCIPvaluehistoryGetValues, libscip), Ptr{Cdouble}, (Ptr{SCIP_VALUEHISTORY},), valuehistory)
end


text:
# Some utility functions

#================================ Compensated Summation =================================#
"""
    comp_sum(sum_value,sum_error,addend)

Kahan (1965) compensated summation for `T<:Real`.
...
# Arguments
- `sum_value::Real` : Current value of the sum.
- `sum_error::Real` : Error from truncation/rounding accumulated from prior steps.
- `addend::Real` : New value to be added to the sum.
...
"""
function comp_sum(sum_value::T,sum_error::T,addend::T) where {T <: Real}
    tmp = zero(T)
    sum_error += addend
    tmp = sum_value + sum_error
    sum_error = (sum_value - tmp) + sum_error
    sum_value = tmp
    return sum_value,sum_error
end

"""
    comp_sum_matrix!(sum_value,sum_error,addend)

Kahan (1965) compenstated summation for `::Matrix{<:Real}`.
...
# Arguments
- `sum_value::Matrix{<:Real}` : Current value of the sum.
- `sum_error::Matrix{<:Real}` : Error from truncation/rounding accumulated from prior steps.
- `addend::Matrix{<:Real}` : New value to be added to the sum.
...
"""
function comp_sum_matrix!(sum_value::Array{T,2},sum_error::Array{T,2},addend::Array{T,2}) where {T <: Real}
    tmp = zero(T)
    @inbounds for i in eachindex(sum_value)
        sum_error[i] += addend[i]
        tmp = sum_value[i] + sum_error[i]
#        sum_error[i] = (sum_value[i] - tmp) + sum_error[i]
        sum_error[i] += sum_value[i] - tmp
        sum_value[i] = tmp
    end
    return
end

"""
  Copies portion of Jacobian for multiplication
"""
     function copy_submatrix!(s::State{T},d::Derivatives{T},indi::Int64,indj::Int64,sevn::Int64)  where T <: AbstractFloat
       # Pick out indices for bodies i & j:
       @inbounds for k2=1:sevn, k1=1:7
           d.jac_tmp1[k1,k2] = s.jac_step[ indi+k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           d.jac_err1[k1,k2] = s.jac_error[indi+k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           d.jac_tmp1[7+k1,k2] = s.jac_step[ indj+k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           d.jac_err1[7+k1,k2] = s.jac_error[indj+k1,k2]
       end
       # Copy current time derivatives for multiplication purposes:
       @inbounds for k1=1:7
           d.dqdt_tmp1[  k1] = s.dqdt[indi+k1]
       end
       @inbounds for k1=1:7
           d.dqdt_tmp1[7+k1] = s.dqdt[indj+k1]
       end
     return
     end

"""
  Copies back:
"""
     function ypoc_submatrix!(s::State{T},d::Derivatives{T},indi::Int64,indj::Int64,sevn::Int64)  where T <: AbstractFloat
       # Copy back to the Jacobian:
       @inbounds for k2=1:sevn, k1=1:7
           s.jac_step[ indi+k1,k2]=d.jac_tmp1[k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           s.jac_error[indi+k1,k2]=d.jac_err1[k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           s.jac_step[ indj+k1,k2]=d.jac_tmp1[7+k1,k2]
       end
       @inbounds for k2=1:sevn, k1=1:7
           s.jac_error[indj+k1,k2]=d.jac_err1[7+k1,k2]
       end
       # Copy back time derivatives:
       @inbounds for k1=1:7
           s.dqdt[indi+k1] = d.dqdt_ij[  k1]
       end
       @inbounds for k1=1:7
           s.dqdt[indj+k1] = d.dqdt_ij[7+k1]
       end
    return
    end
#================================ Integrator Utilities ==================================#

function cubic1(a::T, b::T, c::T) where {T <: Real}
    a3 = a*third
    Q = a3^2 - b*third
    R = a3^3 + 0.5*(-a3*b + c)
    R2 = R^2; Q3=Q^3
    if R2 < Q3
        # println("Error in cubic solver ",R^2," ",Q^3)
        return -c/b
    else
        A = -sign(R)*cbrt(abs(R) + sqrt(R2 - Q3))
        if A == 0.0
            B = 0.0
        else
            B = Q/A
        end
        x1 = A + B - a3
        return x1
    end
    return
end

function G3(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    if gamma < gc
        return G3_series(gamma,beta,sqb)
    else
#        sqb = sqrt(abs(beta))
        if beta >= 0
            return (gamma-sin(gamma))/(sqb*beta)
        else
            return (gamma-sinh(gamma))/(sqb*beta)
        end
    end
end

function G3_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    #epsilon = eps(gamma)
    # Computes G_3(\beta,s) using a series tailored to the precision of s.
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = one(T)
    g3 = one(T)
    g31 = 2g3
    g32 = 2g3
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(g3)
    while true
        g32 = g31
        g31 = g3
        n += 1
        term *= x2/((2n+3)*(2n+2))
        g3 += term
        iter +=1
        if iter >= ITMAX || g3 == g32 || g3 == g31
            break
        end
    end
#    g3 *= gamma^3/(6*sqrt(abs(beta^3)))
    g3 *= -x2*gamma/(6*beta*sqb)
    return g3::T
end

function H1(gamma::T,beta::T;gc=convert(T,0.5)) where {T <: Real}
    #x=sqrt(abs(beta))*s
    if gamma < gc
        return H1_series(gamma,beta)
    else
        if beta >= 0
            return (4sin(0.5*gamma)^2 -gamma*sin(gamma))/beta^2
        else
            return (-4sinh(0.5*gamma)^2+gamma*sinh(gamma))/beta^2
        end
    end
end

function H1_series(gamma::T,beta::T) where {T <: Real}
    # Computes H_1(\beta,s) using a series tailored to the precision of s.
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = one(T)
    h1 = one(T)
    h11 = 2h1
    h12 = 2h1
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h1)
    while true
        h12 = h11
        h11 = h1
        n += 1
        term *= x2*(n+1)
        term /= (2n+4)*(2n+3)*n
        h1 += term
        iter +=1
        if iter >= ITMAX || h1 == h12 || h1 == h11
            break
        end
    end
#    h1 *= gamma^4/(12*beta^2)
    h1 *= x2^2/(12*beta^2)
    return h1::T
end

function H2(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    #x=sqb*s
    if gamma < gc
        return H2_series(gamma,beta,sqb)
    else
#        sqb = sqrt(abs(beta))
        if beta >= 0
            return (sin(gamma)-gamma*cos(gamma))/(sqb*beta)
        else
            return (sinh(gamma)-gamma*cosh(gamma))/(sqb*beta)
        end
    end
end

function H2_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    # Computes H_2(\beta,s) using a series tailored to the precision of s.
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = one(T)
    h2 = one(T)
    h21 = 2h2
    h22 = 2h2
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h2)
    while true
        h22 = h21
        h21 = h2
        n += 1
        term *= x2
        term /= (4n+6)*n
        h2 += term
        iter += 1
        if iter >= ITMAX || h2 == h22 || h2 == h21
            break
        end
    end
#    h2 *= gamma^3/(3*sqrt(abs(beta^3)))
    h2 *= -x2*gamma/(3*beta*sqb)
    return h2::T
end

function H3(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    # This is H_3 = G_1 G_2 - 3 G_3:
    if gamma < gc
        return H3_series(gamma,beta,sqb)
    else
        if beta >= 0
#            return (4*sin(gamma)-sin(gamma)*cos(gamma)-3*gamma)/(beta*sqrt(beta))
            return (4*sin(gamma)-sin(gamma)*cos(gamma)-3*gamma)/(beta*sqb)
        else
#            return (4*sinh(gamma)-sinh(gamma)*cosh(gamma)-3*gamma)/(beta*sqrt(-beta))
            return (4*sinh(gamma)-sinh(gamma)*cosh(gamma)-3*gamma)/(beta*sqb)
        end
    end
end

function H3_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    # Computes H_3(\beta,s) using a series tailored to the precision of gamma:
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = convert(T,1//30)
    h3 = convert(T,1//10)
    h31 = 2h3
    h32 = 2h3
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h3)
    four2n = one(T)*4
    while true
        h32 = h31
        h31 = h3
        n += 1
        term *= x2
        term /= (2n+4)*(2n+5)
        four2n *= 4
#        h3 += term*(4^(n+1)-1)
        h3 += term*(four2n-1)
        iter += 1
        if iter >= ITMAX || h3 == h32 || h3 == h31
            break
        end
    end
#    h3 *= -gamma^5/(beta*sqrt(abs(beta)))
    h3 *= -x2^2*gamma/(beta*sqb)
    return h3::T
end

function H5(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    # This is G_1 G_2 -(2+G_0) G_3 = H_2 - 2 G_3:
    if gamma < gc
        return H5_series(gamma,beta,sqb)
    else
        if beta >= 0
#            return (3sin(gamma)-2gamma-gamma*cos(gamma))/(beta*sqrt(beta))
            return (3sin(gamma)-2gamma-gamma*cos(gamma))/(beta*sqb)
        else
#            return (3sinh(gamma)-2gamma-gamma*cosh(gamma))/(beta*sqrt(-beta))
            return (3sinh(gamma)-2gamma-gamma*cosh(gamma))/(beta*sqb)
        end
    end
end

function H5_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    # Computes H_5(\beta,s) using a series tailored to the precision of gamma:
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = one(T)/60
    h5 = copy(term)
    h51 = 2h5
    h52 = 2h5
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h5)
    while true
        h52 = h51
        h51 = h5
        n += 1
        term *= x2*(n+1)
        term /= (2n+5)*(2n+4)*n
        h5 += term
        iter += 1
        if iter >= ITMAX || h5 == h52 || h5 == h51
            break
        end
    end
#    h5 *= -gamma^5/(beta*sqrt(abs(beta)))
    h5 *= -x2^2*gamma/(beta*sqb)
    return h5::T
end

function H6(gamma::T,beta::T;gc=convert(T,0.5)) where {T <: Real}
    # This is 2 G_2^2 -3 G_1 G_3:
    if gamma < gc
        return H6_series(gamma,beta)
    else
        if beta >= 0
            return (9-8cos(gamma) -cos(2gamma) -6gamma*sin(gamma))/(2beta^2)
        else
            return (9-8cosh(gamma)-cosh(2gamma)+6gamma*sinh(gamma))/(2beta^2)
        end
    end
end

function H6_series(gamma::T,beta::T) where {T <: Real}
    # Computes H_6(\beta,s) using a series tailored to the precision of gamma:
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = convert(T,1//360)
    h6 = convert(T,1//40)
    h61 = 2h6
    h62 = 2h6
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h6)
    four2n = one(T)*16
    while true
        h62 = h61
        h61 = h6
        n += 1
        term *= x2
        term /= (2n+5)*(2n+6)
        four2n *= 4
#        h6 += term*(4^(n+2)-3*n-7)
        h6 += term*(four2n-3*n-7)
        iter +=1
        if iter >= ITMAX || h6 == h62 || h6 == h61
            break
        end
    end
#    h6 *= gamma^6/(beta*abs(beta))
    h6 *= -x2^3/(beta^2)
    return h6::T
end

function H7(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    # This is H_7 = G_1 G_2 (1-2 G_0) + 3 G_0^2 G_3:
    if gamma < gc
        return H7_series(gamma,beta,sqb)
    else
        if beta >= 0
#            return (3*cos(gamma)*(gamma*cos(gamma)-sin(gamma))+sin(gamma)^3)/(beta*sqrt(beta))
            return (3*cos(gamma)*(gamma*cos(gamma)-sin(gamma))+sin(gamma)^3)/(beta*sqb)
        else
#            return (3*cosh(gamma)*(gamma*cosh(gamma)-sinh(gamma))-sinh(gamma)^3)/(beta*sqrt(-beta))
            return (3*cosh(gamma)*(gamma*cosh(gamma)-sinh(gamma))-sinh(gamma)^3)/(beta*sqb)
        end
    end
end

function H7_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    # Computes H_7(\beta,s) using a series tailored to the precision of gamma:
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = -convert(T,3//20160)*x2
    h7 = convert(T,1//10-11//840*x2)
    h71 = 2h7
    h72 = 2h7
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h7)
    four2n = one(T)*128
    nine2n = one(T)*729
    while true
        h72 = h71
        h71 = h7
        n += 1
        term *= x2
        term /= (2n+6)*(2n+7)
        four2n *= 4
        nine2n *= 9
#        h7 += term*(9^(3+n)-1-(5+2*n)*2^(7+2*n))
        h7 += term*(nine2n-1-(5+2*n)*four2n)
        iter += 1
        if iter >= ITMAX || h7 == h72 || h7 == h71
            break
        end
    end
#    h7 *= gamma^5/(beta*sqrt(abs(beta)))
    h7 *= x2^2*gamma/(beta*sqb)
    return h7::T
end

function H8(gamma::T,beta::T,sqb::T;gc=convert(T,0.5)) where {T <: Real}
    # This is H_8 = G_1 G_2 - 3 G_0 G_3:
    if gamma < gc
        return H8_series(gamma,beta,sqb)
    else
        if beta >= 0
#            return (-3gamma*cos(gamma) +sin(gamma) +sin(2gamma))/(beta*sqrt(beta))
            return (-3gamma*cos(gamma) +sin(gamma) +sin(2gamma))/(beta*sqb)
        else
#            return (-3gamma*cosh(gamma)+sinh(gamma)+sinh(2gamma))/(beta*sqrt(-beta))
            return (-3gamma*cosh(gamma)+sinh(gamma)+sinh(2gamma))/(beta*sqb)
        end
    end
end

function H8_series(gamma::T,beta::T,sqb::T) where {T <: Real}
    # Computes H_8(\beta,s) using a series tailored to the precision of gamma:
    #epsilon = eps(gamma)
    #x2 = -beta*s^2
    x2 = -sign(beta)*gamma^2
    term = convert(T,1//120)
    h8 = convert(T,3//20)
    h81 = 2h8
    h82 = 2h8
    n=0
    iter = 0
    ITMAX = 100
    # Terminate series when required precision reached:
    #while abs(term) > epsilon*abs(h8)
    four2n = one(T)*32
    while true
        h82 = h81
        h81 = h8
        n += 1
        term *= x2
        term /= (2n+4)*(2n+5)
        four2n *= 4
#        h8 += term*(2^(5+2n)-14-6n)
        h8 += term*(four2n-14-6n)
        iter +=1
#        println("H8: ",iter," ",h8," ",gamma," ",beta)
        if iter >= ITMAX || h8 == h82 || h8 == h81
            break
        end
    end
#    h8 *= gamma^5/(beta*sqrt(abs(beta)))
    h8 *= x2^2*gamma/(beta*sqb)
    return h8::T
end

# Faster dot product; assumes 3D vector
@inline function dot_fast(a::Vector{T}, b::Vector{T}) where T<:Real
    a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
end

@inline function dot_fast(a::Vector{T}) where T<:Real
    a[1]*a[1] + a[2]*a[2] + a[3]*a[3]
end


text:
# ============================================================
# Structs of Control Volume
# Solver stores data in arrays of struct (AoS)
# ============================================================

# ------------------------------------------------------------
# General
# ------------------------------------------------------------

"""
$(TYPEDEF)

Control volume with no distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolume{T1,T2,T3,ND} <: AbstractControlVolume
    w::T1
    prim::T2
    sw::T3
end

function ControlVolume(w, prim, sw)
    n = size(w, 1) - 2
    return ControlVolume{typeof(w),typeof(prim),typeof(sw),n}(w, prim, sw)
end

function Base.show(io::IO, ctr::ControlVolume{A,B,C,N}) where {A,B,C,N}
    print(
        io,
        "ControlVolume$(N)D{$A,$B,$C}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
    )
end

"""
$(TYPEDEF)

Control volume with 1 distribution function

# Fields

$(FIELDS)
"""
struct ControlVolume1F{T1,T2,T3,T4,ND} <: AbstractControlVolume
    w::T1
    prim::T1
    sw::T2
    f::T3
    sf::T4
end

function ControlVolume1F(w, prim, sw, f, sf)
    n = size(w, 1) - 2
    return ControlVolume1F{typeof(w),typeof(sw),typeof(f),typeof(sf),n}(w, prim, sw, f, sf)
end

function Base.show(io::IO, ctr::ControlVolume1F{A,B,C,D,N}) where {A,B,C,D,N}
    print(
        io,
        "ControlVolume$(N)D1F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

"""
$(TYPEDEF)

Control volume with 2 distribution functions

# Fields

$(FIELDS)
"""
struct ControlVolume2F{T1,T2,T3,T4,ND} <: AbstractControlVolume
    w::T1
    prim::T1
    sw::T2
    h::T3
    b::T3
    sh::T4
    sb::T4
end

function ControlVolume2F(w, prim, sw, h, b, sh, sb)
    n = size(w, 1) - 2
    return ControlVolume2F{typeof(w),typeof(sw),typeof(h),typeof(sh),n}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end

function Base.show(io::IO, ctr::ControlVolume2F{A,B,C,D,N}) where {A,B,C,D,N}
    print(
        io,
        "ControlVolume$(N)D2F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

struct ControlVolume3F{T1,T2,T3,T4,ND} <: AbstractControlVolume end

struct ControlVolume4F{T1,T2,T3,T4,ND} <: AbstractControlVolume end

"""
$(SIGNATURES)

Construct control volume...
"""
function ControlVolume(W, PRIM, ND::Integer)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = begin
        if ND == 1
            zero(w)
        elseif ND == 2
            slope_array(w)
        end
    end

    return ControlVolume{typeof(w),typeof(prim),typeof(sw),ND}(w, prim, sw)
end

"""
$(SIGNATURES)
"""
function ControlVolume(W, PRIM, F, ND::Integer)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    f = deepcopy(F)
    sw, sf = begin
        if ND == 1
            zero(w), zero(f)
        elseif ND == 2
            slope_array(w), slope_array(f)
        end
    end

    return ControlVolume1F{typeof(w),typeof(sw),typeof(f),typeof(sf),ND}(w, prim, sw, f, sf)
end

"""
$(SIGNATURES)
"""
function ControlVolume(W, PRIM, H, B, ND::Integer)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    h = deepcopy(H)
    b = deepcopy(B)
    sw, sh, sb = begin
        if ND == 1
            zero(w), zero(h), zero(b)
        elseif ND == 2
            slope_array(w), slope_array(h), slope_array(b)
        end
    end

    return ControlVolume2F{typeof(w),typeof(sw),typeof(h),typeof(sh),ND}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end

# ------------------------------------------------------------
# The dimension-dependent structures are in archive mode only
# ------------------------------------------------------------
# 1D
# ------------------------------------------------------------

"""
$(TYPEDEF)

1D control volume with no distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolume1D{A,B} <: AbstractControlVolume1D
    w::A
    prim::B
    sw::A
end

function ControlVolume1D(W, PRIM)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    return ControlVolume1D{typeof(w),typeof(prim)}(w, prim, sw)
end

function Base.show(io::IO, ctr::ControlVolume1D{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
    )
end

"""
$(TYPEDEF)

1D control volume with 1 distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolume1D1F{A,B} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    f::B
    sf::B
end

function ControlVolume1D1F(W::T, PRIM::T, F) where {T}
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    f = deepcopy(F)
    sf = zero(f)

    return ControlVolume1D1F{typeof(w),typeof(f)}(w, prim, sw, f, sf)
end

function Base.show(io::IO, ctr::ControlVolume1D1F{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D1F{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

"""
$(TYPEDEF)

1D control volume with 2 distribution functions

# Fields

$(FIELDS)
"""
struct ControlVolume1D2F{A,B} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h::B
    b::B
    sh::B
    sb::B
end

function ControlVolume1D2F(W::T1, PRIM::T1, H::T2, B::T2) where {T1<:AA,T2<:AA}
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h = deepcopy(H)
    b = deepcopy(B)
    sh = zero(h)
    sb = zero(b)

    return ControlVolume1D2F{typeof(w),typeof(h)}(w, prim, sw, h, b, sh, sb)
end

function Base.show(io::IO, ctr::ControlVolume1D2F{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D2F{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h, b\n",
        "pdf slopes: sh, sb\n",
    )
end

"""
$(TYPEDEF)

1D control volume with 3 distribution functions

# Fields

$(FIELDS)
"""
mutable struct ControlVolume1D3F{A,B,C,D,E} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h0::B
    h1::B
    h2::B
    sh0::B
    sh1::B
    sh2::B

    E::C
    B::C
    ϕ::D
    ψ::D
    lorenz::E
end

# deterministic
function ControlVolume1D3F(
    W::AA{<:Real,2},
    PRIM::AA{<:Real,2},
    H0::AA{<:FN,3},
    H1::AA{<:FN,3},
    H2::AA{<:FN,3},
    E0::AA{<:FN,1},
    B0::AA{<:FN,1},
    L::AA{<:FN,2},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

# Rykov
function ControlVolume1D3F(
    W::AA{<:Real,1},
    PRIM::AA{<:Real,1},
    H0::AA{<:FN,1},
    H1::AA{<:FN,1},
    H2::AA{<:FN,1},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = nothing
    B = nothing
    ϕ = nothing
    ψ = nothing
    lorenz = nothing

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

# stochastic
function ControlVolume1D3F(
    W::AA{<:Real,3},
    PRIM::AA{<:Real,3},
    H0::AA{<:FN,4},
    H1::AA{<:FN,4},
    H2::AA{<:FN,4},
    E0::AA{<:FN,2},
    B0::AA{<:FN,2},
    L::AA{<:FN,3},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = zero(E[1, :]) # here is difference
    ψ = zero(B[1, :])
    lorenz = deepcopy(L)

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

function Base.show(io::IO, ctr::ControlVolume1D3F{A,B,C,D,E}) where {A,B,C,D,E}
    print(
        io,
        "ControlVolume1D3F{$A,$B,$C,$D,$E}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2\n",
        "pdf slopes: sh0, sh1, sh2\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

"""
$(TYPEDEF)

1D control volume with 4 distribution functions

# Fields

$(FIELDS)
"""
mutable struct ControlVolume1D4F{A,B,C,D,E} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h0::B
    h1::B
    h2::B
    h3::B
    sh0::B
    sh1::B
    sh2::B
    sh3::B

    E::C
    B::C
    ϕ::D
    ψ::D
    lorenz::E
end

#--- deterministic ---#
function ControlVolume1D4F(
    W::AA{<:Real,2},
    PRIM::AA{<:Real,2},
    H0::AA{<:FN,2},
    H1::AA{Float64,2},
    H2::AA{Float64,2},
    H3::AA{Float64,2},
    E0::AA{Float64,1},
    B0::AA{Float64,1},
    L::AA{Float64,2},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    h3 = deepcopy(H3)
    sh0 = deepcopy(H0)
    sh1 = deepcopy(H1)
    sh2 = deepcopy(H2)
    sh3 = deepcopy(H3)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        h3,
        sh0,
        sh1,
        sh2,
        sh3,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

#--- uncertainty quantification ---#
function ControlVolume1D4F(
    W::AA{<:Real,3},
    PRIM::AA{<:Real,3},
    H0::AA{<:FN,3},
    H1::AA{Float64,3},
    H2::AA{Float64,3},
    H3::AA{Float64,3},
    E0::AA{Float64,2},
    B0::AA{Float64,2},
    L::AA{Float64,3},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    h3 = deepcopy(H3)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)
    sh3 = zero(H3)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = zero(E[1, :])
    ψ = zero(B[1, :])
    lorenz = deepcopy(L)

    return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        h3,
        sh0,
        sh1,
        sh2,
        sh3,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

function Base.show(io::IO, ctr::ControlVolume1D4F{A,B,C,D,E}) where {A,B,C,D,E}
    print(
        io,
        "ControlVolume1D4F{$A,$B,$C,$D,$E}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2, h3\n",
        "pdf slopes: sh0, sh1, sh2, sh3\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

# ------------------------------------------------------------
# 2D
# ------------------------------------------------------------

"""
$(TYPEDEF)

2D control volume with no distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolume2D{A,B} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B
end

function ControlVolume2D(W, PRIM)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    return ControlVolume2D{typeof(w),typeof(sw)}(w, prim, sw)
end

function Base.show(io::IO, ctr::ControlVolume2D{A,B}) where {A,B}
    print(
        io,
        "ControlVolume2D{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
    )
end

"""
$(TYPEDEF)

2D control volume with 1 distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolume2D1F{A,B,C,D} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B

    f::C
    sf::D
end

function ControlVolume2D1F(W, PRIM, F::AA)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    f = deepcopy(F)
    #sf = zeros(eltype(F), (axes(F)..., Base.OneTo(2)))
    sf = slope_array(F)

    return ControlVolume2D1F{typeof(w),typeof(sw),typeof(f),typeof(sf)}(w, prim, sw, f, sf)
end

function Base.show(io::IO, ctr::ControlVolume2D1F{A,B,C,D}) where {A,B,C,D}
    print(
        io,
        "ControlVolume2D1F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

"""
$(TYPEDEF)

2D control volume with 2 distribution functions

# Fields

$(FIELDS)
"""
struct ControlVolume2D2F{A,B,C,D} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B

    h::C
    b::C
    sh::D
    sb::D
end

function ControlVolume2D2F(W::AA, PRIM::AA, H::AA, B::AA)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    h = deepcopy(H)
    b = deepcopy(B)
    #sh = zeros(eltype(H), (axes(H)..., Base.OneTo(2)))
    #sb = zeros(eltype(B), (axes(B)..., Base.OneTo(2)))
    sh = slope_array(H)
    sb = slope_array(B)

    return ControlVolume2D2F{typeof(w),typeof(sw),typeof(h),typeof(sh)}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end

function Base.show(io::IO, ctr::ControlVolume2D2F{A,B,C,D}) where {A,B,C,D}
    print(
        io,
        "ControlVolume2D2F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h, b\n",
        "pdf slopes: sh, sb\n",
    )
end

"""
$(TYPEDEF)

2D control volume with 3 distribution functions

# Fields

$(FIELDS)
"""
mutable struct ControlVolume2D3F{T2,T3,T4,T5,T6,T7,T8} <: AbstractControlVolume2D
    w::T2
    prim::T2
    sw::T3

    h0::T4
    h1::T4
    h2::T4
    sh0::T5
    sh1::T5
    sh2::T5

    E::T6
    B::T6
    ϕ::T7
    ψ::T7
    lorenz::T8
end

#--- deterministic & stochastic ---#
function ControlVolume2D3F(W::AA, PRIM::AA, H0::AA, H1::AA, H2::AA, E0::AA, B0::AA, L::AA)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D
    sw = slope_array(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    #sh0 = zeros(eltype(H0), (axes(H0)..., Base.OneTo(2)))
    #sh1 = zeros(eltype(H1), (axes(H1)..., Base.OneTo(2)))
    #sh2 = zeros(eltype(H2), (axes(H2)..., Base.OneTo(2)))
    sh0 = slope_array(H0)
    sh1 = slope_array(H1)
    sh2 = slope_array(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume2D3F{
        typeof(w),
        typeof(sw),
        typeof(h0),
        typeof(sh2),
        typeof(E),
        typeof(ϕ),
        typeof(lorenz),
    }(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

#--- Rykov ---#
function ControlVolume2D3F(
    W::AA{<:Real,1},
    PRIM::AA{<:Real,1},
    H0::AA{<:FN,2},
    H1::AA{<:FN,2},
    H2::AA{<:FN,2},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D
    sw = slope_array(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    #sh0 = zeros(eltype(H0), (axes(H0)..., Base.OneTo(2)))
    #sh1 = zeros(eltype(H1), (axes(H1)..., Base.OneTo(2)))
    #sh2 = zeros(eltype(H2), (axes(H2)..., Base.OneTo(2)))
    sh0 = slope_array(H0)
    sh1 = slope_array(H1)
    sh2 = slope_array(H2)

    E = nothing
    B = nothing
    ϕ = nothing
    ψ = nothing
    lorenz = nothing

    return ControlVolume2D3F{
        typeof(w),
        typeof(sw),
        typeof(h0),
        typeof(sh2),
        typeof(E),
        typeof(ϕ),
        typeof(lorenz),
    }(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

function Base.show(
    io::IO,
    ctr::ControlVolume2D3F{T1,T2,T3,T4,T5,T6,T7},
) where {T1,T2,T3,T4,T5,T6,T7}
    print(
        io,
        "ControlVolume2D3F{$T1,$T2,$T3,$T4,$T5,$T6,$T7}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2\n",
        "pdf slopes: sh0, sh1, sh2\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

# ------------------------------------------------------------
# Unstructured cell
# ------------------------------------------------------------

"""
$(TYPEDEF)

Unstructured control volume with no distribution function

# Fields

$(FIELDS)
"""
mutable struct ControlVolumeUS{E,F,A,B,C} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::B
    sw::C
end

function ControlVolumeUS(N, X::T, DX::T, W, PRIM) where {T<:Union{Real,AV}}
    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    return ControlVolumeUS{typeof(n),typeof(x),typeof(w),typeof(prim),typeof(sw)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
    )
end

"""
$(TYPEDEF)

Unstructured control volume with 1 distribution function

# Fields

$(FIELDS)
"""
struct ControlVolumeUS1F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::A
    sw::B

    f::C
    sf::D
end

function ControlVolumeUS1F(N, X, DX, W, PRIM, F::AA{T}) where {T}
    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    f = deepcopy(F)
    #sf = zeros(eltype(F), axes(F)..., length(N[1]))
    sf = ifelse(
        length(N[1]) == 2,
        slope_array(F; reduction = true),
        slope_array(F; reduction = false),
    )

    return ControlVolumeUS1F{typeof(n),typeof(x),typeof(w),typeof(sw),typeof(f),typeof(sf)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
        f,
        sf,
    )
end

"""
$(TYPEDEF)

Unstructured control volume with 2 distribution functions

# Fields

$(FIELDS)
"""
struct ControlVolumeUS2F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::A
    sw::B

    h::C
    b::C
    sh::D
    sb::D
end

function ControlVolumeUS2F(N, X, DX, W::T1, PRIM::T1, H::T2, B::T2) where {T1<:AA,T2<:AA}
    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    h = deepcopy(H)
    b = deepcopy(B)
    #sh = zeros(eltype(H), axes(H)..., length(N[1]))
    #sb = zeros(eltype(B), axes(B)..., length(N[1]))
    sh = ifelse(
        length(N[1]) == 2,
        slope_array(H; reduction = true),
        slope_array(H; reduction = false),
    )
    sb = ifelse(
        length(N[1]) == 2,
        slope_array(B; reduction = true),
        slope_array(B; reduction = false),
    )

    return ControlVolumeUS2F{typeof(n),typeof(x),typeof(w),typeof(sw),typeof(h),typeof(sh)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end


text:
function copper_plate(ps_m::CanonicalModel, expression::Symbol, bus_count::Int64, time_range::UnitRange{Int64})

    devices_netinjection = _remove_undef!(ps_m.expressions[expression])

    ps_m.constraints[:CopperPlateBalance] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, time_range)

    for t in time_range
        ps_m.constraints[:CopperPlateBalance][t] = JuMP.@constraint(ps_m.JuMPmodel, sum(ps_m.expressions[expression].data[1:bus_count,t]) == 0)
    end

    return

end


text:
# Generate 3D data for the overthrust model
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: May 2017
#

using PyCall, PyPlot, HDF5, opesciSLIM.TimeModeling, opesciSLIM.SLIM_optim, SeisIO

# Load starting model
n,d,o,m0 = h5open("../data/overthrust_2D_initial_model.h5","r") do file
	read(file, "n", "d", "o", "m0")
end

# Set up model structure
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)	# need n,d,o as tuples and m0 as array

# Bound constraints
vmin = ones(Float32,model0.n) + 0.3f0
vmax = ones(Float32,model0.n) + 5.5f0

# Slowness squared [s^2/km^2]
mmin = vec((1f0./vmax).^2)
mmax = vec((1f0./vmin).^2)

# Load data
block = segy_read("../data/overthrust_2D.segy")
dobs = joData(block)

# Set up wavelet
src_geometry = Geometry(block; key="source", segy_depth_key="SourceDepth")
wavelet = source(src_geometry.t[1],src_geometry.dt[1],0.008f0)	# 8 Hz wavelet
q = joData(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
srand(1)	# set seed of random number generator
fevals = 20
batchsize = 20
fvals = zeros(21)

count = 0
function objective_function(x)
	model0.m = reshape(x,model0.n);

	# select batch
	idx = randperm(dobs.nsrc)[1:batchsize]
	dsub = subsample(dobs,idx)
	qsub = subsample(q,idx)

	# fwi function value and gradient
	fval,grad = fwi_objective(model0,qsub,dsub)
	grad = reshape(grad,model0.n); grad[:,1:21] = 0f0	# reset gradient in water column to 0.
	grad = .125f0*grad/maximum(abs.(grad))

	global count; count += 1; fvals[count] = fval
    return fval, vec(grad)
end

# Bound projection
ProjBound(x) = median([mmin x mmax],2)

# FWI with SPG
options = spg_options(verbose=3, maxIter=fevals, memory=3)
x, fsave, funEvals= minConf_SPG(objective_function, vec(model0.m), ProjBound, options)

# Save final velocity model, function value and history
h5open("fwi_overthrust_2D_SPG.h5", "w") do file
	write(file, "x", sqrt.(1f0./reshape(x,model0.n)), "fsave", fsave, "fhistory",fvals)
end

# Plot results
figure(); imshow(sqrt.(1f0./m0)',cmap="jet",vmin=1.5,vmax=6.); title("Initial model")
figure(); imshow(sqrt.(1f0./reshape(x,model0.n))',cmap="jet",vmin=1.5,vmax=6.); title("FWI")
figure(); plot(fvals); title("Function value")

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
#Special named numbers and symbols

export bellnum,
    catalannum,
    lobbnum,
    narayana,
    fibonaccinum,
    jacobisymbol,
    lassallenum,
    legendresymbol,
    lucasnum,
    stirlings1,
    stirlings2,
    dealnnoy

"""
    bellnum(n)

Compute the ``n``th Bell number.
"""
function bellnum(n::Integer)
    if n < 0
        throw(DomainError(n))
    elseif n < 2
        return 1
    end
    list = Vector{BigInt}(undef, n)
    list[1] = 1
    for i = 2:n
        for j = 1:i - 2
            list[i - j - 1] += list[i - j]
        end
        list[i] = list[1] + list[i - 1]
    end
    return list[n]
end

"""
    catalannum(n)

Compute the ``n``th Catalan number.
"""
function catalannum(bn::Integer)
    if bn < 0
        throw(DomainError(bn, "n must be nonnegative"))
    else
        n = BigInt(bn)
    end
    div(binomial(2*n, n), (n + 1))
end

"""
    lobbnum(m,n)

Compute the Lobb number `L(m,n)`, or the generalised Catalan number given by ``\\frac{2m+1}{m+n+1} \\binom{2n}{m+n}``.
Wikipedia : https://en.wikipedia.org/wiki/Lobb_number
"""
function lobbnum(bm::Integer,bn::Integer)
    if !(0 <= bm <= bn)
        throw(DomainError("m and n must be non-negative"))
    else
        m = BigInt(bm)
        n = BigInt(bn)
    end
    div(binomial(2*n, m + n)*(2*m + 1), (m + n + 1))
end

"""
    narayana(n,k)

Compute the Narayana number `N(n,k)`` given by ``\\frac{1}{n}\\binom{n}{k}\\binom{n}{k-1}``
Wikipedia : https://en.wikipedia.org/wiki/Narayana_number
"""
function narayana(bn::Integer,bk::Integer)
    if !(1 <= bk <= bn)
        throw(DomainError("Domain is 1 <= k <= n"))
    else
        n = BigInt(bn)
        k = BigInt(bk)
    end
    div(binomial(n, k)*binomial(n, k - 1) , n)
end

function fibonaccinum(n::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    end
    z = Ref{BigInt}(0)
    ccall((:__gmpz_fib_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
    return z[]
end

function jacobisymbol(a::Integer, b::Integer)
    ba = Ref{BigInt}(a)
    bb = Ref{BigInt}(b)
    return ccall((:__gmpz_jacobi, :libgmp), Cint, (Ref{BigInt}, Ref{BigInt}), ba, bb)
end

"""
    lassallenum(n)

Compute the ``n``th entry in Lassalle's sequence, OEIS entry A180874.
"""
function lassallenum(m::Integer)
    A = ones(BigInt, m)
    for n = 2:m
        A[n] = (-1)^(n-1) * (catalannum(n) + sum(j->(-1)^j*binomial(2n-1, 2j-1)*A[j]*catalannum(n-j), 1:n-1))
    end
    A[m]
end

function legendresymbol(a::Integer, b::Integer)
    ba = Ref{BigInt}(a)
    bb = Ref{BigInt}(b)
    return ccall((:__gmpz_legendre, :libgmp), Cint, (Ref{BigInt}, Ref{BigInt}), ba, bb)
end

function lucasnum(n::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    end
    z = Ref{BigInt}(0)
    ccall((:__gmpz_lucnum_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
    return z[]
end

function stirlings1(n::Int, k::Int, signed::Bool=false)
    if signed == true
        return (-1)^(n - k) * stirlings1(n, k)
    end

    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == k == 0
        return 1
    elseif n == 0 || k == 0
        return 0
    elseif n == k
        return 1
    elseif k == 1
        return factorial(n-1)
    elseif k == n - 1
        return binomial(n, 2)
    elseif k == n - 2
        return div((3 * n - 1) * binomial(n, 3), 4)
    elseif k == n - 3
        return binomial(n, 2) * binomial(n, 4)
    end

    return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
end

function stirlings2(n::Int, k::Int)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == k == 0
        return 1
    elseif n == 0 || k == 0
        return 0
    elseif k == n - 1
        return binomial(n, 2)
    elseif k == 2
        return 2^(n-1) - 1
    end

    return k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
end

"""
Delannoy number  describes the number of paths from the southwest corner (0, 0) of a rectangular grid
to the northeast corner (m, n), using only single steps north, northeast, or east
"""
function dealnnoy(m::Integer,n::Integer)
    if m<0
        throw(DomainError(m, "m must be nonnegative"))
    end
    if n<0
        throw(DomainError(n, "n must be nonnegative"))
    end
    A = ones(BigInt, m+1, n+1)
    for i = 1:m+1
        A[i,1] = 1
    end
    for i = 1:n+1
        A[1,i] = 1
    end
    for i = 2:m+1
        for j = 2:n+1
            A[i,j] = A[i-1,j] + A[i-1,j-1] + A[i,j-1]
        end
    end
    return(A[m+1,n+1])
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 98c9bbd2-aac5-4c90-ac0c-d8d935f5cdaf
# For WGLMakie
begin
	using Pkg
	Pkg.add("JSServe")
	using JSServe
	Page()
end

# ╔═╡ 0a622047-c238-49c3-bdf0-3248d0f0d261
# Just for me
begin
	using Revise
	cd("../")
	Pkg.activate("@v1.6")
	Pkg.develop(path="./AllenNeuropixels.jl")
	cd("./AllenNeuropixels.jl")
	Pkg.activate("./")
end

# ╔═╡ 2e3d037c-932d-4d5a-afdb-351e836bdfb2
using WGLMakie

# ╔═╡ 79b92c17-cc5f-4ca2-8b08-f6015729a9a9
using DataFrames

# ╔═╡ 1a3eaeb7-db19-4bc3-b61a-84dd9caad009
using Statistics

# ╔═╡ c1397e8e-3ac2-4f49-8d88-cc1400a3b93e
using FileIO

# ╔═╡ 766a8af5-4c89-4fe7-883d-d960ef91edfd
md"""
# Allen Neuropixels
_Accessing the Allen Neuropixels visual coding dataset in Julia_
"""

# ╔═╡ bea5de79-6e8a-42d8-ab76-bae8e3c23747
md"""
## Background

Details on neuropixels and the visual coding dataset can be found in the Allen SDK [docs](https://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html), the [white-paper](https://dl.dropbox.com/s/tav6rft6iyd173k/neuropixels_visual_coding_-_white_paper_v10.pdf) or the [cheat-sheet](https://dl.dropbox.com/s/9v1x9eibmtb8fwg/neuropixels_cheat_sheet_nov_2019.pdf)
"""

# ╔═╡ f9ac9f6e-f129-4542-81a8-36e6cef9857a
md"## Required packages"

# ╔═╡ c48bd78e-aab0-49c0-b137-567c208b8aa1
md"""
Interfacing with the Python [Allen SDK](https://github.com/AllenInstitute/AllenSDK) is handled by the [AllenNeuropixels](https://github.com/brendanjohnharris/AllenNeuropixels.jl) package, which uses the [Makie](https://github.com/JuliaPlots/Makie.jl) package for plotting (with GPU acceleration and interactivity). We will use the WGLMakie as the backend, as well as DataFrames and Statistics for working with tables.
"""

# ╔═╡ c445ccf4-cf10-43b9-9c01-4051abc400ba
#begin
#	Pkg.add(url="https://github.com/brendanjohnharris/AllenNeuropixels.jl#main")
#	Pkg.add("WGLMakie")
#	Pkg.add("DataFrames")
#	Pkg.add("Statistics")
#	Pkg.add("FileIO")
#end

# ╔═╡ d8fd4541-08a5-4ace-a998-769771c976e8
import AllenNeuropixels as AN

# ╔═╡ 5bfaefae-11a5-4567-8b83-d583f03a75a8
md"""
## Choosing a session
The Allen neuropixels visual coding dataset is subdivided into sessions. A session contains data from one murine subject across a recording interval during which it was shown a variety of visual stimuli (such as natural scenes, drift gratings and gabors). These data include local field potentials (LFP) and spike times from each of the 374 data channels on, usually, six neuropixels probes inserted around the visual cortex. The LFP data is also downsampled by a half in time and a quarter over channels.

The entire neuropixels visual coding dataset contains dozens of sessions and is hundreds of megabytes in size, so we will first choose one session (a few gigabytes of data) instead of performing group-level analyses. To produce a dataframe of session details:
"""

# ╔═╡ 754f94e0-ccb2-4dc0-a534-ae94dd02bc02
session_table = AN.getsessiontable()

# ╔═╡ a2f81f74-36ed-42d4-89d5-828327c67318
md"""We'll pick a session that has six probes:"""

# ╔═╡ efd3ef52-4acb-4ffd-b794-1dfc4d9819c8
function numprobes(sessionid)
	session_table[session_table[!, :id].==sessionid, :probe_count][1]
end

# ╔═╡ e30e349c-6ad7-402b-90d1-7720b85a9c2c
md"And look for sessions that have probes intersecting the target regions listed in the white paper:"

# ╔═╡ d2d20884-0fcd-4ac7-8a14-dbf936445c3b
function targetintersections(sessionid)
	s = session_table[session_table.id.==sessionid, :ecephys_structure_acronyms][1]
	targets = ["VISp", "VISl", "VISrl", "VISal", "VISpm", "VISam", "CA1", "CA3", "DG", "SUB", "ProS", "LGd", "LP", "APN"];
	structures = eachmatch(r"'(.*?)'", s)
	structures = [replace(x.match, r"'"=>s"") for x ∈ structures]
	targetsintersected = [x ∈ structures for x ∈ targets]
	return mean(targetsintersected)
end;

# ╔═╡ 99fe41d6-661d-4b9f-a853-2f32ace53d72
md"""
We can filter down sessions in a few more ways. Firstly, the visaul coding data is divided into two stimulus sets. To summuraise the white paper:
- The **Brain Observatory 1.1** stimulus set:
    - Gabor patches appearing randomly on a 9 x 9 grid
    - Dark or light flashes
    - Drift gratings in 8 directions and 5 temporal frequencies, with 15 repeats per condition
    - Static gratings at 6 orientations, 5 spatial frequencies, and 4 phases
    - 118 natural images
    - Two natural movies from Touch of Evil; a 30 second clips repeated 20 times and a 120 second clip repeated 10 times

- The **Functional Connectivity** stimulus set:
    - Gabor patches
    - Dark or light flashes
    - Drift gratings in 4 directions and one temporal frequency with 75-80 repeats
    - Drift gratings in 4 directions and 9 contrasts
    - One natural movie shown 60 times plus 20 repeats of a temporally shuffled version
    - A dot motion stimulus of approximately 200 white dots on a gray background moving at one of 7 speeds in four different directions

In summary, the Brain Observatory dataset contains a greater variety of stimuli but a smaller number of repeats than the Functional Connectivity dataset. **We'll look at sessions in the Brain Observatory dataset.**
"""

# ╔═╡ 82ce6d0f-b60b-41f2-bdce-c6ecf545bf65
md"""
Next, we can inspect the unit quality metrics of each session. Three metric criteria have recommended values:
- `amplitude_cutoff_maximum = 0.1`: Limits the approximate false negative rate in spike detection, calculated by assuming the distribution of amplitudes is symmetric.
- `presence_ratio_minimum = 0.9`: Ensures proportion of 100 recording sub-intervals that have at least one spike detection is above 90% (i.e. the unit has not drifted and the sorting algorithm is correctly identifying spikes).
- `isi_violations_maximum = 0.5`: Limits the number of units that record from multiple neurons. Inter-spike interval (ISI) violations are detections of spikes that occur during the typical refractory period of a single neuron, so are likely to originate from a second neuron.

To access a dataframe of metrics for the brain observatory dataset with the default metric filters:
"""

# ╔═╡ 5d47cb91-d8c7-41df-9778-c9a77266ba93
# This can take a while
metrics = AN.getunitanalysismetricsbysessiontype("brain_observatory_1.1")

# ╔═╡ 7d00cd72-a2c0-45a9-b777-b347286f7390
md"Some session-level metrics are:"

# ╔═╡ fd03139c-4bd1-4f26-a8c3-46c3feefd9c5
session_metrics = combine(groupby(metrics, :ecephys_session_id),
		:ecephys_session_id => numprobes∘unique => :num_probes,
		:ecephys_session_id => targetintersections∘unique => :target_intersections,
		:has_lfp_data=>all,
		:genotype => (x->all(isequal.(x, ("wt/wt",)))) => :is_wt,
		:max_drift=>median∘skipmissing,
		:d_prime=>median∘skipmissing,
		:isolation_distance=>median∘skipmissing,
		:silhouette_score=>median∘skipmissing,
		:snr=>median∘skipmissing)

# ╔═╡ 51eccd32-d4e9-4330-99f5-fca0f8534e43
md"""
Of particular importance are the first five columns. The ideal session will have 6 probes, LFP data for all units and for now we wish to choose a wild type mouse. It should also contain data for probes that intersect with the major regions of interest highlighted in the white paper. The maximum drift, which measures mean (and possibly variance) nonstationarity unrelated to the visual stimulus (such as probe motion relative to the brain) should also be minimised. We can also look for probes that intersect all of the target regions, given in the white paper as the visual areas, the hippocampal formation, the thalamus and the midbrain:
"""

# ╔═╡ 927605f4-0b59-4871-a13f-420aadedd487
oursession = subset(session_metrics,
						:num_probes => ByRow(==(6)),
						:target_intersections => ByRow(>(0.9)),
						:is_wt => ByRow(==(true)),
						:has_lfp_data_all => ByRow(==(true)),
						:max_drift_median_skipmissing => ByRow(<(25.0)))

# ╔═╡ 23d8fad0-5c51-4056-8df3-fd850db7b560
md"""
For a smaller quantity of data we can also pick just one probe from this session: `probeC: 769322749` intersects a few interesting regions, like the primary visual cortex, the hippocampus and the thalamus.
"""

# ╔═╡ 99e2e888-c1cf-4370-b1b4-eccc6cbda7de
probeid = 769322749

# ╔═╡ d552f132-6f90-4aba-9e05-1e7e20929756
md"""
## Inspecting data

To mimic the Allen SDK's interface we can use a custom type wrapping a session obect:
"""

# ╔═╡ c3093ce3-7b73-49d4-8ce8-aaea4b49b685
session = AN.Session(oursession.ecephys_session_id[1])

# ╔═╡ 3ec0b209-2287-4fcf-be04-688bd4fb327a
sessionid = AN.getid(session)

# ╔═╡ ebe8631b-4248-4b5a-94e2-071c0b560747
md"The six probes for this session are:"

# ╔═╡ c77034a4-a6b4-4970-bdc8-9a6b30cbb687
probes = AN.getprobes(session)

# ╔═╡ 12d66773-c567-4e5a-a119-4fa5e06ec98c
md"We are only interested in channels located within the brain, so the `missing` structure ids are removed and the channels for this session become:"

# ╔═╡ 5f944e90-5232-4508-9a3a-c678e6804104
channels = subset(AN.getchannels(session), :ecephys_structure_id=>ByRow(!ismissing))

# ╔═╡ 423edd0f-60ed-4d9a-b84a-47e47e560ae2
md"""
We can then plot the probe locations on the reference atlas. Note that the colored regions correspond to the target structures, and the straight probes have been warped into the common coordinate framework:
"""

# ╔═╡ 12f8e03b-4ea3-4211-a9b1-f8432dfae3a9
#s = AN.Plots.plotreferencevolume(session; dostructures=true,
#								ids=:targets,
#								show_axis=false,
#								shading=true); rotate_cam!(s, Vec3f0(0, 2.35, 0)); s

# To save you from waiting for the structure masks to download:
@html_str read(download("https://dl.dropbox.com/s/se2doygr56ox8hs/probelocations.html?dl=0"), String)
# This renders best in firefox

# ╔═╡ b9c2e65c-03cf-45d7-9309-b601f486c62b
md"""
The LFP data for our probe can be accessed as (warning, it is slow):
"""

# ╔═╡ f0241126-912b-4863-9d4e-917af1426602
LFP = AN.getdownsampledlfp(session, probeid)

# ╔═╡ d0a301a0-038c-4827-8683-e9d8807186ea
md"And sorted by depth with:"

# ╔═╡ a5fca30e-531e-4b09-97e9-a762059dc66c
sortedLFP = AN.sortbydepth(session, probeid, LFP)[:, end:-1:1]

# ╔═╡ 39ae3db9-b799-4389-80e7-e898e4e88a84
fig = AN.Plots.neuroslidingcarpet(sortedLFP[1:5000, :]; resolution=(800, 1600))

# ╔═╡ 87ab5ec5-67a9-413d-8789-6a8b7113de65
md"""
There are still channels with little signal outside of the brain and at the very surface of the isocortex. To remove these, downsample by 10x and use data from only one stimulus epoch:
"""

# ╔═╡ da9fa8b0-afcf-4f3d-bd6d-856b69ab6d28
begin
	depthcutoff = 200 # μm
	times = AN.getepochs(session, "natural_movie_three")[1, :]
	times = [times.start_time, times.stop_time]
	inbrainLFP = AN.getlfp(session, probeid; inbrain=depthcutoff, times)
end

# ╔═╡ Cell order:
# ╠═98c9bbd2-aac5-4c90-ac0c-d8d935f5cdaf
# ╠═0a622047-c238-49c3-bdf0-3248d0f0d261
# ╟─766a8af5-4c89-4fe7-883d-d960ef91edfd
# ╟─bea5de79-6e8a-42d8-ab76-bae8e3c23747
# ╟─f9ac9f6e-f129-4542-81a8-36e6cef9857a
# ╟─c48bd78e-aab0-49c0-b137-567c208b8aa1
# ╠═c445ccf4-cf10-43b9-9c01-4051abc400ba
# ╠═d8fd4541-08a5-4ace-a998-769771c976e8
# ╠═2e3d037c-932d-4d5a-afdb-351e836bdfb2
# ╠═79b92c17-cc5f-4ca2-8b08-f6015729a9a9
# ╠═1a3eaeb7-db19-4bc3-b61a-84dd9caad009
# ╠═c1397e8e-3ac2-4f49-8d88-cc1400a3b93e
# ╟─5bfaefae-11a5-4567-8b83-d583f03a75a8
# ╠═754f94e0-ccb2-4dc0-a534-ae94dd02bc02
# ╟─a2f81f74-36ed-42d4-89d5-828327c67318
# ╠═efd3ef52-4acb-4ffd-b794-1dfc4d9819c8
# ╟─e30e349c-6ad7-402b-90d1-7720b85a9c2c
# ╠═d2d20884-0fcd-4ac7-8a14-dbf936445c3b
# ╟─99fe41d6-661d-4b9f-a853-2f32ace53d72
# ╟─82ce6d0f-b60b-41f2-bdce-c6ecf545bf65
# ╠═5d47cb91-d8c7-41df-9778-c9a77266ba93
# ╟─7d00cd72-a2c0-45a9-b777-b347286f7390
# ╠═fd03139c-4bd1-4f26-a8c3-46c3feefd9c5
# ╟─51eccd32-d4e9-4330-99f5-fca0f8534e43
# ╠═927605f4-0b59-4871-a13f-420aadedd487
# ╟─23d8fad0-5c51-4056-8df3-fd850db7b560
# ╟─99e2e888-c1cf-4370-b1b4-eccc6cbda7de
# ╟─d552f132-6f90-4aba-9e05-1e7e20929756
# ╠═c3093ce3-7b73-49d4-8ce8-aaea4b49b685
# ╟─3ec0b209-2287-4fcf-be04-688bd4fb327a
# ╟─ebe8631b-4248-4b5a-94e2-071c0b560747
# ╠═c77034a4-a6b4-4970-bdc8-9a6b30cbb687
# ╟─12d66773-c567-4e5a-a119-4fa5e06ec98c
# ╠═5f944e90-5232-4508-9a3a-c678e6804104
# ╟─423edd0f-60ed-4d9a-b84a-47e47e560ae2
# ╟─12f8e03b-4ea3-4211-a9b1-f8432dfae3a9
# ╟─b9c2e65c-03cf-45d7-9309-b601f486c62b
# ╠═f0241126-912b-4863-9d4e-917af1426602
# ╟─d0a301a0-038c-4827-8683-e9d8807186ea
# ╠═a5fca30e-531e-4b09-97e9-a762059dc66c
# ╠═39ae3db9-b799-4389-80e7-e898e4e88a84
# ╠═87ab5ec5-67a9-413d-8789-6a8b7113de65
# ╠═da9fa8b0-afcf-4f3d-bd6d-856b69ab6d28


text:
using Printf

include("INC.jl")

function create_lattice_vars()
    LatVecs = zeros(3,3)
    A = 5.13
    LatVecs[1,:] = [A, A, 0.0]
    LatVecs[2,:] = [A, 0.0, A]
    LatVecs[3,:] = [0.0, A, A]
    lattice_vars = LatticeVars( LatVecs )
    return lattice_vars
end

function create_atomic_vars(lattice_vars)

    maxatoms = 200
    maxspecies = 8
    atposl = zeros(3,maxatoms,maxspecies)

    nspecies = 2
    natoms = [1,1]

    # species 1, atom 1
    atposl[:,1,1] = [0.0, 0.0, 0.0]
    # species 2, atom 1
    atposl[:,1,2] = [0.25, 0.25, 0.25]

    atomic_vars = AtomicVars(nspecies, natoms, atposl, lattice_vars)
end

function main()

    latt_vars = create_lattice_vars()
    atm_vars = create_atomic_vars(latt_vars)

    Nspecies = 2
    atsp_vars = AtomicSpeciesVars(Nspecies)
    mtr_vars = MuffinTinRadialVars(Nspecies)
    apwlo_vars = APWLOVars(Nspecies, mtr_vars.lmaxapw)

    readspecies!(1, "DATA_species/Si.in", atsp_vars, mtr_vars, apwlo_vars)
    readspecies!(2, "DATA_species/Pt.in", atsp_vars, mtr_vars, apwlo_vars)

    init_zero!( mtr_vars )

    println("before nrsp = ", atsp_vars.nrsp)
    println("atsp_vars.rsp = ", size(atsp_vars.rsp))

    checkmt!( latt_vars, atm_vars, atsp_vars.spsymb, mtr_vars )
    genrmesh!( atm_vars, atsp_vars, mtr_vars )
    init_packed_mtr!(mtr_vars)

    println("after nrsp = ", atsp_vars.nrsp)
    println("atsp_vars.rsp = ", size(atsp_vars.rsp))

    for is in 1:Nspecies

        println()
        println("species: ", is)
        println()

        println("rsp = ")
        println(atsp_vars.rsp[1:5,is])

        println("wmrt = ")
        println(mtr_vars.wrmt[1:5,is])

        println("wcmrt = ")
        println(mtr_vars.wrcmt[1:5,is])

        #println("wprmt 1 = ")
        #println(mtr_vars.wprmt[1,1:5,is])

        println("wprcmt 1 = ")
        println(mtr_vars.wprcmt[1,1:5,is])
        println("wprcmt 2 = ")
        println(mtr_vars.wprcmt[2,1:5,is])
        println("wprcmt 3 = ")
        println(mtr_vars.wprcmt[3,1:5,is])
        println("wprcmt 4 = ")
        println(mtr_vars.wprcmt[4,1:5,is])

        #println("wprmt 2 = ")
        #println(mtr_vars.wprmt[2,1:5,1])

        #println("wprmt 3 = ")
        #println(mtr_vars.wprmt[3,1:5,1])

        #println("wprmt 4 = ")
        #println(mtr_vars.wprmt[4,1:5,1])
    end

    println("lmmaxi   = ", mtr_vars.lmmaxi)
    println("npmtmax  = ", mtr_vars.npmtmax)
    println("npmt     = ", mtr_vars.npmt)
    println("npcmtmax = ", mtr_vars.npcmtmax)
    println("npcmt    = ", mtr_vars.npcmt)

end

@time main()