text:
function rhs!(du::AbstractArray{<:Any,5}, u, t,
    mesh::CurvedMesh, equations,
    initial_condition, boundary_conditions, source_terms,
    dg::DG, cache)
  # Reset du
  @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, mesh, equations,
                                                                dg.volume_integral, dg, cache)

  # Calculate interface fluxes
  @timeit_debug timer() "interface flux" calc_interface_flux!(u, mesh, equations, dg, cache)

  # Calculate surface integrals
  @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  @timeit_debug timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache)

  # Calculate source terms
  @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, mesh::CurvedMesh, equations,
                               volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis
  @unpack contravariant_vectors = cache.elements

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)

      flux1 = flux(u_node, 1, equations)
      flux2 = flux(u_node, 2, equations)
      flux3 = flux(u_node, 3, equations)

      # Compute the contravariant flux by taking the scalar product of the
      # first contravariant vector Ja^1 and the flux vector
      contravariant_flux1 = (contravariant_vectors[1, 1, i, j, k, element] * flux1 +
                             contravariant_vectors[1, 2, i, j, k, element] * flux2 +
                             contravariant_vectors[1, 3, i, j, k, element] * flux3)

      for ii in eachnode(dg)
        integral_contribution = derivative_dhat[ii, i] * contravariant_flux1
        add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, k, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # second contravariant vector Ja^2 and the flux vector
      contravariant_flux2 = (contravariant_vectors[2, 1, i, j, k, element] * flux1 +
                             contravariant_vectors[2, 2, i, j, k, element] * flux2 +
                             contravariant_vectors[2, 3, i, j, k, element] * flux3)

      for jj in eachnode(dg)
        integral_contribution = derivative_dhat[jj, j] * contravariant_flux2
        add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, k, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # third contravariant vector Ja^3 and the flux vector
      contravariant_flux3 = (contravariant_vectors[3, 1, i, j, k, element] * flux1 +
                             contravariant_vectors[3, 2, i, j, k, element] * flux2 +
                             contravariant_vectors[3, 3, i, j, k, element] * flux3)

      for kk in eachnode(dg)
        integral_contribution = derivative_dhat[kk, k] * contravariant_flux3
        add_to_node_vars!(du, integral_contribution, equations, dg, i, j, kk, element)
      end
    end
  end

  return nothing
end

function calc_interface_flux!(u::AbstractArray{<:Any,5}, mesh::CurvedMesh{3},
                              equations, dg::DG, cache)
  @unpack elements = cache
  @unpack surface_flux = dg

  @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2, 3)"

    # Interfaces in x-direction (`orientation` = 1)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, mesh, equations, dg, cache)

    # Interfaces in x-direction (`orientation` = 2)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, mesh, equations, dg, cache)

    # Interfaces in x-direction (`orientation` = 3)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[3, element],
                         element, 3, u, mesh, equations, dg, cache)
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u,
                                      mesh::CurvedMesh{3}, equations, dg::DG, cache)
  @unpack surface_flux = dg
  @unpack contravariant_vectors = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for j in eachnode(dg), i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, j, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, j, right_element)

      # First contravariant vector Ja^1 as SVector
      normal = get_contravariant_vector(1, contravariant_vectors, 1, i, j, right_element)
    elseif orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), j, left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          j, right_element)

      # Second contravariant vector Ja^2 as SVector
      normal = get_contravariant_vector(2, contravariant_vectors, i, 1, j, right_element)
    else # orientation == 3
      u_ll = get_node_vars(u, equations, dg, i, j, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, j, 1,          right_element)

      # Third contravariant vector Ja^3 as SVector
      normal = get_contravariant_vector(3, contravariant_vectors, i, j, 1, right_element)
    end

    flux = surface_flux(u_ll, u_rr, normal, equations)

    for v in eachvariable(equations)
      surface_flux_values[v, i, j, right_direction, left_element] = flux[v]
      surface_flux_values[v, i, j, left_direction, right_element] = flux[v]
    end
  end

  return nothing
end

function apply_jacobian!(du::AbstractArray{<:Any,5}, mesh::CurvedMesh, equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      factor = -cache.elements.inverse_jacobian[i, j, k, element]

      for v in eachvariable(equations)
        du[v, i, j, k, element] *= factor
      end
    end
  end

  return nothing
end


text:
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Imports
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

using JuliaSub # if it's a registered package

# for working without installing JuliaSub it as a package
#include("../src/JuliaSub.jl")
#using Main.JuliaSub

using Test

using MacroTools
using Multisets

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Aux values and functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

const TEST_FILES_DIR = "test-files"

# String → String
testFilePath(path :: String) = joinpath(@__DIR__, TEST_FILES_DIR, path)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

include("utils.jl")
include("lb-analysis.jl")

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

ext: jl

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
using Darwin
import Darwin: uniform_mutation
import Base.isless
import JSON.json

mutable struct NeurodevoInd <: Individual
    genes::Array{Array{Float64}}
    fitness::Array{Float64}
    func::Function
    seed::Int64
end

function NeurodevoInd(cfg::Dict)
    neuro_cfg = Config("cfg/evo.yaml")
    chromos = make_chromos!(neuro_cfg; cinds=[4, 5, 9, 10])
    genes = Array{Array{Float64}}(undef, 0)
    push!(genes, rand(4))
    for c in chromos
        push!(genes, c.genes)
    end
    fitness = -Inf*ones(cfg["n_fitness"])
    NeurodevoInd(genes, fitness)
end

function isless(i1::NeurodevoInd, i2::NeurodevoInd)
    minimum(i1.fitness) < minimum(i2.fitness)
end

function json(i::NeurodevoInd)
    json(Dict("genes"=>i.genes, "fitness"=>i.fitness))
end

function uniform_mutation(parent::NeurodevoInd; m_rate=0.1)
    genes = Array{Array{Float64}}(undef, 0)
    for gene in parent.genes
        cgene = rand(length(gene))
        inds = rand(length(gene)) .>= m_rate
        cgene[inds] = gene[inds]
        push!(genes, cgene)
    end
    NeurodevoInd(genes, deepcopy(parent.fitness))
end


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: jl_mfcc2sdc.jl
#
#  Decription:
#       Transform MFCCs into SDCs
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-24   wm              Initial version
#
################################################################################

using ArgParse
using HDF5
using DSP.TFFilter

function deltas{T<:FloatingPoint}(x::Matrix{T}, w::Int=9)
    """
    Code by David van Leeuwen (davidavdav) taken from https://github.com/davidavdav/MFCC.jl/blob/master/src/mfccs.jl
    """
    (nr, nc) = size(x)

    if nr == 0 || w <= 1
        return x
    end

    hlen = ifloor(w / 2)
    w = 2hlen + 1                 # make w odd
    win = [convert(T, hlen):-1:-hlen]
    xx = vcat(repmat(x[1, :], hlen, 1), x, repmat(x[end, :], hlen, 1)) ## take care of boundaries
    norm = 3 / (hlen * w * (hlen + 1))
    return norm * filt(TFFilter(win, [1.]), xx)[2hlen + (1:nr), :]
end

function sdc{T<:FloatingPoint}(x::Matrix{T}, n::Int=7, d::Int=1, p::Int=3, k::Int=7)
    """
    Code by David van Leeuwen (davidavdav) taken from https://github.com/davidavdav/MFCC.jl/blob/master/src/mfccs.jl
    """
    nx, nfea = size(x)
    n <= nfea || error("Not enough features for n")

    hnfill = (k - 1) * p / 2
    nfill1, nfill2 = ifloor(hnfill), iceil(hnfill)

    xx = vcat(zeros(eltype(x), nfill1, n), deltas(x[:,1:n], 2d + 1), zeros(eltype(x), nfill2, n))
    y = zeros(eltype(x), nx, n * k)
    for (dt, offset) = zip(0:p:p * k - 1, 0:n:n * k - 1)
        y[:, offset + (1:n)] = xx[dt + (1:nx), :]
    end

    return y
end

function transform(
    NMFCC::Int,
    N::Int, D::Int, P::Int, K::Int,
    MFCC_H5_FNAME::String, MFCC_H5_DSET::String, OUT_H5_FNAME::String, OUT_H5_DSET::String)

    fid = h5open(MFCC_H5_FNAME, "r")
    fullX = fid[MFCC_H5_DSET]

    const NCOLS = size(fullX, 2)
    @assert length(fullX[:, 1]) % NMFCC == 0
    const NFRAMES = div(length(fullX[:, 1]), NMFCC)

    ofid = h5open(OUT_H5_FNAME, "w")
    odset = d_create(ofid, OUT_H5_DSET, datatype(Float64), dataspace(N * K * NFRAMES, NCOLS))

    for coli = 1:NCOLS
        Xi = reshape(fullX[:, coli], (NFRAMES, NMFCC))
        sdc_ = sdc(Xi, N, D, P, K)
        odset[:, coli] = vec(sdc_)

        if coli % 100 == 0
            println(coli)
        end
    end

    close(fid)
    close(ofid)

end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        #"--opt1"
        #    help = "an option with an argument"
        "--nmfcc", "-m"
            help = "number of MFCCs"
            arg_type = Int
            default = 7
        "--sdc-n", "-n"
            help = "SDC 'N' parameter"
            arg_type = Int
            default = 7
        "--sdc-d", "-d"
            help = "SDC 'D' parameter"
            arg_type = Int
            default = 1
        "--sdc-p", "-p"
            help = "SDC 'P' parameter"
            arg_type = Int
            default = 3
        "--sdc-k", "-k"
            help = "SDC 'K' parameter"
            arg_type = Int
            default = 7
        "--h5-mfcc-file"
            help = "input H5 store for MFCC data"
            arg_type = String
            required = true
        "--h5-mfcc-dset"
            help = "input H5 dataset for MFCC data"
            arg_type = String
            required = true
        "--h5-sdc-file"
            help = "output H5 store for SDC data"
            arg_type = String
            required = true
        "--h5-sdc-dset"
            help = "output H5 dataset for SDC data"
            arg_type = String
            required = true

        #"arg1"
        #    help = "a positional argument"
        #    required = true
    end

    return parse_args(s)
end

function main(args)

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

    transform(
        parsed_args["nmfcc"],
        parsed_args["sdc-n"], parsed_args["sdc-d"], parsed_args["sdc-p"], parsed_args["sdc-k"],
        parsed_args["h5-mfcc-file"], parsed_args["h5-mfcc-dset"], parsed_args["h5-sdc-file"], parsed_args["h5-sdc-dset"])

end

if ~isinteractive()
    main(ARGS)
end


text:
function get_devices(sys::PowerSystems.PowerSystem,category::Type{PowerSystems.ThermalGen})
    return sys.generators.thermal
end
function get_devices(sys::PowerSystems.PowerSystem,category::Type{PowerSystems.RenewableGen})
    return sys.generators.renewable
end
function get_devices(sys::PowerSystems.PowerSystem,category::Type{PowerSystems.HydroGen})
    return sys.generators.hydro
end
function get_devices(sys::PowerSystems.PowerSystem,category::Type{PowerSystems.PowerSystems.ElectricLoad})
    return sys.loads
end

function constructservice!(m::JuMP.AbstractModel, service::PowerSystems.StaticReserve, category_formulation::Type{PS.RampLimitedReserve},devices::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, sys::PowerSystems.PowerSystem; args...)

    dev_set = Array{NamedTuple{(:device,:formulation),Tuple{PowerSystems.PowerSystemDevice,DataType}}}([])

    for category in devices
        if category != nothing
            D = get_devices(sys,category.device)
            for d in D
                if d in service.contributingdevices
                    push!(dev_set,(device=d,formulation=category.formulation))
                end
            end
        end
    end

    if !isempty(dev_set)

        p_rsv =  reservevariables(m, dev_set, sys.time_periods)
        m =  PS.reserves(m, dev_set, service, sys.time_periods)

    end

    return m

end


text:
module LazyDistances

using Reexport
using Distances

include("CachedViews/CachedViews.jl")
@reexport using .CachedViews
using .CachedViews: CachedView, AbstractCacheStrategy, NullCacheArray, LocalWindowCacheArray, make_cache

export PairwiseDistance

# extend Distances.result_type
# https://github.com/JuliaStats/Distances.jl/pull/185
result_type(f, a::Type, b::Type) = typeof(f(oneunit(a), oneunit(b)))
result_type(f::Distances.PreMetric, a::Type, b::Type) = Distances.result_type(f, a, b) # ambiguity patch
result_type(f, A::AbstractArray, b::AbstractArray) = result_type(f, eltype(A), eltype(B))
result_type(dist::Distances.PreMetric, a, b) = Distances.result_type(dist, a, b)

include("utilities.jl")
include("pairwise.jl")

end #module


text:
export A2CGAELearner

"""
    A2CGAELearner(;kwargs...)
# Keyword arguments
- `approximator`, an [`ActorCritic`](@ref) based [`NeuralNetworkApproximator`](@ref)
- `γ::Float32`, reward discount rate.
- `λ::Float32`, lambda for GAE-lambda
- `actor_loss_weight::Float32`
- `critic_loss_weight::Float32`
- `entropy_loss_weight::Float32`
"""
Base.@kwdef mutable struct A2CGAELearner{A<:ActorCritic} <: Any
    approximator::A
    γ::Float32
    λ::Float32
    max_grad_norm::Union{Nothing,Float32} = nothing
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    update_freq::Int
    update_step::Int = 0
    # for logging
    actor_loss::Float32 = 0.0f0
    critic_loss::Float32 = 0.0f0
    entropy_loss::Float32 = 0.0f0
    loss::Float32 = 0.0f0
    norm::Float32 = 0.0f0
end

Functors.functor(x::A2CGAELearner) =
    (app = x.approximator,), y -> @set x.approximator = y.app

(learner::A2CGAELearner)(env::MultiThreadEnv) =
    learner.approximator.actor(send_to_device(device(learner), state(env))) |> send_to_host

function RLBase.update!(learner::A2CGAELearner, t::CircularArraySARTTrajectory)
    length(t) == 0 && return  # in the first update, only state & action is inserted into trajectory
    learner.update_step += 1
    if learner.update_step % learner.update_freq == 0
        _update!(learner, t)
    end
end

function _update!(learner::A2CGAELearner, t::CircularArraySARTTrajectory)
    n = length(t)

    AC = learner.approximator
    to_device(x) = send_to_device(device(AC), x)
    γ = learner.γ
    λ = learner.λ
    w₁ = learner.actor_loss_weight
    w₂ = learner.critic_loss_weight
    w₃ = learner.entropy_loss_weight

    S = t[:state] |> to_device
    # (state_size..., n_thread * update_step)
    states_flattened = select_last_dim(S, 1:n) |> flatten_batch

    actions =
        t[:action] |>
        x ->
            select_last_dim(x, 1:n) |> flatten_batch |> a -> CartesianIndex.(a, 1:length(a))

    rollout_values =
        S |> flatten_batch |> AC.critic |> x -> reshape(x, :, n + 1) |> send_to_host

    advantages = generalized_advantage_estimation(
        t[:reward],
        rollout_values,
        γ,
        λ;
        dims = 2,
        terminal = t[:terminal],
    )

    gains = to_device(advantages + select_last_dim(rollout_values, 1:n))

    advantages = advantages |> flatten_batch |> to_device

    ps = Flux.params(AC)
    gs = gradient(ps) do
        logits = AC.actor(states_flattened)
        probs = softmax(logits)
        log_probs = logsoftmax(logits)
        log_probs_select = log_probs[actions]
        values = AC.critic(states_flattened)
        advantage = vec(gains) .- vec(values)
        actor_loss = -mean(log_probs_select .* advantages)
        critic_loss = mean(advantage .^ 2)
        entropy_loss = -sum(probs .* log_probs) * 1 // size(probs, 2)
        loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss
        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
            learner.entropy_loss = entropy_loss
            learner.loss = loss
        end
        loss
    end

    if !isnothing(learner.max_grad_norm)
        learner.norm = clip_by_global_norm!(gs, ps, learner.max_grad_norm)
    end

    update!(AC, gs)
end

RLCore.check(::QBasedPolicy{<:A2CGAELearner}, ::MultiThreadEnv) = nothing


text:
using AugmentedGaussianProcesses
using Distributions
using Distances
using Random
ν = 3.0

testfunction(x) = (1 + x / ν)^(-0.5 * (ν + 1))
testfunction(x) = sech(sqrt(x / 2))
# testfunction(x) = (1+sqrt(3*x))*exp(-sqrt(3*x))
# testfunction(x) = (1+sqrt(5*x)+5/3*x)*exp(-sqrt(5*x))
testdist = AGP.LaplaceTransformDistribution(
    testfunction, 0.0, AGP.BromwichInverseLaplace(19, 1, 10, 5)
)
@btime rand($testdist)
@btime rand(Gamma((ν + 1) / 2, 1 / ν))
rand(testdist)
# @btime AGP.draw(AugmentedGaussianProcesses.PolyaGammaDist(),1,0)
evaluate(KLDivergence(), rand(testdist, 10000), rand(Gamma((ν + 1) / 2, 1 / ν), 10000))
# laptrans(testdist,n=1)[1]
p = display(histogram(rand(testdist, 10000); bins=range(0, 3; length=100), normalize=true))
histogram!(
    rand(Gamma((ν + 1) / 2, 1 / ν), 10000);
    alpha=0.5,
    bins=range(0, 10; length=100),
    normalize=true,
)
# histogram!([AGP.draw(AugmentedGaussianProcesses.PolyaGammaDist(),1,0) for i in 1:10000],alpha=0.5,bins=range(0,10,length=100),normalize=true)


text:
import MechanicalSketch
import MechanicalSketch: color_with_lumin, empty_figure, background, sethue
import MechanicalSketch: PALETTE, O, EM, FS, finish, x_y_iterators_at_pixels, latexify
import MechanicalSketch: dimension_aligned, settext, arrow, setfont
import MechanicalSketch: ComplexQuantity, generate_complex_potential_source, @import_expand, string_polar_form
import MechanicalSketch: place_image, normalize_binwise, ∙, fontface, fontsize, @layer
import MechanicalSketch: bin_bounds, BinLegend, draw_legend, StaticArrays, ColorSchemes
import ColorSchemes:     PuOr_8, ColorScheme
import StaticArrays:     SVector

using Test

let
if !@isdefined s⁻¹
    @import_expand(m,s,°)
end

BACKCOLOR = color_with_lumin(PALETTE[8], 80)
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_20.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
end
restart()

"Source position"
p_source = (-7.5, 0.0)m |> ComplexQuantity

"Flow rate, 2d flow"
q_source = 1.0m²/s

"2d velocity potential function. Complex quantity domain, real quantity range."
ϕ_source = generate_complex_potential_source(; pos = p_source, massflowout = q_source)

physwidth = 20m
physheight = physwidth / 3
xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
A = [ϕ_source(complex(x, y)) for y in ys, x in xs] # 0.189560 seconds (1.17 M allocations: 34.138 MiB, 64.04% compilation time)

# Prepare a plot
ma = maximum(A)
mi = minimum(A)
binwidth3 = (ma - mi) / 3
legend3 = BinLegend(;maxlegend = ma, minlegend = mi, binwidth = binwidth3)
tfoo(x) = normalize_binwise(length(bin_bounds(legend3)), bin_bounds(legend3), x)

@test tfoo(mi) == 0
@test tfoo(ma) == 1
@test tfoo(mi + binwidth3) == 1 / 2
@test tfoo(mi + 2binwidth3) == 1
@test tfoo(ma) == 1
@test tfoo(mi + 0.5 * binwidth3) == 0
@test tfoo(mi + 1.5 * binwidth3) == 1 / 2
@test tfoo(mi + 2.5 * binwidth3) == 1

@test tfoo(mi + 0.6 * binwidth3) == 0.0
@test isnan(tfoo( mi + 3.001 * binwidth3))
@test isnan(tfoo(mi - 0.001 * binwidth3))
bibo = bin_bounds(legend3)
@test bibo[1] == mi
@test bibo[2] == mi + binwidth3
@test bibo[4] == ma
@test legend3(mi) != legend3.nan_color
@test legend3(ma) != legend3.nan_color
legend3.(bibo)

# Since we have values on both sides of zero, create a legend including zero.
binwidth = 0.1m²∙s⁻¹
binbounds = Tuple((-5:5)binwidth)
# Distinguishing between positive and negative is often important. We
# truncate a colorscheme to make our own.
cv(x) = get(PuOr_8, x, :clamp)
nomiddle = [[cv(x) for x in range(0, 0.3, length = 10)];
        [cv(1 - x) for x in range(0.3, 0, length = 11)]]
nomiddlescheme = ColorScheme(SVector(nomiddle...), "testschemes", "emphasize plus / minus")
legend = BinLegend(;binbounds, colorscheme = nomiddlescheme, name = :Potential)

colormat = legend.(A);
ulp, lrp = place_image(O, colormat);

# Top left of legend
legendpos = O + ((lrp[1] - ulp[1]) / 2 + EM, (ulp[2] - lrp[2]) / 2)
@layer begin
    # Font for the 'toy' text interface
    fontface("Calibri")
    fontsize(0.9FS)
    draw_legend(legendpos, legend)
end

# Add some decoration to the plot
sethue(color_with_lumin(PALETTE[5], 30))
dimension_aligned(O + (-physwidth / 2, physheight / 2), O + (physwidth / 2, physheight / 2))
dimension_aligned(O + p_source, O + (-7.5 + 1.0im)m)
dimension_aligned(O + (-physwidth / 2, - physheight / 2 ),  O +  (-physwidth / 2, physheight / 2 ))

@layer begin # Leaders
    sethue(color_with_lumin(PALETTE[5], 90))
    p0 = (-7.5 + 1.0im)m
    arrow(O + p0 + (2EM, -EM), O + p0)
    fvalue = ϕ_source(p0)
    strvalue = string(round(typeof(fvalue ), fvalue , digits = 3))
    strargument = string(p0)
    txt =  "ϕsource(" * strargument * ") = " * strvalue
    settext(txt, O + p0 + (3EM, -EM), markup=true)
    p0 = (-4.0 - 2.0im)m
    arrow(O + p0 + (2EM, -EM), O + p0)
    fvalue = ϕ_source(p0)
    strargument = string_polar_form(p0)
    strvalue = string(round(typeof(fvalue), fvalue , digits = 3))
    txt =  "ϕsource( " * strargument * " ) = " * strvalue
    settext(txt, O + p0 + (3EM, -EM), markup=true)
end

# Convert the function to latex and display
la1 = latexify(quote
    ϕ_source(r:: Length)::Potential = $q_source * log(r / oneunit(r)) / (2π)
end)
pt = O + (-14.0, -7 / 2)m
ptul, ptlr, scalefactor = place_image(pt, la1, height = 3EM)

# Following text
setfont("DejaVu Sans", FS)
str = "exists for irrotational flows only"
settext(str, ptlr + (EM, -0.5FS), markup = true)

# Convert the outer function to latex
la2 = latexify(quote
    ϕ_source(p::ComplexQuantity)::Potential = ϕ(hypot(p - $p_source))
end)
place_image(pt + (0, 5EM), la2; scalefactor)

finish()

end # let


text:
let
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 1, 4) == [2,3,4,1]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 2, 4) == [3,4,1,2]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 1, 4) == [2,3,4,1,5,6,7]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 2, 4) == [3,4,1,2,5,6,7]

    td = get_test_trajdata()
    scene = Scene()

    col = get_first_collision(get!(scene, td, 1))
    @test col.is_colliding
    @test col.A == 1
    @test col.B == 2

    @test get_first_collision(get!(scene, td, 2), CPAMemory()).is_colliding == true
end


text:
# spherical bessel functions
sbesselj(nu, x) = √(π./2x)*besselj.(nu+1/2, x)
sbessely(nu, x) = √(π./2x)*bessely.(nu+1/2, x)
shankelh1(nu,x) = sbesselj.(nu,x) + 1im.*sbessely.(nu,x)

function cart2d2sph2d(xyz,ax_horz,ax_vert)
# generate 2D meshgrid in spherical coordinates system
# corresponding to Cartesian 2D profile (xz, yz, xy)
# for calculation of field profiles
# Input:
# xyz: 1->xz, 2->yz, 3->xy
# ax_horz, ax_vert: 1D array of axes
    mat1 = ax_horz' .* ones(length(ax_vert))
    mat2 = ones(length(ax_horz))' .* ax_vert

    if xyz == 1
        # xz plot range
        mat_r,mat_theta,mat_phi = cart2sph(mat1,0,mat2)
    elseif xyz == 2
        # yz plot range
        mat_r,mat_theta,mat_phi = cart2sph(0,mat1,mat2)
    elseif xyz == 3
        # xy plot range
        mat_r,mat_theta,mat_phi = cart2sph(mat1,mat2,0)
    end
    return mat_r,mat_theta,mat_phi
end

function cart2sph(x,y,z)
# convert Cartesian coordinates(x,y,z) to spherical coordinates(r,θ,φ)
# θ is a polar angle (θ=0 is z-axis), φ is an azimuthal angle (φ=0 is x-axis)
    r = sqrt.(x.^2 .+ y.^2 .+ z.^2)
    φ = mod.(atan.(y,x) .+ 2*π, 2*π)
    θ = acos.(z./r)
    return r,θ,φ
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

board = []
@test """
""" == table(board) |> string

board = [11 12 13; 21 22 23]
@test """
+---+----+----+----+
|   |  1 |  2 |  3 |
+===+====+====+====+
| 1 | 11 | 12 | 13 |
+---+----+----+----+
| 2 | 21 | 22 | 23 |
+---+----+----+----+""" == table(board) |> string

board = ([1 2], [5 6;7 8], [9 10; 11 12])
@test """
+---+-------+
|   |     1 |
+===+=======+
| 1 |   1 2 |
+---+-------+
|   |   5 6 |
| 2 |   7 8 |
+---+-------+
|   |  9 10 |
| 3 | 11 12 |
+---+-------+""" == table(board) |> string

@test """
+---+--------+
|   | result |
+===+========+
| x |    1 2 |
+---+--------+
|   |    5 6 |
| y |    7 8 |
+---+--------+
|   |   9 10 |
| z |  11 12 |
+---+--------+""" ==
  table(board, colnames=["result"], rownames=["x", "y", "z"]) |> string

board = ([1 2], [5 6;7 8], [9 10; 11 12])
@test """
+---+-------+
|   |     1 |
+===+=======+
| 1 |   1 2 |
+---+-------+
|   |   5 6 |
| 2 |   7 8 |
+---+-------+
|   |  9 10 |
| 3 | 11 12 |
+---+-------+""" == table(board) |> string

board = ([], [5 6;7 8], [9 10; 11 12])
@test """
+---+-------+
|   |     1 |
+===+=======+
| 1 |       |
+---+-------+
|   |   5 6 |
| 2 |   7 8 |
+---+-------+
|   |  9 10 |
| 3 | 11 12 |
+---+-------+""" == table(board) |> string

board = ([1 2;3 4], [5 6;7 8], [9 10; 11 12])
@test """
+---+-------+
|   |     1 |
+===+=======+
|   |   1 2 |
| 1 |   3 4 |
+---+-------+
|   |   5 6 |
| 2 |   7 8 |
+---+-------+
|   |  9 10 |
| 3 | 11 12 |
+---+-------+""" == table(board) |> string

Base.eval(:(have_color = true))


text:
# ---
# title: 957. Prison Cells After N Days
# id: problem957
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/prison-cells-after-n-days/description/>
# hidden: true
# ---
#
# There are 8 prison cells in a row, and each cell is either occupied or vacant.
#
# Each day, whether the cell is occupied or vacant changes according to the
# following rules:
#
#   * If a cell has two adjacent neighbors that are both occupied or both vacant, then the cell becomes occupied.
#   * Otherwise, it becomes vacant.
#
# (Note that because the prison is a row, the first and the last cells in the
# row can't have two adjacent neighbors.)
#
# We describe the current state of the prison in the following way: `cells[i] ==
# 1` if the `i`-th cell is occupied, else `cells[i] == 0`.
#
# Given the initial state of the prison, return the state of the prison after
# `N` days (and `N` such changes described above.)
#
#
#
# **Example 1:**
#
#
#
#     Input: cells = [0,1,0,1,1,0,0,1], N = 7
#     Output: [0,0,1,1,0,0,0,0]
#     Explanation: The following table summarizes the state of the prison on each day:
#     Day 0: [0, 1, 0, 1, 1, 0, 0, 1]
#     Day 1: [0, 1, 1, 0, 0, 0, 0, 0]
#     Day 2: [0, 0, 0, 0, 1, 1, 1, 0]
#     Day 3: [0, 1, 1, 0, 0, 1, 0, 0]
#     Day 4: [0, 0, 0, 0, 0, 1, 0, 0]
#     Day 5: [0, 1, 1, 1, 0, 1, 0, 0]
#     Day 6: [0, 0, 1, 0, 1, 1, 0, 0]
#     Day 7: [0, 0, 1, 1, 0, 0, 0, 0]
#
#
#
# **Example 2:**
#
#
#
#     Input: cells = [1,0,0,1,0,0,1,0], N = 1000000000
#     Output: [0,0,1,1,1,1,1,0]
#
#
#
#
# **Note:**
#
#   1. `cells.length == 8`
#   2. `cells[i]` is in `{0, 1}`
#   3. `1 <= N <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
################################################################################
#
#  Ringmorphisms of orders onto finite fields
#
################################################################################

# NfOrd -> FqNmod
mutable struct NfOrdToFqNmodMor <: Map{NfOrd, FqNmodFiniteField, HeckeMap, NfOrdToFqNmodMor}
  header::MapHeader{NfOrd, FqNmodFiniteField}
  poly_of_the_field::gfp_poly
  P::NfOrdIdl
  powers::Vector{nf_elem}

  function NfOrdToFqNmodMor()
    r = new()
    r.header = MapHeader{NfOrd, FqNmodFiniteField}()
    return r
  end

  function NfOrdToFqNmodMor(O::NfOrd, F::FqNmodFiniteField, g::gfp_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    d = degree(nf(O))
    p = characteristic(F)
    a = gen(nf(O))
    tmp_gfp_poly = parent(g)()
    z.poly_of_the_field = g

    powers = Vector{nf_elem}(undef, d)

    powers[1] = a
    for i in 2:d
      powers[i] = powers[i - 1] * a
    end
    z.powers = powers

    function _image(x::NfOrdElem)
      u = F()
      gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      fmpq_poly_to_gfp_poly_raw!(tmp_gfp_poly, gg)
      ccall((:nmod_poly_rem, :libflint), Nothing,
            (Ref{gfp_poly}, Ref{gfp_poly}, Ref{gfp_poly}, Ptr{Nothing}),
            tmp_gfp_poly, tmp_gfp_poly, g, pointer_from_objref(F)+sizeof(fmpz))
      ccall((:fq_nmod_set, :libflint), Nothing,
            (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}),
            u, tmp_gfp_poly, F)
      ccall((:fq_nmod_reduce, :libflint), Nothing,
            (Ref{fq_nmod}, Ref{FqNmodFiniteField}), u, F)
      return u
    end

    # The lift is even simpler!
    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      zz.elem_in_nf = mod(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end

  function NfOrdToFqNmodMor(O::NfOrd, P::NfOrdIdl)
    z = NfOrdToFqNmodMor()
    z.P = P
    p = minimum(P)
    a, g, b = get_residue_field_data(P)
    psmall = Int(p)
    R = GF(psmall, cached = false)
    Rx, x = PolynomialRing(R, "_\$", cached = false)
    F = FqNmodFiniteField(Rx(g), Symbol("_\$"), false)
    d = degree(g)
    n = degree(O)
    imageofbasis = Vector{fq_nmod}(undef, n)
    powers = Vector{nf_elem}(undef, d)
    c = Rx()
    for i in 1:n
      ib = F()
      @assert d == ncols(b[i])
      for j in 1:d
        setcoeff!(c, j - 1, b[i][1, j])
      end
      ccall((:fq_nmod_set, :libflint), Nothing, (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), ib, c, F)
      imageofbasis[i] = ib
    end

    powers[1] = a.elem_in_nf
    for i in 2:d
      powers[i] = powers[i - 1] * a.elem_in_nf
    end
    z.powers = powers

    tempF = F()

    function _image(x::NfOrdElem)
      v = coordinates(x, copy = false)
      zz = zero(F)
      for i in 1:n
        ccall((:fq_nmod_mul_fmpz, :libflint), Nothing,
              (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}),
              tempF, imageofbasis[i], v[i], F)
        add!(zz, zz, tempF)
      end
      return zz
    end

    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      zz.elem_in_nf = mod(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end
end

function preimage(f::NfOrdToFqNmodMor, y::fq_nmod)
  O = domain(f)
  p = minimum(f.P)
  powers = f.powers
  d = length(powers)
  zz = O()
  zz.elem_in_nf = nf(O)(coeff(y, 0))
  for i in 2:d
    add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
  end
  zz.elem_in_nf = mod(zz.elem_in_nf, p)
  return zz
end

# S is the type of the order, T the type of the ideal and U the elem_type of the order, which define the quotient ring
mutable struct AbsOrdQuoMap{S, T, U} <: Map{S, AbsOrdQuoRing{S, T}, HeckeMap, AbsOrdQuoMap}
  header::MapHeader{S, AbsOrdQuoRing{S, T}}

  function AbsOrdQuoMap{S, T, U}(O::S, Q::AbsOrdQuoRing{S, T}) where {S, T, U}
    z = new()

    _image = function (x::U)
      return Q(x)
    end

    _preimage = function (x::AbsOrdQuoRingElem{S, T, U})
      return x.elem
    end

    z.header = MapHeader(O, Q, _image, _preimage)
    return z
  end
end

function AbsOrdQuoMap(O::S, Q::AbsOrdQuoRing{S, T}) where {S, T}
  U = elem_type(O)
  return AbsOrdQuoMap{S, T, U}(O, Q)
end

const NfOrdQuoMap = AbsOrdQuoMap{NfOrd, NfOrdIdl, NfOrdElem}

function Mor(O::NfOrd, F::FqNmodFiniteField, y::fq_nmod)
  return NfOrdToFqNmodMor(O, F, y)
end

function Mor(O::NfOrd, F::FqNmodFiniteField, h::gfp_poly)
  return NfOrdToFqNmodMor(O, F, h)
end

function evaluate(f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  if length(f) == 0
    return parent(r)()
  end

  l = f.length-1
  s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    s = s*r + parent(r)(coeff(f, i))
  end
  return s
end

function evaluate!(z::fq_nmod, f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  zero!(z)

  if length(f) == 0
    return z
  end

  l = f.length-1
  set!(z, parent(r)(coeff(f, l)))
  #s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    mul!(z, z, r)
    add!(z, z, parent(r)(coeff(f, i)))
    #s = s*r + parent(r)(coeff(f, i))
  end
  return z
end

function _get_coeff_raw(x::fq_nmod, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ref{fq_nmod}, Int), x, i)
  return u
end

function _get_coeff_raw(x::fq, i::Int)
  t = FlintZZ
  ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Nothing, (Ref{fmpz}, Ref{fq}, Int), t, x, i)
  return t
end

# this is super slow
# improve
function (f::NfOrdQuoMap)(I::NfOrdIdl)
  O = domain(f)
  Q = codomain(f)
  B = Q.ideal + I
  nB = norm(B)
  b = basis(B, copy = false)

  z = O()

  nQ = norm(Q.ideal)
  OnQ = ideal(O, nQ)
  range1nQ2 = fmpz(1):nQ^2

  while true
    z = rand!(z, b, range1nQ2)
    #z = sum(rand(range1nQ2) * b[i] for i in 1:degree(O))
    if norm(ideal(O, z) + OnQ) == nB
      break
    end
  end

  return Q(z)
end

function (f::NfOrdQuoMap)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

mutable struct NfOrdToFqMor <: Map{NfOrd, FqFiniteField, HeckeMap, NfOrdToFqMor}
  header::MapHeader{NfOrd, FqFiniteField}
  poly_of_the_field::gfp_fmpz_poly
  P::NfOrdIdl
  fastpath::Bool
  # Some temporary variables
  tmp_gfp_fmpz_poly::gfp_fmpz_poly
  t_fmpz_poly::fmpz_poly
  t_fmpz::fmpz
  a::nf_elem

  function NfOrdToFqMor()
    z = new()
    return z
  end

  function NfOrdToFqMor(O::NfOrd, F::FqFiniteField, g::gfp_fmpz_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    z.fastpath = true
    p = characteristic(F)
    z.tmp_gfp_fmpz_poly = parent(g)()
    z.t_fmpz_poly = fmpz_poly()
    z.t_fmpz = fmpz()

    z.a = gen(nf(O))
    z.poly_of_the_field = g

    z.header = MapHeader{NfOrd, FqFiniteField}(O, F)# _image, _preimage)

    return z
  end
end

function NfOrdToFqMor(O::NfOrd, P::NfOrdIdl)#, g::fmpz_poly, a::NfOrdElem, b::Vector{fmpz_mat})
  z = NfOrdToFqMor()
  z.fastpath = false
  z.P = P
  a, g, b = get_residue_field_data(P)
  p = minimum(P)
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "_\$", cached = false)
  F = FqFiniteField(Rx(g), Symbol("_\$"), false)
  d = degree(g)
  n = degree(O)
  imageofbasis = Vector{fq}(undef, n)
  powers = Vector{nf_elem}(undef, d)
  c = Rx()

  for i in 1:n
    ib = F()
    @assert d == ncols(b[i])
    for j in 1:d
      setcoeff!(c, j - 1, b[i][1, j])
    end
    ccall((:fq_set, :libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), ib, c, F)
    imageofbasis[i] = ib
  end

  powers[1] = a.elem_in_nf
  for i in 2:d
    powers[i] = powers[i - 1] * a.elem_in_nf
  end

  tempF = F()

  function _image(x::NfOrdElem)
    v = coordinates(x, copy = false)
    zz = zero(F)
    for i in 1:n
      ccall((:fq_mul_fmpz, :libflint), Nothing,
            (Ref{fq}, Ref{fq}, Ref{fmpz}, Ref{FqFiniteField}),
            tempF, imageofbasis[i], v[i], F)
      add!(zz, zz, tempF)
    end
    return zz
  end

  function _preimage(y::fq)
    zz = O()
    zz.elem_in_nf = nf(O)(coeff(y, 0))
    for i in 2:d
      add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
    end
    zz.elem_in_nf = mod(zz.elem_in_nf, p)
    return zz
  end

  z.header = MapHeader{NfOrd, FqFiniteField}(O, F, _image, _preimage)
  return z
end

function image(f::NfOrdToFqMor, x::NfOrdElem)
  if f.fastpath
    F = codomain(f)
    O = domain(f)
    u = F()
    gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
    fmpq_poly_to_gfp_fmpz_poly_raw!(f.tmp_gfp_fmpz_poly, gg, f.t_fmpz_poly, f.t_fmpz)
    ccall((:fmpz_mod_poly_rem, :libflint), Nothing, (Ref{gfp_fmpz_poly}, Ref{gfp_fmpz_poly}, Ref{gfp_fmpz_poly}), f.tmp_gfp_fmpz_poly, f.tmp_gfp_fmpz_poly, f.poly_of_the_field)
    ccall((:fq_set, :libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), u, f.tmp_gfp_fmpz_poly, F)
    ccall((:fq_reduce, :libflint), Nothing, (Ref{fq}, Ref{FqFiniteField}), u, F)
    return u
  else
    return f.header.image(x)
  end
end

function preimage(f::NfOrdToFqMor, x::fq)
  if f.fastpath
    O = domain(f)
    F = codomain(f)
    zz = nf(O)()

    a = f.a
    # TODO: Do something more clever here
    for i in 0:degree(F)-1
      zz = zz + coeff(x, i)*a^i
    end

    return O(zz, false)
  else
    @assert isdefined(f.header, :preimage)
    return f.header.preimage(x)
  end
end

function Mor(O::NfOrd, F::FqFiniteField, h::gfp_fmpz_poly)
  return NfOrdToFqMor(O, F, h)
end

function sub(M::Nemo.MatElem{T}, r::UnitRange{<:Integer}, c::UnitRange{<:Integer}) where {T}
  z = similar(M, length(r), length(c))
  for i in 1:length(r)
    for j in 1:length(c)
      z[i, j] = M[r[i], c[j]]
    end
  end
  return z
end

mutable struct NfToFqNmodMor <: Map{AnticNumberField, FqNmodFiniteField, HeckeMap, NfToFqNmodMor}
  header::MapHeader{AnticNumberField, FqNmodFiniteField}

  function NfToFqNmodMor()
    r = new()
    r.header = MapHeader{AnticNumberField, FqNmodFiniteField}()
    return r
  end
end

mutable struct NfToFqMor <: Map{AnticNumberField, FqFiniteField, HeckeMap, NfToFqMor}
  header::MapHeader{AnticNumberField, FqFiniteField}

  function NfToFqMor()
    r = new()
    r.header = MapHeader{AnticNumberField, FqFiniteField}()
    return r
  end
end

function extend(f::Union{NfOrdToFqNmodMor, NfOrdToFqMor}, K::AnticNumberField)
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  if f isa NfOrdToFqNmodMor
    z = NfToFqNmodMor()
  elseif f isa NfOrdToFqMor
    z = NfToFqMor()
  end

  z.header.domain = K
  z.header.codomain = f.header.codomain

  p = characteristic(z.header.codomain)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  y = f(NfOrdElem(domain(f), gen(K)))
  pia = anti_uniformizer(f.P)
  O = domain(f)
  P = f.P

  #function _image(x::nf_elem)
  #  g = parent(K.pol)(x)
  #  u = inv(z.header.codomain(denominator(g)))

  #  g = Zx(denominator(g)*g)
  #  return u*evaluate(g, y)
  #end
  function _image(x::nf_elem)
    m = denominator(x, domain(f))
    l = valuation(m, P)
    if l == 0
      return f(O(m*x))//f(O(m))
    else
      return f(O(pia^l * m * x))//f(O(pia^l * m))
    end
  end

  function _preimage(x::Union{fq, fq_nmod})
    return elem_in_nf(preimage(f, x))
  end

  z.header.image = _image
  z.header.preimage = _preimage

  return z
end

function (f::NfOrdToFqNmodMor)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

function (f::NfOrdToFqMor)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

mutable struct AbsOrdToAlgAssMor{S, T} <: Map{S, AlgAss{T}, HeckeMap, AbsOrdToAlgAssMor}
  header::MapHeader

  function AbsOrdToAlgAssMor{S, T}(O::S, A::AlgAss{T}, _image::Function, _preimage::Function) where {S <: Union{ NfAbsOrd, AlgAssAbsOrd }, T}
    z = new{S, T}()
    z.header = MapHeader(O, A, _image, _preimage)
    return z
  end
end

function AbsOrdToAlgAssMor(O::Union{ NfAbsOrd, AlgAssAbsOrd }, A::AlgAss{T}, _image, _preimage) where {T}
  return AbsOrdToAlgAssMor{typeof(O), T}(O, A, _image, _preimage)
end


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

import Pkg.TOML
import LibGit2

const STDLIBS = sort!(by = lowercase, [
    "Base64"
    "CRC32c"
    "Dates"
    "DelimitedFiles"
    "Distributed"
    "FileWatching"
    "Future"
    "InteractiveUtils"
    "LibGit2"
    "Libdl"
    "LinearAlgebra"
    "Logging"
    "Markdown"
    "Mmap"
    "Pkg"
    "Printf"
    "Profile"
    "REPL"
    "Random"
    "Sockets"
    "Serialization"
    "SHA"
    "SharedArrays"
    "SparseArrays"
    "SuiteSparse"
    "Statistics"
    "Test"
    "UUIDs"
    "Unicode"
])

function uses(repo::String, tree::String, lib::String)
    pattern = string(raw"\b(import|using)\s+((\w|\.)+\s*,\s*)*", lib, raw"\b")
    success(`git -C $repo grep -Eq $pattern $tree -- '*.jl'`)
end

function gitmeta(pkgs::Dict{String,Package})
    @assert length(STDLIBS) ≤ 64 # use 64 bits to encode usage
    fd = joinpath(@__DIR__, "sha1map.toml")
    fs = joinpath(@__DIR__, "stdlib.toml")
    d = ispath(fd) ? TOML.parsefile(fd) : Dict()
    s = Dict()
    if ispath(fs)
        s = TOML.parsefile(fs)
        get(s, "STDLIBS", nothing) == STDLIBS || empty!(s)
        mv(fs, "$fs.old", force=true)
    end
    s["STDLIBS"] = STDLIBS
    io = open(fs, "w")
    println(io, "STDLIBS = [")
    for lib in STDLIBS
        println(io, "    ", repr(lib), ",")
    end
    println(io, "]")
    println(io)
    for (pkg, p) in sort!(collect(pkgs), by=first)
        (pkg == "julia" || isempty(p.versions)) && continue
        uuid = string(p.uuid)
        @info "Package [$uuid] $pkg"
        haskey(d, uuid) || (d[uuid] = Dict())
        haskey(s, uuid) || (s[uuid] = Dict())
        updated = false
        repo_path = joinpath(homedir(), ".julia", "clones", uuid)
        repo = nothing
        for (ver, v) in sort!(collect(p.versions), by=first, rev=true)
            haskey(d[uuid], v.sha1) &&
            ((v"0.7" ∉ v.requires["julia"].versions &&
              v"1.0" ∉ v.requires["julia"].versions &&
              v"1.1" ∉ v.requires["julia"].versions) ||
            haskey(s[uuid], v.sha1)) && continue
            if repo == nothing
                repo = ispath(repo_path) ? LibGit2.GitRepo(repo_path) : begin
                    updated = true
                    @info "Cloning [$uuid] $pkg $(p.url)"
                    LibGit2.clone(p.url, repo_path, isbare=true)
                end
            end
            git_commit_hash = LibGit2.GitHash(v.sha1)
            if !updated
                try LibGit2.GitObject(repo, git_commit_hash)
                catch err
                    err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
                    @info "Updating $pkg from $(p.url)"
                    LibGit2.fetch(repo, remoteurl=p.url, refspecs=["+refs/*:refs/remotes/cache/*"])
                end
            end
            failed = false
            git_commit = try LibGit2.GitObject(repo, git_commit_hash)
            catch err
                failed = true
                err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
                @error("$pkg: git object $(v.sha1) could not be found")
            end
            failed && continue
            git_commit isa LibGit2.GitCommit || git_commit isa LibGit2.GitTag ||
                error("$pkg: git object $(v.sha1) not a commit – $(typeof(git_commit))")
            git_tree = LibGit2.peel(LibGit2.GitTree, git_commit)
            @assert git_tree isa LibGit2.GitTree
            git_tree_hash = string(LibGit2.GitHash(git_tree))
            d[uuid][v.sha1] = git_tree_hash
            # scan for stdlib dependencies
            v"0.7" in v.requires["julia"].versions ||
            v"1.0" in v.requires["julia"].versions ||
            v"1.1" in v.requires["julia"].versions || continue
            haskey(s[uuid], v.sha1) && continue
            libs = [uses(repo_path, git_tree_hash, lib) for lib in STDLIBS]
            s[uuid][v.sha1] = sum(d*2^(i-1) for (i,d) in enumerate(libs))
        end
        isempty(s[uuid]) && continue
        println(io, "[$uuid]")
        for (sha1, n) in sort!(collect(s[uuid]), by=string∘first)
            println(io, "$sha1 = $n")
        end
        println(io)
        flush(io)
    end
    open(fd, "w") do io
        TOML.print(io, d, sorted=true)
    end
    close(io)
    rm("$fs.old", force=true)
    return d, s
end


text:
"""
    runfilter(fil, b0, us, ys)

Run the particle filter `fil` starting with belief `b0` with control sequence `us` and measurement sequence `ys`. Return the resulting sequence of particle collections.

This function is intended for post-processing pre-recorded data with a particle filter.

# Arguments
- `fil`: A particle filter (e.g. constructed with `BasicParticleFilter`)
- `b0`: The initial belief. Can be any distribution over initial states.
- `us`: A sequence of control inputs for the system.
- `ys`: A sequence of measurements (observations) received after each control input.

# Return value
A sequence of particle collections. This will not include `b0` and will be the same length as `us` and `ys`. Each particle collection in this sequence is the one after being updated with the corresponding measurement in `ys`.
"""
function runfilter(f::Updater, b0, us::AbstractVector, ys::AbstractVector)
    b = initialize_belief(f, b0)
    bs = Any[]
    @assert length(ys) >= length(us)
    for i in 1:length(us)
        b = update(f, b, us[i], ys[i])
        push!(bs, b)
    end
    return bs
end


text:
module VehicleModels

using Media, DifferentialEquations, Dierckx, Atom, Plots, Parameters, Interpolations

macro def(name, definition)
  return quote
    macro $name()
      esc($(Expr(:quote,definition)))
    end
  end
end

# Three DOF Vehicle Model
include("Three_DOF/F_YR.jl")
include("Three_DOF/F_YF.jl")

# other functions to export
include("Three_DOF/FZ_RL.jl")
include("Three_DOF/FZ_RR.jl")
include("Three_DOF/Ax_max.jl")
include("Three_DOF/Ax_min.jl")

# funcitons in the VehicleModels.jl package
include("utils.jl")

export
  # Objects
  Vpara,

  # Functions
  Three_DOF,
  Linear_Spline,

  # Macros and support functions
  @F_YF,
  @F_YR,
  @FZ_RL,
  @FZ_RR,
  @Ax_min,
  @Ax_max,
  @unpack_Vpara,
  @pack_Vpara
  # MAKE SURE YOU REMOVE THE FINAL COMMA!!

#############################
# types/functions/constants #
#############################
@with_kw immutable Vpara @deftype Float64
    # define model parameters
    m          = 2.6887e+03
    Izz        = 4.1101e+03
    la         = 1.5775        # distance from CoG to front axle
    lb         = 1.7245        # distance from CoG to rear axle
    FzF0       = 1.3680e+04
    FzR0       = 1.2696e+04
    dFzx_coeff = 705.
    KZX        = 705.
    KZYR       = 1040.

    # defines polynominal for acceleration bounds
    AXC::Array{Float64,1} = [-0.000128015180401862,	0.00858618595422724,	-0.225657108071454,	3.08283259993589,	-0.000138537090018958,	0.00684702729623608,	-0.120391102052425,	-3.55886697370079]

    # vehicle Limits
    x_min    = 0.
    x_max    = 400.
    y_min    = 0.
    y_max    = 400.
    sa_min   = -30*pi/180
    sa_max   = 30*pi/180
    psi_min  = -2*pi
    psi_max  = 2*pi
    u_min    = 5.   #5.
    u_max    = 29.
    sr_min   = -5*pi/180
    sr_max   = 5*pi/180
    jx_min   = -5.
    jx_max   = 5.

    # tire parameters
    FZ0     = 35000.0;
    PCY1    = 1.5874;               #Shape factor Cfy for lateral forces
    PDY1    = 0.73957;              #Lateral friction Muy
    PDY2    = -0.075004;            #Variation of friction Muy with load
    PEY1    = 0.37562;              #Lateral curvature Efy at Fznom
    PEY2    = -0.069325;            #Variation of curvature Efy with load
    PEY3    = 0.29168;              #Zero order camber dependency of curvature Efy
    PKY1    = -10.289;              #Maximum value of stiffness Kfy/Fznom
    PKY2    = 3.3343;               #Load at which Kfy reaches maximum value
    PHY1    = 0.0056509;            #Horizontal shift Shy at Fznom
    PHY2    = -0.0020257;           #Variation of shift Shy with load
    PVY1    = 0.015216;             #Vertical shift in Svy/Fz at Fznom
    PVY2    = -0.010365;            #Variation of shift Svy/Fz with load
    PC1     = PCY1;
    PD1     = PDY1 - PDY2;
    PD2     = PDY2/FZ0;
    PE1     = PEY1 - PEY2;
    PE2     = PEY2/FZ0;
    PE3     = PEY3;
    PK1     = PKY1*FZ0;
    PK2     = 1/(PKY2*FZ0);
    PH1     = PHY1 - PHY2;
    PH2     = PHY2/FZ0;
    PV1     = PVY1 - PVY2;
    PV2     = PVY2/FZ0;

    # constrained initial states
    x0_     = 200.;
    y0_     = 0.;
    psi0_   = pi/2;
    v0_     = 0.;
    u0_     = 15.;
    sa0_    = 0.; # adding additional initial states
    sr0_    = 0.;
    ax0_    = 0.;
    jx0_    = 0.;
    r0_     = 0.;

    ##TODO: Move some of these parameters somewhere else
    # weights
    w_goal = 1.; # should be zero when vehicle is not within goal distance during prediction horizon
    w_psi  = 0.01;
    w_time = 0.05;
    w_haf  = 1.0e-5;
    w_Fz   = 0.1;  #0.5 in paper
    w_ce   = 1.;
    w_sa   = 0.1;
    w_sr   = 1.;
    w_jx   = 0.01;

    sm      = 5. # (m) distance to make sure we don't hit obstacle
    Fz_min  = 1000.;
    Fz_off  = 100.;
    a_t     = Fz_min + 3*Fz_off;  # soft tire force constraint constants
    b_t     = Fz_off;
    EP      = 0.001;

    # other parameters
    L_Rd      = 5.;       # relaxation constant for LiDAR range (m)
    L_R       = 100.;     # LiDAR range (m)
    sigma     = 1.;     # 0.05 (m)small margin, if the vehicle is within this margin, then the target is considered to be reached
end

function Three_DOF(pa::Vpara,
                   x0::Vector,
                   t::Vector,
                   SR::Vector,
                   Jx::Vector,
                   t0::Float64,
                   tf::Float64)
    @unpack_Vpara pa

    # create splines
    sp_SR=Linear_Spline(t,SR);
    sp_Jx=Linear_Spline(t,Jx);

    f = (t,x,dx) -> begin
    # states
    X	  = x[1];  # 1. X position
    Y	  = x[2];  # 2. Y position
    V   = x[3];  # 3. Lateral Speed
    R   = x[4];  # 4. Yaw Rate
    SA  = x[5];  # 5. Steering Angle
    PSI = x[6];  # 6. Yaw angle
    U   = x[7];  # 7. Longitudinal Speed
    Ax  = x[8];  # 8. Longitudinal Acceleration

    # controls
    #SR  = sp_SR(t); # using Interpolations.jl
    #Jx  = sp_Jx(t); # using Interpolations.jl
    SR  = sp_SR[t];
    Jx  = sp_Jx[t];

    # diff eqs.
    dx[1]   = U*cos(PSI) - (V + la*R)*sin(PSI);    # X position
    dx[2] 	= U*sin(PSI) + (V + la*R)*cos(PSI);    # Y position
    dx[3]   = (@F_YF() + @F_YR())/m - R*U;         # Lateral Speed
    dx[4]  	= (la*@F_YF()-lb*@F_YR())/Izz;         # Yaw Rate
    dx[5]   = SR;                                  # Steering Angle
    dx[6]  	= R;                                   # Yaw Angle
    dx[7]  	= Ax;                                  # Longitudinal Speed
    dx[8]  	= Jx;                                  # Longitudinal Acceleration
  end
  tspan = [t0,tf]
  prob = ODEProblem(f, x0)
  solve(prob::ODEProblem,tspan,alg=:RK4)
end

# this vehicle model is controlled using speed and steering angle TODO finish this!!
function Three_DOF_2(pa::Vpara,
                   x0::Vector,
                   t::Vector,
                   SA::Vector,
                   U::Vector,
                   t0::Float64,
                   tf::Float64)
error("fix x0! and finish this before using!")
    @unpack_Vpara pa

    # create splines
    sp_SA=Linear_Spline(t,SA);
    sp_U=Linear_Spline(t,U);

    f = (t,x,dx) -> begin
    # states
    X	  = x[1];  # 1. X position
    Y	  = x[2];  # 2. Y position
    V   = x[3];  # 3. Lateral Speed
    r   = x[4];  # 4. Yaw Rate
    #SA  = x[5];  # 5. Steering Angle
    PSI = x[5];  # 6. Yaw angle
    #U   = x[7];  # 7. Longitudinal Speed
    #Ax  = x[6];  # 8. Longitudinal Acceleration

    # controls
    SA  = sp_SA(t);
    U  = sp_U(t);

    # diff eqs.
    dx[1]   = U*cos(PSI) - (V + la*r)*sin(PSI);    # X position
    dx[2] 	= U*sin(PSI) + (V + la*r)*cos(PSI);    # Y position
    dx[3]   = (@F_YF() + @F_YR())/m - r*U;         # Lateral Speed
    dx[4]  	= (la*@F_YF()-lb*@F_YR())/Izz;         # Yaw Rate
    #dx[5]   = SR;                                  # Steering Angle
    dx[5]  	= r;                                   # Yaw Angle
    #dx[6]  	= Ax;                                  # Longitudinal Speed
    #dx[8]  	= Jx;                                  # Longitudinal Acceleration
  end
  tspan = [t0,tf]
  prob = ODEProblem(f, x0)
  solve(prob::ODEProblem,tspan,alg=:RK4)
end

end # module


text:
using HDF5

function test_hdf5_output_layer(backend::Backend, T, eps)
  println("-- Testing HDF5 Output Layer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  tensor_dim = abs(rand(Int)) % 4 + 2
  dims = tuple((abs(rand(Int, tensor_dim)) % 8 + 1)...)
  println("    > $dims")

  input = rand(T, dims)
  input_blob = make_blob(backend, input)

  output_fn = string(Mocha.temp_filename(), ".hdf5")
  open(output_fn, "w") do file
    # create an empty file
  end
  layer = HDF5OutputLayer(bottoms=[:input], datasets=[:foobar], filename=output_fn)
  @test_throws ErrorException setup(backend, layer, Blob[input_blob], Blob[NullBlob()])

  layer = HDF5OutputLayer(bottoms=[:input], datasets=[:foobar],
      filename=output_fn, force_overwrite=true)
  state = setup(backend, layer, Blob[input_blob], Blob[NullBlob()])

  # repeat 3 times
  forward(backend, state, Blob[input_blob])
  forward(backend, state, Blob[input_blob])
  forward(backend, state, Blob[input_blob])

  shutdown(backend, state)

  expected_output = cat(tensor_dim, input, input, input)
  got_output = h5open(output_fn, "r") do h5
    read(h5, "foobar")
  end

  @test size(expected_output) == size(got_output)
  @test eltype(expected_output) == eltype(got_output)
  @test all(abs(expected_output-got_output) .< eps)

  rm(output_fn)
end

function test_hdf5_output_layer(backend::Backend)
  test_hdf5_output_layer(backend, Float32, 1e-5)
  test_hdf5_output_layer(backend, Float64, 1e-10)
end

if test_cpu
  test_hdf5_output_layer(backend_cpu)
end
if test_gpu
  test_hdf5_output_layer(backend_gpu)
end


text:
"""
    SComplex{I,R}
"""
struct SComplex{R,I} <: Number
    SComplex(r::T, i::T) where T<:Real = new{r,i}()
end

@pure Base.real(::SComplex{R,I}) where {R,I} = R()
@pure Base.imag(::SComplex{R,I}) where {R,I} = I()
@pure Base.reim(::SComplex{R,I}) where {R,I} = (R(), I())


text:
using Pkg; Pkg.activate("/Users/mbabar/Desktop/PhD/Analysis/Gerischer/")

loc = "/Users/mbabar/Desktop/PhD/Analysis/Gerischer/ElectrochemicalKinetics.jl/src/"

#using ElectrochemicalKinetics
using MAT
using CSV
using DelimitedFiles
using Printf
using Interpolations
using Trapz

# Compute QC interp
include(loc*"dos.jl")
export DOS
using .DOS: DOSData, get_dos
export DOSData, get_dos
include(loc*"quantum_capacitance.jl")
include(loc*"kinetic_models.jl")
export MarcusHushChidseyDOS

# Testing
#mhcd = MarcusHushChidseyDOS(20, 0.2, "/Users/mbabar/Desktop/PhD/Analysis/Gerischer/ElectrochemicalKinetics.jl/data/DOSes/bernal_graphene.txt")

"""
function calculate_Vdl_interp(dos_f, Vq_min, Vq_max, C_dl)
    Vq_range = range(Vq_min, Vq_max, step=0.001)
    E_min = dos_f.ranges[1][1]
    E_max = dos_f.ranges[1][end]
    CQ = [compute_cq(E_min, E_max, V, dos_f)*1.602*10 for V in Vq_range]
    Vdl_data = []
    Vappl_data = []
    for i = 1:length(Vq_range)
        V_app_i =  Vq_range[i]*((CQ[i]/C_dl) + 1)
        V_dl_i = V_app_i - Vq_range[i]
        push!(Vdl_data, V_dl_i)
        push!(Vappl_data, V_app_i)
    end
    return Vappl_data, Vdl_data
end

function QC_integrand(E, eVq, dos_func; kT=.026)
    dos_func(E) * sech((E-eVq)/(2*kT))^2/(4*kT)
end

function compute_cq(E_min, E_max, eVq, dos_func; kT=.026)
    fcn = E->QC_integrand(E, eVq, dos_func; kT=kT)
    quadgk(fcn, E_min, E_max)[1]
end
"""

Vq_min=-0.45
Vq_max=0.45
C_dl=10.0
#Vappl_data, Vdl_data = calculate_Vdl_interp(mhcd.dos.interp_func, Vq_min, Vq_max, C_dl)
#v_interp = LinearInterpolation(Vappl_data, Vdl_data)

## Debug Integrand limits

function integrand(
    mhcd::MarcusHushChidseyDOS,
    Eo,
    V_dl,
    ox::Bool;
    kT::Real = 0.026,
    V_q = 0.0,
)
    function marcus_term(E)
        local exp_arg
        if ox
            exp_arg = -(( E .+ mhcd.λ) .^ 2) ./ (4 * mhcd.λ * kT)
        else
            exp_arg = -(( E .- mhcd.λ) .^ 2) ./ (4 * mhcd.λ * kT)
        end
        exp.(exp_arg)
    end
    fd(E) = ox ? 1 .- fermi_dirac(E; kT = kT) : fermi_dirac(E; kT = kT)
    E -> mhcd.A .* ((mhcd.dos.interp_func.(E .- V_q)).^ 1) .* marcus_term(E .- (Eo .+ V_q .+ V_dl)) .* fd(E)
end

function compute_k_cq(
    V_app,
    model::MarcusHushChidseyDOS,
    ox::Bool;
    Eo = -0.07, # E_f,red (solvent) - E_f,vac (bilayer)
    C_dl = 10.0,
    Vq_min = -0.5,
    Vq_max = 0.5,
    kT = 0.026,
    E_min = model.dos.E_min,
    E_max = model.dos.E_max,
)
    V_dl_interp = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl)
    #Vappl_data, Vdl_data = calculate_Vdl_interp(model.dos.interp_func, Vq_min, Vq_max, C_dl)
    #v_interp = LinearInterpolation(Vappl_data, Vdl_data)

    V_dl = V_dl_interp(V_app)
    #V_dl = v_interp(V_app)

    V_q = V_app - V_dl
    if V_q < 0
           E_max = E_max .+ V_q
    elseif V_q > 0
    	   E_min = E_min .+ V_q
    end
    k_rate = quadgk(integrand(model, Eo, V_dl, ox; kT = kT, V_q = V_q), E_min, E_max)[1]
    return k_rate, V_q
end

# Modify string
function chop_str(str::String)
         while str[length(str)] == '0'
               str = chop(str)
         end
         if str[length(str)] == '.'
            str = chop(str)
         end
         return str
end

## Calculate rates
# Load matlab ldos data
dos_loc = "/Users/mbabar/Google Drive/My Drive/Analysis/TBLG_Carr/kp_tblg/MATLAB_vers/Data_24kpts/"
theta_list = [0.42, 0.77, 1.15, 1.34, 1.5, 2.39, 2.6, 3, 5];

#theta = theta_list[3];
k_ox = zeros(Float64, size(theta_list));
k_red = zeros(Float64, size(theta_list));
dos_ar = zeros(Float64, size(theta_list)); # ldos area till fermi level

# Interpolate Fermi level
ef_data=[[0.22 -2.7201e-18];[0.43 -0.0017];[0.57 -0.0042];[0.8 -0.0057];[1.1 -0.0078];[1.2 -0.0084];[1.36 -0.0091];[1.55 -0.0098];[1.96 -0.0108];[2.36 -0.0113];[2.67 -0.0116];[2.8 -0.0117];[3.2 -0.0119];[3.73 -0.0121];[4.4 -0.0122];[5 -0.0123]];

ef_func = LinearInterpolation(ef_data[:,1], ef_data[:,2])
kT = 0.026 #eV
lim = 1*kT # limit of dos integration
#print(ef_func(3.4))

## Identify ldos index for AA/AB
region = "AB"
# AA ldos at (1,1) in data matrix
# AB ldos at (12,12) in data matrix
if region=="AA"
   id_x = 1
   id_y = 1
elseif region=="AB"
   id_x = 12
   id_y	= 12
end

Vapp = 0.07 #V
for i in 1:size(theta_list)[1]
    local theta, k, ef, ef_id
    theta = chop_str(string(theta_list[i]));
    file = matopen(dos_loc*"ldos-"*string(theta)*"_90x90.mat")
    data = read(file, "data");
    rscx = read(file, "rscx");
    rscy = read(file, "rscy");

    # Read energy and dos csv
    d = readdlm(dos_loc*"dos-"*string(theta)*".csv", ',', Float64);
    Elist = d[:,1];
    tdos = d[:,2];

    ldos = [Elist data[:,id_x,id_y]]
    ef = ef_func(theta_list[i])
    ef_id = findmin(abs.(ldos[:,1] .- ef))[2] # index closest to fermi level
    lim_id = findmin(abs.(ldos[:,1] .- ef .+ lim))[2] # index closest to integration limit

    dos_ar[i] = trapz(ldos[lim_id:ef_id,1], ldos[lim_id:ef_id,2]) # ldos area from ef-kT to ef
    mhcd = MarcusHushChidseyDOS(1.0, 0.82, ldos, Ef=ef)
    k_ox[i], V_q = compute_k_cq(Vapp, mhcd, true; Eo=-0.07,  Vq_min=-0.45, Vq_max=0.45)
    k_red[i], V_q = compute_k_cq(Vapp, mhcd, false; Eo=-0.07,  Vq_min=-0.45, Vq_max=0.45)
    print(i," ",theta," ",k_ox[i]," ",k_red[i]," ", V_q, "\n")
end

#print(dos_ar,"\n","\n")
print(k_ox,"\n")
print(k_red,"\n")
print(k_ox-k_red,"\n")

#print(k)

# Write rate data as .mat
#file = matopen("k_data.mat", "w")
#write(file, "k_data", k_mat)
#write(file, "rscx", rscx)
#write(file, "rscy", rscy)
#close(file)

##
#@time begin
#     mhcd = MarcusHushChidseyDOS(20.0, 0.82, dos)
#     k = compute_k(0.4, mhcd; calc_cq=true, Vq_min=-0.45, Vq_max=0.45)
#end


text:
## Code for JuliaPlanners repository ##

"""
    JuliaPlannersRepo

Planning repository provided by the JuliaPlanners organization.

- Collections: None
- Domains: Specified as a string (e.g. "blocksworld")
- Problems: Formatted as "problem-i", organized by increasing difficulty.
"""
struct JuliaPlannersRepo <: PlanningRepository end

const JULIA_PLANNERS_DIR =
    normpath(dirname(pathof(@__MODULE__)), "..", "repositories", "julia-planners")

function load_domain(repo::JuliaPlannersRepo, domain::AbstractString)
    domain_dir = joinpath(JULIA_PLANNERS_DIR, domain)
    return load_domain(joinpath(domain_dir, "domain.pddl"))
end

function load_domain(repo::JuliaPlannersRepo, collection, domain)
    error("Collections not supported by $(typeof(repo))")
end

function load_problem(repo::JuliaPlannersRepo,
                      domain::AbstractString, problem::AbstractString)
    domain_dir = joinpath(JULIA_PLANNERS_DIR, domain)
    return load_problem(joinpath(domain_dir, "$problem.pddl"))
end

function load_problem(repo::JuliaPlannersRepo, domain::AbstractString, idx::Int)
    return load_problem(repo, domain, "problem-$idx")
end

function load_problem(repo::JuliaPlannersRepo, collection, domain, problem)
    error("Collections not supported by $(typeof(repo))")
end

function list_collections(repo::JuliaPlannersRepo)
    return String[]
end

function list_domains(repo::JuliaPlannersRepo)
    return readdir(JULIA_PLANNERS_DIR)
end

function list_domains(repo::JuliaPlannersRepo, collection)
    error("Collections not supported by $(typeof(repo))")
end

function list_problems(repo::JuliaPlannersRepo, domain::AbstractString)
    # Filter out non-problem files
    fns = filter(readdir(joinpath(JULIA_PLANNERS_DIR, domain))) do fn
        match(r".*\.pddl", fn) !== nothing && fn != "domain.pddl"
    end
    # Strip .pddl suffix
    return map(fn -> fn[1:end-5], fns)
end

function list_problems(repo::JuliaPlannersRepo, collection, domain)
    error("Collections not supported by $(typeof(repo))")
end


text:
# abstracttensor.jl
#
# Abstract Tensor type
#----------------------
"""
    abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end

Abstract supertype of all tensor maps, i.e. linear maps between tensor products
of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from
an input space of type `ProductSpace{S, N₂}` to an output space of type
`ProductSpace{S, N₁}`.
"""
abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end
"""
    AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0}

Abstract supertype of all tensors, i.e. elements in the tensor product space
of type `ProductSpace{S, N}`, built from elementary spaces of type `S<:IndexSpace`.

An `AbstractTensor{S, N}` is actually a special case `AbstractTensorMap{S, N, 0}`,
i.e. a tensor map with only a non-trivial output space.
"""
const AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0}

# tensor characteristics
Base.eltype(T::Type{<:AbstractTensorMap}) = eltype(storagetype(T))
similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} =
    Core.Compiler.return_type(similar, Tuple{storagetype(TT), Type{T}})

storagetype(t::AbstractTensorMap) = storagetype(typeof(t))
similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T)
Base.eltype(t::AbstractTensorMap) = eltype(typeof(t))
spacetype(t::AbstractTensorMap) = spacetype(typeof(t))
sectortype(t::AbstractTensorMap) = sectortype(typeof(t))
field(t::AbstractTensorMap) = field(typeof(t))
numout(t::AbstractTensorMap) = numout(typeof(t))
numin(t::AbstractTensorMap) = numin(typeof(t))
numind(t::AbstractTensorMap) = numind(typeof(t))

spacetype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = S
sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S)
field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S)
numout(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₁
numin(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₂
numind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₁ + N₂

const order = numind

# tensormap implementation should provide codomain(t) and domain(t)
codomain(t::AbstractTensorMap, i) = codomain(t)[i]
domain(t::AbstractTensorMap, i) = domain(t)[i]
source(t::AbstractTensorMap) = domain(t) # categorical terminology
target(t::AbstractTensorMap) = codomain(t) # categorical terminology
space(t::AbstractTensorMap) = HomSpace(codomain(t), domain(t))
space(t::AbstractTensorMap, i::Int) = space(t)[i]
dim(t::AbstractTensorMap) = dim(space(t))

# some index manipulation utilities
codomainind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} =
    ntuple(n->n, N₁)
domainind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} =
    ntuple(n-> N₁+n, N₂)
allind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} =
    ntuple(n->n, N₁+N₂)

codomainind(t::AbstractTensorMap) = codomainind(typeof(t))
domainind(t::AbstractTensorMap) = domainind(typeof(t))
allind(t::AbstractTensorMap) = allind(typeof(t))

adjointtensorindex(t::AbstractTensorMap{<:IndexSpace, N₁, N₂}, i) where {N₁, N₂} =
    ifelse(i<=N₁, N₂+i, i-N₁)

adjointtensorindices(t::AbstractTensorMap, indices::IndexTuple) =
    map(i->adjointtensorindex(t, i), indices)

# Equality and approximality
#----------------------------
function Base.:(==)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    (codomain(t1) == codomain(t2) && domain(t1) == domain(t2)) || return false
    for c in blocksectors(t1)
        block(t1, c) == block(t2, c) || return false
    end
    return true
end
function Base.hash(t::AbstractTensorMap, h::UInt)
    h = hash(codomain(t), h)
    h = hash(domain(t), h)
    for (c, b) in blocks(t)
        h = hash(c, hash(b, h))
    end
    return h
end

function Base.isapprox(t1::AbstractTensorMap, t2::AbstractTensorMap;
                atol::Real=0, rtol::Real=Base.rtoldefault(eltype(t1), eltype(t2), atol))
    d = norm(t1 - t2)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(t1), norm(t2)))
    else
        return false
    end
end

# Conversion to Array:
#----------------------
# probably not optimized for speed, only for checking purposes
function Base.convert(::Type{Array}, t::AbstractTensorMap{S, N₁, N₂}) where {S, N₁, N₂}
    I = sectortype(t)
    if I === Trivial
        convert(Array, t[])
    else
        cod = codomain(t)
        dom = domain(t)
        local A
        for (f1, f2) in fusiontrees(t)
            F1 = convert(Array, f1)
            F2 = convert(Array, f2)
            sz1 = size(F1)
            sz2 = size(F2)
            d1 = TupleTools.front(sz1)
            d2 = TupleTools.front(sz2)
            F = reshape(reshape(F1, TupleTools.prod(d1), sz1[end])*reshape(F2, TupleTools.prod(d2), sz2[end])', (d1..., d2...))
            if !(@isdefined A)
                if eltype(F) <: Complex
                    T = complex(float(eltype(t)))
                elseif eltype(F) <: Integer
                    T = eltype(t)
                else
                    T = float(eltype(t))
                end
                A = fill(zero(T), (dims(cod)..., dims(dom)...))
            end
            Aslice = StridedView(A)[axes(cod, f1.uncoupled)..., axes(dom, f2.uncoupled)...]
            axpy!(1, StridedView(_kron(convert(Array, t[f1, f2]), F)), Aslice)
        end
        return A
    end
end


text:
function test(classifier::Array{CategoricalSpiker,1}, inputs, fs;
			  mindelay = 0.001,
			  maxdelay = 0.020,
			  weight=0.5,
			  maxduration=0.150,
			  verbose=false)
	dt = 1.0/fs
	spikers = classifier
	numcategories = length(spikers)
	numchannels = length(spikers[1].dls)

	results = zeros(Int64, length(inputs))

	for (i,input) ∈ enumerate(inputs)
		unspikedidcs = collect(1:numcategories)
		spiketimes = maxdelay .* ones(Float64, numcategories) .+ 0.1
		numcontributors = zeros(Int64, numcategories)
		contributoridcs = [zeros(Int64, numchannels) for _ ∈ 1:numcategories]
		for t ∈ 0.0:dt:maxduration
			# advance buffers and update neurons
			numcontributors[unspikedidcs] .*= 0
			for spikeridx ∈ unspikedidcs
				spiker = spikers[spikeridx]
				# Advance buffers
				for (k,dl) ∈ enumerate(spiker.dls)
					inputspike = advance(dl)
					if inputspike != 0
						numcontributors[spikeridx] += 1
						contributoridcs[spikeridx][numcontributors[spikeridx]] = k
					end
				end
				# Update neuron (if any inputs)
				spiked = false
				if numcontributors[spikeridx] > 0
					spiked = srmneuronupdate(spiker.neuron, t,
										spiker.weight*numcontributors[spikeridx])
				end
				# Remove neuron from processing if fired
				if spiked
					spiketimes[spikeridx] = t
					filter!(x -> x != spikeridx, unspikedidcs)
				end
			end
			# push new input into buffers
			inputspikes = filter(s -> t <= s.t < t+dt, input)
			for spikeridx ∈ unspikedidcs
				spiker = spikers[spikeridx]
				for spike ∈ inputspikes
					pushevent(spiker.dls[spike.n])
				end
			end
		end
		firingorder = sortperm(spiketimes)
		results[i] = firingorder[1]

		# Reset classifier state
		for spiker ∈ spikers
			resetsrmneuron(spiker.neuron)
			for dl ∈ spiker.dls
				clearbuffer(dl)
			end
		end
	end

	return results
end

function trainclassifier(inputs, labels, numcategories, numchannels, fs;
						 weight=0.5,
						 μ = 0.005,
						 mindelay = 0.001,
						 maxdelay = 0.020,
						 maxduration=0.150,
						 numtrainingcycles=1,
						 τ_membrane=3.0/1000.0,
						 verbose=false)
	dt = 1.0/fs
	channels = 1:numchannels
	categories = 1:numcategories
	d_min = Int(round(mindelay*fs))
	d_max = Int(round(maxdelay*fs))
	possibledelays = d_min:d_max
	spikers = [CategoricalSpiker(1:numchannels, possibledelays, weight;
								   	maxevents = Int(round(d_max)),
									τ_membrane = τ_membrane)
			   for _ ∈ 1:numcategories]

	for i ∈ 1:numtrainingcycles
	println("Training cycle $i of $numtrainingcycles")
	numcorrect = 0
	numties = 0
	numwrong = 0
	numtested = 0
	for (label,input) ∈ zip(labels,inputs)
		unspikedidcs = collect(1:numcategories)
		spiketimes = maxdelay .* ones(Float64, numcategories) .+ 0.1
		numcontributors = zeros(Int64, numcategories)
		contributoridcs = [zeros(Int64, numchannels) for _ ∈ 1:numcategories]
		for t ∈ 0.0:dt:maxduration
			# advance buffers and update neurons
			numcontributors[unspikedidcs] .*= 0
			for spikeridx ∈ unspikedidcs
				spiker = spikers[spikeridx]
				for (k,dl) ∈ enumerate(spiker.dls)
					inputspike = advance(dl)
					if inputspike != 0
						numcontributors[spikeridx] += 1
						contributoridcs[spikeridx][numcontributors[spikeridx]] = k
					end
				end
				spiked = false
				if numcontributors[spikeridx] > 0
					spiked = srmneuronupdate(spiker.neuron, t,
										spiker.weight*numcontributors[spikeridx])
				end
				if spiked
					verbose && println("$spikeridx spiked at $t (correct: $label)")
					spiketimes[spikeridx] = t
					filter!(x -> x != spikeridx, unspikedidcs)
				end
			end

			# push events
			inputspikes = filter(s -> t <= s.t < t+dt, input)
			for spikeridx ∈ unspikedidcs
				spiker = spikers[spikeridx]
				for spike ∈ inputspikes
					pushevent(spiker.dls[spike.n])
				end
			end
		end
		verbose && println("----------")

		# Stochastic update of delay times
		firingorder = sortperm(spiketimes)
		if spiketimes[firingorder[1]] < spiketimes[firingorder[2]] && firingorder[1] == label
			if !(spiketimes[firingorder[2]] -
				 spiketimes[label] > μ)
				# Shorten random delay on correct-label spiker
				if numcontributors[label] > 0
					#println("Yay, contributors!")
					didx = rand(contributoridcs[label][1:numcontributors[label]])
				else
					didx = rand(1:numchannels)
				end
				if spikers[label].dls[didx].delaylen > d_min
					spikers[label].dls[didx].delaylen -= 1
				end
				# Lengthen random delay on random incorrect-label spiker
				while (sidx = rand(categories)) == label end
				if numcontributors[sidx] > 0
					#println("Yay, contributors!")
					didx = rand(contributoridcs[sidx][1:numcontributors[sidx]])
				else
					didx = rand(1:numchannels)
				end
				if spikers[sidx].dls[didx].delaylen < d_max
					spikers[sidx].dls[didx].delaylen += 1
				end
			end
			numcorrect += 1
		else
			# Shorten random delay on correct-label spiker
			if numcontributors[label] > 0
				#println("Yay, contributors!")
				didx = rand(contributoridcs[label][1:numcontributors[label]])
			else
				didx = rand(channels)
			end
			if spikers[label].dls[didx].delaylen > d_min
				spikers[label].dls[didx].delaylen -= 1
			end
			# Lengthen delay on random incorrect-label spiker that spikes before
			# desired category
			#prespikeridces = filter(x -> x <= label, firingorder)
			#prespikeridces = firingorder[1:label-1]
			prespikeridcs = findall(x -> x <= spiketimes[label],
								   spiketimes)
			filter!(x -> x != label, prespikeridcs)
			sidx = rand(prespikeridcs)
			if numcontributors[sidx] > 0
				#println("Yay, contributors!")
				didx = rand(contributoridcs[sidx][1:numcontributors[sidx]])
			else
				didx = rand(1:numchannels)
			end
			if spikers[sidx].dls[didx].delaylen < d_max
				spikers[sidx].dls[didx].delaylen += 1
			end
			#spikers[idx].dls[rand(channels)].delaylen += 1
			if spiketimes[firingorder[1]] == spiketimes[firingorder[2]]
				numties += 1
			else
				numwrong += 1
			end
		end
		numtested += 1
		# Reset neurons between inputs
		for spiker ∈ spikers
			resetsrmneuron(spiker.neuron)
			for dl ∈ spiker.dls
				clearbuffer(dl)
			end
		end
	end
	println("Correct: $(numcorrect/numtested)")
	println("Wrong: $(numwrong/numtested)")
	println("Ties: $(numties/numtested)")
	println()
	end

	return spikers
end

function trainmulticlassifier(inputs, labels, numcategories,
                              numchannels, numneurons, fs;
                              weight=0.5,
                              μ = 0.005,
                              mindelay = 0.001,
                              maxdelay = 0.020,
                              maxduration=0.150,
                              numtrainingcycles=1,
                              τ_membrane=3.0/1000.0,
                              verbose=false)
	dt = 1.0/fs
	channels = 1:numchannels
	categories = 1:numcategories
	d_min = Int(round(mindelay*fs))
	d_max = Int(round(maxdelay*fs))
	possibledelays = d_min:d_max
	# spikers = [[CategoricalSpiker(1:numchannels, possibledelays, weight;
	# 							   	maxevents = Int(round(d_max)),
    #                                 τ_membrane = τ_membrane) for _ ∈ 1:numneurons]
    #             for _ ∈ 1:numcategories]
	# spikers = [[CategoricalSpiker(1:numchannels, possibledelays, weight;
	# 							   	maxevents = Int(round(d_max)),
    #                                 τ_membrane = τ_membrane) for _ ∈ 1:numneurons]
    #             for _ ∈ 1:numcategories]
    spikers = Array{CategoricalSpiker,2}(undef,(numcategories,numneurons))
    for r ∈ 1:numcategories
        for c ∈ 1:numneurons
            spikers[r,c] = CategoricalSpiker(1:numchannels, possibledelays, weight;
                                             maxevents = Int(round(d_max)),
                                             τ_membrane = τ_membrane)
        end
    end

	for i ∈ 1:numtrainingcycles
        println("Training cycle $i of $numtrainingcycles")
        numcorrect = 0
        numties = 0
        numwrong = 0
        numtested = 0
        for (label,input) ∈ zip(labels,inputs)
            unspikedcats = [collect(1:numcategories) for _ ∈ 1:numneurons]
            spiketimes = ones(Float64, numcategories, numneurons)
            numcontributors = zeros(Int64, numcategories, numneurons)
		    # contributoridcs = [zeros(Int64, numchannels) for _ ∈ 1:numcategories]
            contributoridcs = zeros(Int64, numchannels, numcategories, numneurons)

            # Run input through classifier neurons
            for t ∈ 0.0:dt:maxduration
                for n ∈ 1:numneurons
                    # advance buffers and update neurons
                    numcontributors[unspikedcats[n],n] .= 0     # <--- here's the problem
                    for spikercat ∈ unspikedcats[n]
                        spiker = spikers[spikercat,n]
                        for (k,dl) ∈ enumerate(spiker.dls)
                            inputspike = advance(dl)
                            if inputspike != 0
                                numcontributors[spikercat,n] += 1
						        contributoridcs[numcontributors[spikercat,n], spikercat, n] = k
                            end
                        end
                        spiked = false
                        if numcontributors[spikercat,n] > 0
                            spiked = srmneuronupdate(spiker.neuron, t,
                                                     spiker.weight*numcontributors[spikercat,n])
                        end
                        if spiked
                            verbose && println("$spikercat ($(n)) spiked at $t (correct: $label)")
                            spiketimes[spikercat,n] = t
                            filter!(x -> x != spikercat, unspikedcats[n])
                        end
                    end
                    # push events
                    inputspikes = filter(s -> t <= s.t < t+dt, input)
                    for spikercat ∈ unspikedcats[n]
                        spiker = spikers[spikercat,n]
                        for spike ∈ inputspikes
                            pushevent(spiker.dls[spike.n])
                        end
                    end
                end
            end
            verbose && println("----------")

		    # Stochastic update of delay times
            incatorders = zeros(Int64, size(spiketimes))
            for k ∈ 1:numcategories
                incatorders[k,:] = sortperm(spiketimes[k,:])
            end
            besttimes = [spiketimes[c,incatorders[c,1]] for c ∈ 1:numcategories]
            crosscatorders = sortperm(besttimes)
            firsttospike = findfirst(i->i==1, crosscatorders)
            secondtospike = findfirst(i->i==2, crosscatorders)

            if firsttospike != label
                # If wrong category spiked first, shorten the correct-category
                # spiker that spiked with the shortest latency with the shortest
                # latency & lengthen the wrong-category spiker with the shortest
                # latency
                numwrong += 1

                # shorten shortest latency right-category spiker
                if numcontributors[label,incatorders[label,1]] > 0
                    didx = rand(contributoridcs[1:numcontributors[label,incatorders[label,1]],label,incatorders[label,1]])
 				else
 					didx = rand(1:numchannels)
 				end
                if spikers[label,incatorders[label,1]].dls[didx].delaylen > d_min
                    spikers[label,incatorders[label,1]].dls[didx].delaylen -= 1
 				end

                # lengthen shortest latency wrong-category spiker
                if numcontributors[firsttospike,incatorders[firsttospike,1]] > 0
                    didx = rand(contributoridcs[1:numcontributors[firsttospike,incatorders[firsttospike,1]],firsttospike,incatorders[firsttospike,1]])
                else
                    didx = rand(1:numchannels)
                end
                if spikers[firsttospike,incatorders[firsttospike,1]].dls[didx].delaylen < d_max
                    spikers[firsttospike,incatorders[firsttospike,1]].dls[didx].delaylen += 1
                end
            else
                # Otherwise, check if margin is large enough. If not, do same as
                # above, otherwise do nothing.
                if besttimes[1] == besttimes[2]
                    numties += 1
                else
                    numcorrect += 1
                end
                if besttimes[2] - besttimes[1] < μ
                    if numcontributors[label,incatorders[label,1]] > 0
                        didx = rand(contributoridcs[1:numcontributors[label,incatorders[label,1]],label,incatorders[label,1]])
                    else
                        didx = rand(1:numchannels)
                    end
                    if spikers[label,incatorders[label,1]].dls[didx].delaylen > d_min
                        spikers[label,incatorders[label,1]].dls[didx].delaylen -= 1
                    end

                    # lengthen shortest latency wrong-category spiker
                    if numcontributors[firsttospike,incatorders[firsttospike,1]] > 0
                        didx = rand(contributoridcs[1:numcontributors[firsttospike,incatorders[firsttospike,1]],firsttospike,incatorders[firsttospike,1]])
                    else
                        didx = rand(1:numchannels)
                    end
                    if spikers[firsttospike,incatorders[firsttospike,1]].dls[didx].delaylen < d_max
                        spikers[firsttospike,incatorders[firsttospike,1]].dls[didx].delaylen += 1
                    end
                end
            end

     		# Reset neurons between inputs
            for spiker ∈ spikers
                resetsrmneuron(spiker.neuron)
                for dl ∈ spiker.dls
                    clearbuffer(dl)
                end
            end
     		numtested += 1
        end
        correctrate = numcorrect/numtested
        println("Correct: $(correctrate)")
        println("Wrong: $(numwrong/numtested)")
        println("Ties: $(numties/numtested)")
        println()
        if correctrate == 1.0
            break
        end
    end

    return spikers
end

function testmulti(spikers::Array{CategoricalSpiker,2}, inputs, fs;
                              weight=0.5,
                              maxduration=0.150,
                              verbose=false)
	dt = 1.0/fs
    numcategories, numneurons = size(spikers)
	numchannels = length(spikers[1,1].dls)

	results = zeros(Int64, length(inputs))

    for (i,input) ∈ enumerate(inputs)
        unspikedcats = [collect(1:numcategories) for _ ∈ 1:numneurons]
        spiketimes = ones(Float64, numcategories, numneurons)
        numcontributors = zeros(Int64, numcategories, numneurons)
        # contributoridcs = [zeros(Int64, numchannels) for _ ∈ 1:numcategories]
        contributoridcs = zeros(Int64, numchannels, numcategories, numneurons)

        # Run input through classifier neurons
        for t ∈ 0.0:dt:maxduration
            for n ∈ 1:numneurons
                # advance buffers and update neurons
                numcontributors[unspikedcats[n],n] .= 0
                for spikercat ∈ unspikedcats[n]
                    spiker = spikers[spikercat,n]
                    for (k,dl) ∈ enumerate(spiker.dls)
                        inputspike = advance(dl)
                        if inputspike != 0
                            numcontributors[spikercat,n] += 1
                            contributoridcs[numcontributors[spikercat,n], spikercat, n] = k
                        end
                    end
                    spiked = false
                    if numcontributors[spikercat,n] > 0
                        spiked = srmneuronupdate(spiker.neuron, t,
                                                 spiker.weight*numcontributors[spikercat,n])
                    end
                    if spiked
                        verbose && println("$spikercat ($(n)) spiked at $t (correct: $label)")
                        spiketimes[spikercat,n] = t
                        filter!(x -> x != spikercat, unspikedcats[n])
                    end
                end
                # push events
                inputspikes = filter(s -> t <= s.t < t+dt, input)
                for spikercat ∈ unspikedcats[n]
                    spiker = spikers[spikercat,n]
                    for spike ∈ inputspikes
                        pushevent(spiker.dls[spike.n])
                    end
                end
            end
        end

        # Find first firerer
        incatorders = zeros(Int64, size(spiketimes))
        for k ∈ 1:numcategories
            incatorders[k,:] = sortperm(spiketimes[k,:])
        end
        besttimes = [spiketimes[c,incatorders[c,1]] for c ∈ 1:numcategories]
        crosscatorders = sortperm(besttimes)
        results[i] = findfirst(x->x==1, crosscatorders)

        # Reset classifiers
        for spiker ∈ spikers
            resetsrmneuron(spiker.neuron)
            for dl ∈ spiker.dls
                clearbuffer(dl)
            end
        end
    end

    return results
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
module AdvancedHMC

const DEBUG = convert(Bool, parse(Int, get(ENV, "DEBUG_AHMC", "0")))

using Statistics: mean, var, middle
using LinearAlgebra: Symmetric, UpperTriangular, mul!, ldiv!, dot, I, diag, cholesky, UniformScaling
using StatsFuns: logaddexp, logsumexp
using Random: GLOBAL_RNG, AbstractRNG
using ProgressMeter: ProgressMeter
using Parameters: @unpack, reconstruct
using ArgCheck: @argcheck

import StatsBase: sample
import Parameters: reconstruct

include("utilities.jl")

# Notations
# ℓπ: log density of the target distribution
# θ: position variables / model parameters
# ∂ℓπ∂θ: gradient of the log density of the target distribution w.r.t θ
# r: momentum variables
# z: phase point / a pair of θ and r

include("metric.jl")
export UnitEuclideanMetric, DiagEuclideanMetric, DenseEuclideanMetric

include("hamiltonian.jl")
export Hamiltonian

include("integrator.jl")
export Leapfrog, JitteredLeapfrog, TemperedLeapfrog

include("trajectory.jl")
@deprecate find_good_eps find_good_stepsize
export EndPointTS, SliceTS, MultinomialTS,
       StaticTrajectory, HMCDA, NUTS,
       ClassicNoUTurn, GeneralisedNoUTurn,
       find_good_stepsize

include("adaptation/Adaptation.jl")
using .Adaptation
import .Adaptation: StepSizeAdaptor, MassMatrixAdaptor, StanHMCAdaptor, NesterovDualAveraging

# Helpers for initializing adaptors via AHMC structs

StepSizeAdaptor(δ::AbstractFloat, stepsize::AbstractScalarOrVec{<:AbstractFloat}) =
    NesterovDualAveraging(δ, stepsize)
StepSizeAdaptor(δ::AbstractFloat, i::AbstractIntegrator) = StepSizeAdaptor(δ, nom_step_size(i))

MassMatrixAdaptor(m::UnitEuclideanMetric{T}) where {T} =
    UnitMassMatrix{T}()
MassMatrixAdaptor(m::DiagEuclideanMetric{T}) where {T} =
    WelfordVar{T}(size(m); var=copy(m.M⁻¹))
MassMatrixAdaptor(m::DenseEuclideanMetric{T}) where {T} =
    WelfordCov{T}(size(m); cov=copy(m.M⁻¹))

MassMatrixAdaptor(
    m::Type{TM},
    sz::Tuple{Vararg{Int}}=(2,)
) where {TM<:AbstractMetric} = MassMatrixAdaptor(Float64, m, sz)

MassMatrixAdaptor(
    ::Type{T},
    ::Type{TM},
    sz::Tuple{Vararg{Int}}=(2,)
) where {T, TM<:AbstractMetric} = MassMatrixAdaptor(TM(T, sz))

# Deprecations

@deprecate StanHMCAdaptor(n_adapts, pc, ssa) initialize!(StanHMCAdaptor(pc, ssa), n_adapts)
@deprecate NesterovDualAveraging(δ::AbstractFloat, i::AbstractIntegrator) StepSizeAdaptor(δ, i)
@deprecate Preconditioner(args...) MassMatrixAdaptor(args...)

export StepSizeAdaptor, NesterovDualAveraging,
       MassMatrixAdaptor, UnitMassMatrix, WelfordVar, WelfordCov,
       NaiveHMCAdaptor, StanHMCAdaptor

include("diagnosis.jl")

include("sampler.jl")
export sample

### Init

using Requires

function __init__()
    include(joinpath(@__DIR__, "contrib", "diffeq.jl"))
    include(joinpath(@__DIR__, "contrib", "ad.jl"))
end

end # module