text:
"""

```
VMDselect(inputfile::String, selection::String; vmd="vmd" )
```

Select atoms using vmd selection syntax, with vmd in background

Returns the list of index (one-based) and atom names

Function to return the selection from a input file (topology, coordinates, etc),
by calling VMD in the background.

"""
function VMDselect(inputfile::String, selection::String; vmd = "vmd")

    if !isfile(inputfile)
        error("Could not find file: $inputfile")
    end

    vmd_input = Base.open("./VMDINPUT_TMP.VMD", "w")
    Base.write(vmd_input, "mol new \"$inputfile\" \n")
    Base.write(vmd_input, "set sel [ atomselect top \"$selection\" ] \n")
    Base.write(vmd_input, "puts \"INDEXLIST\" \n")
    Base.write(vmd_input, "set indexes [ \$sel get index ] \n")
    Base.write(vmd_input, "puts \"ENDINDEXLIST\" \n")
    Base.write(vmd_input, "puts \"NAMELIST\" \n")
    Base.write(vmd_input, "set names [ \$sel get name ] \n")
    Base.write(vmd_input, "puts \"ENDNAMELIST\" \n")
    Base.write(vmd_input, "exit \n")
    Base.close(vmd_input)

    vmd_output = Base.read(`$vmd -dispdev text -e ./VMDINPUT_TMP.VMD`, String)

    # Read indexes
    local index_list::String
    readnext = false
    for line in split(vmd_output, "\n")
        if readnext
            if line == "ENDINDEXLIST"
                error("ERROR: Selection '$selection' does not contain any atom")
            end
            index_list = line
            break
        end
        if line == "INDEXLIST"
            readnext = true
        end
    end
    index_split = split(index_list)
    nsel = length(index_split)
    selection_indexes = Vector{Int}(undef, nsel)
    for i = 1:nsel
        selection_indexes[i] = parse(Int, index_split[i]) + 1
    end

    # Read atom names
    local name_list::String
    readnext = false
    for line in split(vmd_output, "\n")
        if readnext
            if line == "ENDNAMELIST"
                error("ERROR: Selection '$selection' does not contain any atom")
            end
            name_list = line
            break
        end
        if line == "NAMELIST"
            readnext = true
        end
    end
    name_split = split(name_list)
    nsel = length(name_split)
    selection_names = Vector{String}(undef, nsel)
    for i = 1:nsel
        selection_names[i] = strip(name_split[i])
    end
    run(`\rm -f ./VMDINPUT_TMP.VMD`)

    return selection_indexes, selection_names

end


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

using Markdown
using InteractiveUtils

# ╔═╡ ccccd3d4-5b5c-11eb-2d48-8593366fb942
using PhaseSpaceTools, Plots, LaTeXStrings

# ╔═╡ b45b5d84-5b5c-11eb-2189-371f4b3bf28a
html"<button onclick=present()>Present</button>"

# ╔═╡ 9535e4b6-9308-4077-a6c6-f5860a8da4c2
md"""
# Sampling quantum states
Sampling single mode quantum states for use in quantum phase space simulations using the truncated Wigner and +P distributions. Similar methods can be applied to multimode simulations of Bose-Einstein condensates.

We use the methods presented in the article
- [Olsen, Bradley, Opt. Comm. 282 (2009) 3924-3929](https://doi.org/10.1016/j.optcom.2009.06.033)

with erratum [Olsen, Lewis-Swan, Bradley, Opt. Comm. 370 (2016) 327-328](https://doi.org/10.1016/j.optcom.2016.02.068)

A complete erratum can be found at [SciPost](https://scipost.org/commentaries/10.1016/j.optcom.2009.06.033/).

# New States
We have add sampling for the additional distributions:
- Fock state for +W
- Thermal state for +P
- Bogoliubov states for P,Q,W

# Papers
For more information on Bose-Einstein condensate applications, see the review

_Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques_, P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, [Advances in Phyiscs 57, 363 (2008)](http://dx.doi.org/10.1080/00018730802564254)

For a recent paper using this package, see

_Steady states, squeezing, and entanglement in intracavity triplet down conversion_,
Mathew D. E. Denys, Murray K. Olsen, Luke S. Trainor, Harald G. L. Schwefel, and Ashton S. Bradley,  [Optics Communications 484 (2021) 126699](https://doi.org/10.1016/j.optcom.2020.126699), [arXiv](https://arxiv.org/abs/1907.09572)
"""

# ╔═╡ 2a3a9c5e-5b67-11eb-2250-41d162813a4a
md"# Examples"

# ╔═╡ bf17f9ca-8dcf-483d-a399-42bdb69f79f9
gr(size=(350,350),xlabel=L"\alpha_r",ylabel=L"\alpha_i",
	    xlims=(-20,20),ylims=(-20,20),ms=.2,
	    legend=false,grid=false,aspect_ratio=1,colorbar=true)

# ╔═╡ adce375f-31dd-4bbd-896a-6bdf823a4456
function show_meansW(N,State,p,Dist)
	s = State(p...)
	a,ā = Dist(s,N)
    n̄ = mean(a.*ā)-.5
    Vn= mean(a.^2 .*ā.^2)-mean(a.*ā)-n̄.^2
    Text("
	N = $N samples.\n
    Averages:\n
	<â> = $(mean(a))\n
	<â⁺â> = $(n̄)\n
	V(n̂) = $(Vn)\n
	Relative width = $(sqrt(abs(Vn))/abs(n̄))")
end

# ╔═╡ eba146d8-5b5e-11eb-363c-1d84d1f72610
function sample_plot(N,State,p,Dist)
	s = State(p...)
	a,ā = Dist(s,N)
	scatter(real(a),imag(a),c=:red)
end

# ╔═╡ 86a58536-5b66-11eb-2e94-0f42924c2bee
md"# Coherent State
Let's look at the coherent state in different representations."

# ╔═╡ 2421fc20-5b5f-11eb-35c0-253c6b5420f6
sample_plot(500,Coherent,10,wigner)

# ╔═╡ e5fedb01-83be-4fe0-9072-be338998707a
show_meansW(100000,Coherent,10,wigner) # Poisson number statistics

# ╔═╡ 7b1a6556-5b66-11eb-21c0-53b40b983fab
md"# Thermal State"

# ╔═╡ 9d51cea0-5b5e-11eb-17f8-b9ec85cfd39a
sample_plot(500,Thermal,(0,10),wigner)

# ╔═╡ 130e58fa-72b1-4005-8e3b-c624da6db44f
show_meansW(100000,Thermal,(0,10),wigner)

# ╔═╡ 3a3f2a44-5b66-11eb-02df-ad97df0d48d0
md"# Squeezed state"

# ╔═╡ 10c936c9-452a-4c05-a5d7-6e2d1735d3c2
begin
	β = 10
	ϕ = 0.
	r = 1
	ϵ = r*exp(2*im*ϕ)
	sample_plot(500,Squeezed,(β,ϵ),wigner)
end

# ╔═╡ b163aaf8-47c2-4d87-a92f-8703eb9f2d6d
show_meansW(100000,Squeezed,(β,ϵ),wigner)

# ╔═╡ 2c9b2906-5b66-11eb-2eb8-932c97b5e2b0
md"# Fock state"

# ╔═╡ 472762c2-5b60-11eb-361d-37601eaed6ff
sample_plot(500,Fock,100,wigner)

# ╔═╡ 37f1a96e-5b62-11eb-217f-9dc0fae96dd6
show_meansW(100000,Fock,100,wigner)

# ╔═╡ 4c12caae-5b62-11eb-1ad2-5fc5d8120fba
sample_plot(1000,Fock,100,positiveP)

# ╔═╡ c4a573b8-5b62-11eb-0726-3d42d1187b64
sample_plot(1000,Fock,100,positiveW)

# ╔═╡ e0583d34-5b62-11eb-10a8-f380c764daa7
sample_plot(1000,Fock,321,positiveP) #uses asymptotic expansion for n>320

# ╔═╡ 2df01300-5b63-11eb-0faf-9d8fb6d796c9
md"# Crescent state
A Kerr interaction will create shear for a coherent state. This can be modelled using a crescent state."

# ╔═╡ 3c8a60b4-5b63-11eb-13ef-037243089ec1
sample_plot(500,Crescent,(10,0,.5),wigner)

# ╔═╡ 6d01bc24-5b63-11eb-2fea-a3e993675754
sample_plot(500,Crescent,(10,0,.5),husimiQ)

# ╔═╡ 7befd068-5b63-11eb-15d8-a19df4d9da1b
sample_plot(500,Crescent,(10,0,.5),positiveP)

# ╔═╡ Cell order:
# ╟─b45b5d84-5b5c-11eb-2189-371f4b3bf28a
# ╟─9535e4b6-9308-4077-a6c6-f5860a8da4c2
# ╟─2a3a9c5e-5b67-11eb-2250-41d162813a4a
# ╠═ccccd3d4-5b5c-11eb-2d48-8593366fb942
# ╠═bf17f9ca-8dcf-483d-a399-42bdb69f79f9
# ╠═adce375f-31dd-4bbd-896a-6bdf823a4456
# ╠═eba146d8-5b5e-11eb-363c-1d84d1f72610
# ╟─86a58536-5b66-11eb-2e94-0f42924c2bee
# ╠═2421fc20-5b5f-11eb-35c0-253c6b5420f6
# ╠═e5fedb01-83be-4fe0-9072-be338998707a
# ╟─7b1a6556-5b66-11eb-21c0-53b40b983fab
# ╠═9d51cea0-5b5e-11eb-17f8-b9ec85cfd39a
# ╠═130e58fa-72b1-4005-8e3b-c624da6db44f
# ╟─3a3f2a44-5b66-11eb-02df-ad97df0d48d0
# ╠═10c936c9-452a-4c05-a5d7-6e2d1735d3c2
# ╠═b163aaf8-47c2-4d87-a92f-8703eb9f2d6d
# ╟─2c9b2906-5b66-11eb-2eb8-932c97b5e2b0
# ╠═472762c2-5b60-11eb-361d-37601eaed6ff
# ╠═37f1a96e-5b62-11eb-217f-9dc0fae96dd6
# ╠═4c12caae-5b62-11eb-1ad2-5fc5d8120fba
# ╠═c4a573b8-5b62-11eb-0726-3d42d1187b64
# ╠═e0583d34-5b62-11eb-10a8-f380c764daa7
# ╟─2df01300-5b63-11eb-0faf-9d8fb6d796c9
# ╠═3c8a60b4-5b63-11eb-13ef-037243089ec1
# ╠═6d01bc24-5b63-11eb-2fea-a3e993675754
# ╠═7befd068-5b63-11eb-15d8-a19df4d9da1b

ext: jl


text:
# # Free-vibration solution for a clamped plate in fluid

# ## Description

# Free-vibration solution for a clamped plate fully submerged in inviscid
# fluid. The plate is modeled with energy-stabilized NICE tetrahedra.

# For the dry plate, these are the natural frequencies:

# | Mode | Dry frequencies [rad/s] |
# | ----- | ---------------------- |
# | 1 | 12.45   |
# | 2 | 29.44   |
# | 3 | 75.05   |
# | 4 | 94.27   |
# | 5 | 106.79 |

# For the immersed (wet) plate, these are the natural frequencies:

# | Mode | Wet frequencies [rad/s] |
# | ----- | ---------------------- |
# | 1 | 7.35   |
# | 2 | 20.2   |
# | 3 | 50.45   |
# | 4 | 70.41   |
# | 5 | 78.85 |

# ## References

# [1] Fu, Y., and Price, W. G., 1987, “Interactions Between a Partially
# or Totally Immersed Vibrating Cantilever Plate and the Surrounding Fluid,”
# J. Sound Vib., 118(3), pp. 495–513.

# ## Goals

# - Define model data for an algorithm to solve the free vibration problem.

##
# ## Definitions

# Bring in needed basic facilities.
using LinearAlgebra
using Arpack

# Bring in the finite element tools, and the linear deformation support.
using FinEtools
using FinEtoolsDeforLinear
# The dry-vibration solver is provided by this algorithm module.
using FinEtoolsDeforLinear.AlgoDeforLinearModule

# The wet-vibration solver is provided by this algorithm module.
using FinEtoolsVibInFluids.AlgoVibInFluidsModule

# The basic data corresponds to a steel plate.
E = 206000*phun("MPa");# Young's modulus
nu = 0.3;# Poisson ratio
rho = 7850*phun("KG*M^-3");# mass density
# The fluid properties correspond to water.
rhow = 1000*phun("KG*M^-3");
# Dimensions of the plate: it is quite massive.
Length= 10e3*phun("MM"); Width= 10e3*phun("MM"); Height= 0.238e3*phun("MM");
# The mesh is quite coarse, but still provides engineering accuracy for the
# angular frequencies.
nHeight = 4; nLength  = 2*10; nWidth = nLength;
tolerance = Height/nHeight/100;

# The reference angular frequencies in radians per second. The dry frequencies:
dreffs = [12.45153,   29.44122,   75.04792,   94.26498,  106.7881]*phun("rad/s")
# The reference angular frequencies in radians per second. The wet frequencies:
wreffs = [7.35, 20.2, 50.45, 70.41, 78.85]*phun("rad/s")
# Solve for this many natural frequencies.
neigvs = length(wreffs);

# Generate a tetrahedral-mesh block.
fens,fes  = T4block(Length, Width, Height, nLength, nWidth, nHeight)

# Define the material model.
MR = DeforModelRed3D
material = MatDeforElastIso(MR, rho, E, nu, 0.0)

# The finite element machine for the interior of the domain: we used to
# energy-stabilized tetrahedra which are quite effective for modeling thin
# structures.
femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)

# The plate is clamped along $X = 0$. Select the nodes using the box criterion.
# Set all the displacements components to zero.
nl1 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
ebc1 = FDataDict("node_list"=>nl1, "component"=>1, "displacement"=>0.0)
ebc2 = FDataDict("node_list"=>nl1, "component"=>2, "displacement"=>0.0)
ebc3 = FDataDict("node_list"=>nl1, "component"=>3, "displacement"=>0.0)

# Make the region
region1 = FDataDict("femm"=>femm, "femm_mass"=>femm)

# Make model data
modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "essential_bcs"=>[ebc1 ebc2 ebc3], "omega_shift"=>0.0, "neigvs"=>neigvs)

# Solve for the in-vacuo free vibration modes
modeldata = AlgoDeforLinearModule.modal(modeldata)

println("Dry angular frequencies: $(modeldata["omega"]) [rad/s]")
println("Reference dry ang. fre.: $(dreffs) [rad/s]")

# Determine the wet boundary faces. Both square faces of the plate are exposed
# to the fluid as it is fully immersed. We select the faces inside the boxes.
bfes = meshboundary(fes)
rl = vcat(selectelem(fens, bfes, box = [-Inf Inf -Inf Inf Height Height], inflate = tolerance),
	selectelem(fens, bfes, box = [-Inf Inf -Inf Inf 0.0 0.0], inflate = tolerance));
# These are the wet surface panels. Note that the algorithm assumes that these
# are three-node triangular finite elements.
modeldata["wet_boundary_fes"] = subset(bfes, rl);
# The only property of the fluid that matters is the mass density.
modeldata["rhow"] = rhow

# Solve for the wet-vibration modes. The solution is returned in the model data.
modeldata = AlgoVibInFluidsModule.modal(modeldata)

println("Wet angular frequencies: $(modeldata["wet_omega"]) [rad/s]")
println("Reference wet ang. fre.: $(wreffs) [rad/s]")

# Export the fluid pressure modes, and the vibration modes of the structure. The
# names of the VTK files are left up to the algorithms. We should be able to
# locate them in the current (working) folder.
modeldata = AlgoVibInFluidsModule.exportfluidpressuremode(modeldata)
modeldata["postprocessing"]["mode"] = 1:neigvs
modeldata = AlgoVibInFluidsModule.exportmode(modeldata)

true

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
isbiorthogonal(op::DictionaryOperator) = true
isorthogonal(op::DictionaryOperator) = isdiag(op) || isorthonormal(op)
isorthonormal(op::DictionaryOperator) = false

isorthonormal(::IdentityOperator) = true
isorthonormal(::IndexExtension) = true
isorthonormal(::IndexRestriction) = true
function isorthonormal(op::CompositeOperator)
    @warn "just a rough estimate on orthonormality"
    mapreduce(isorthonormal, &, components(op))
end

function isorthogonal(op::CompositeOperator)
    @warn "just a rough estimate on orthonormality"
    mapreduce(isorthogonal, &, components(op))
end


text:
@adjoint Base.literal_pow(::typeof(^), x::Number, ::Val{p}) where {p} =
  Base.literal_pow(^,x,Val(p)),
  Δ -> (nothing, Δ * conj(p * Base.literal_pow(^,x,Val(p-1))), nothing)

@adjoint Base.convert(T::Type{<:Real}, x::Real) = convert(T, x), ȳ -> (nothing, ȳ)
@adjoint (T::Type{<:Real})(x::Real) = T(x), ȳ -> (nothing, ȳ)

for T in Base.uniontypes(Core.BuiltinInts)
    @adjoint (::Type{T})(x::Core.BuiltinInts) = T(x), Δ -> (Δ,)
end

@adjoint Base.:+(xs::Number...) = +(xs...), Δ -> map(_ -> Δ, xs)

@adjoint a // b = (a // b, c̄ -> (c̄ * 1//b, - c̄ * a // b // b))

@nograd floor, ceil, trunc, round, hash

# Complex Numbers

@adjoint (T::Type{<:Complex})(re, im) = T(re, im), c̄ -> (nothing, real(c̄), imag(c̄))

@adjoint real(x::Number) = real(x), r̄ -> (real(r̄),)
@adjoint conj(x::Number) = conj(x), r̄ -> (conj(r̄),)
@adjoint imag(x::Number) = imag(x), ī -> (real(ī)*im,)

# we intentionally define these here rather than falling back on ChainRules.jl
# because ChainRules doesn't really handle nonanalytic complex functions
@adjoint abs(x::Real) = abs(x), Δ -> (real(Δ)*sign(x),)
@adjoint abs(x::Complex) = abs(x), Δ -> (real(Δ)*x/abs(x),)
@adjoint abs2(x::Number) = abs2(x), Δ -> (real(Δ)*(x + x),)


text:
export BACKSPACE, NULL,TAB, LINE_FEED,NEWLINE, FORM_FEED, RETURN,
       SPACE, STRING_DELIM, PLUS_SIGN, DELIMITER, MINUS_SIGN,
       DECIMAL_POINT, SOLIDUS, DIGIT_ZERO, DIGIT_NINE, SEPARATOR,
       LATIN_UPPER_A, LATIN_UPPER_F, LATIN_UPPER_I, BACKSLASH, LATIN_A, LATIN_B,LATIN_E,
       LATIN_F, LATIN_I, LATIN_L, LATIN_N, LATIN_R, LATIN_S, LATIN_T, LATIN_U,LATIN_Z,
       OBJECT_BEGIN, OBJECT_END, ESCAPES, REVERSE_ESCAPES,
       LEFT_PAREN, RIGHT_PAREN, LESS_THAN,
       GREATER_THAN, LEFT_CB, RIGHT_CB, LEFT_SB, RIGHT_SB,
       PERCENT,PERIOD, NUMBER_SIGN, BANG, TILDE, LATIN_UPPER_D,STREAM,ENDSTREAM,
       LATIN_UPPER_E,LATIN_UPPER_F, LATIN_UPPER_O, LATIN_UPPER_P,
       LATIN_UPPER_R,XREF, TRAILER, STARTXREF, EOF, OBJ, ENDOBJ, ispdfspace,
       ispdfdelimiter, ispdfdigit, ispdfodigit, ispdfxdigit, gethexval, getnumval, is_crorlf

# The following bytes have significant meaning in PDF
const NULL           = UInt8('\0')
const TAB            = UInt8('\t')
const LINE_FEED      = UInt8('\n')
const FORM_FEED      = UInt8('\f')
const RETURN         = UInt8('\r')
const SPACE          = UInt8(' ')

const BACKSPACE      = UInt8('\b')

const LEFT_PAREN     = UInt8('(')
const RIGHT_PAREN    = UInt8(')')
const LESS_THAN      = UInt8('<')
const GREATER_THAN   = UInt8('>')
const LEFT_CB        = UInt8('{')
const RIGHT_CB       = UInt8('}')
const LEFT_SB        = UInt8('[')
const RIGHT_SB       = UInt8(']')

const PERCENT        = UInt8('%')
const PERIOD         = UInt8('.')
const NUMBER_SIGN    = UInt8('#')
const BANG           = UInt8('!')
const TILDE          = UInt8('~')

const STRING_DELIM   = UInt8('"')
const PLUS_SIGN      = UInt8('+')
const DELIMITER      = UInt8(',')
const MINUS_SIGN     = UInt8('-')
const DECIMAL_POINT  = UInt8('.')
const SOLIDUS        = UInt8('/')
const DIGIT_ZERO     = UInt8('0')
const DIGIT_ONE      = UInt8('1')
const DIGIT_SEVEN    = UInt8('7')
const DIGIT_NINE     = UInt8('9')
const SEPARATOR      = UInt8(':')
const LATIN_UPPER_A  = UInt8('A')
const LATIN_UPPER_D  = UInt8('D')
const LATIN_UPPER_E  = UInt8('E')
const LATIN_UPPER_F  = UInt8('F')
const LATIN_UPPER_I  = UInt8('I')
const LATIN_UPPER_M  = UInt8('M')
const LATIN_UPPER_O  = UInt8('O')
const LATIN_UPPER_P  = UInt8('P')
const LATIN_UPPER_R  = UInt8('R')

const BACKSLASH      = UInt8('\\')

const LATIN_A        = UInt8('a')
const LATIN_B        = UInt8('b')
const LATIN_D        = UInt8('d')
const LATIN_E        = UInt8('e')
const LATIN_F        = UInt8('f')
const LATIN_I        = UInt8('i')
const LATIN_J        = UInt8('j')
const LATIN_L        = UInt8('l')
const LATIN_M        = UInt8('m')
const LATIN_N        = UInt8('n')
const LATIN_O        = UInt8('o')
const LATIN_R        = UInt8('r')
const LATIN_S        = UInt8('s')
const LATIN_T        = UInt8('t')
const LATIN_U        = UInt8('u')
const LATIN_X        = UInt8('x')
const LATIN_Z        = UInt8('z')

const ESCAPES = Dict(
    LEFT_PAREN   => LEFT_PAREN,
    RIGHT_PAREN  => RIGHT_PAREN,
    BACKSLASH    => BACKSLASH,
    LATIN_B      => BACKSPACE,
    LATIN_F      => FORM_FEED,
    LATIN_N      => LINE_FEED,
    LATIN_R      => RETURN,
    LATIN_T      => TAB)

#=
const REVERSE_ESCAPES = Dict(map(reverse, ESCAPES))
const ESCAPED_ARRAY = Vector{Vector{UInt8}}(256)
for c in 0x00:0xFF
    ESCAPED_ARRAY[c + 1] =
        if c == SOLIDUS
            [SOLIDUS]  # don't escape this one
        elseif c ≥ 0x80
            [c]  # UTF-8 character copied verbatim
        elseif haskey(REVERSE_ESCAPES, c)
            [BACKSLASH, REVERSE_ESCAPES[c]]
        elseif iscntrl(@compat Char(c)) || !isprint(@compat Char(c))
            UInt8[BACKSLASH, LATIN_U, hex(c, 4)...]
        else
            [c]
        end
end
=#

const XREF     =[LATIN_X,LATIN_R,LATIN_E,LATIN_F]
const TRAILER  =[LATIN_T,LATIN_R,LATIN_A,LATIN_I,LATIN_L,LATIN_E,LATIN_R]
const STARTXREF=[LATIN_S,LATIN_T,LATIN_A,LATIN_R,LATIN_T,LATIN_X,LATIN_R,LATIN_E,LATIN_F]
const EOF      =[PERCENT,PERCENT,LATIN_UPPER_E,LATIN_UPPER_O,LATIN_UPPER_F]
const OBJ      =[LATIN_O,LATIN_B,LATIN_J]
const ENDOBJ   =[LATIN_E,LATIN_N,LATIN_D,LATIN_O,LATIN_B,LATIN_J]
const STREAM   =[LATIN_S,LATIN_T,LATIN_R,LATIN_E,LATIN_A,LATIN_M]
const ENDSTREAM=[LATIN_E,LATIN_N,LATIN_D,LATIN_S,LATIN_T,LATIN_R,LATIN_E,LATIN_A,LATIN_M]

"""
Like `isspace`, but work on bytes and includes only the four whitespace
characters defined by the PDF standard: null, space, tab, line feed,
form feed, and carriage return.
"""
ispdfspace(b::UInt8) = b == NULL || b == TAB || b == LINE_FEED || b == RETURN || b == FORM_FEED || b == SPACE

ispdfdelimiter(b::UInt8) = (b == LEFT_PAREN || b == RIGHT_PAREN || b == LESS_THAN ||
                            b == GREATER_THAN || b == LEFT_CB || b == RIGHT_CB ||
                            b == LEFT_SB || b == RIGHT_SB || b == PERCENT || b == SOLIDUS)

"""
Like `isdigit`, but for bytes.
"""
ispdfdigit(b::UInt8) = (DIGIT_ZERO ≤ b ≤ DIGIT_NINE)

"""
Like `isdigit`, but for bytes.
"""
ispdfodigit(b::UInt8) = (DIGIT_ZERO ≤ b ≤ DIGIT_SEVEN)

ispdfxdigit(b::UInt8) =
    (ispdfdigit(b) || (LATIN_UPPER_A <= b <= LATIN_UPPER_F) || (LATIN_A <= b <= LATIN_F))

gethexval(b::UInt8) =   (DIGIT_ZERO <= b <= DIGIT_NINE) ? b - DIGIT_ZERO :
                        (LATIN_UPPER_A <= b <= LATIN_UPPER_F) ? b - LATIN_UPPER_A + 0xa :
                        (LATIN_A <= b <= LATIN_F) ? b - LATIN_A + 0xa :
                        throw(ErrorException(E_BAD_NUMBER))

getnumval(b::UInt8) =   (DIGIT_ZERO <= b <= DIGIT_NINE) ? b - DIGIT_ZERO :
                        throw(ErrorException(E_BAD_NUMBER))

is_crorlf(b::UInt8) = ((b == RETURN) ||(b == LINE_FEED))


text:
using Random
using DelimitedFiles

"""
    showtypetree(T, level=0)

Show the tree of subtypes for a type.
```julia
showtypetree(Number)
```
"""
function showtypetree(T, level=0)
    println("\t" ^ level, T)
    for t in subtypes(T)
        showtypetree(t, level+1)
    end
end # showtypetree(T, level=0)

"""
    get_cvi_data(data_file::String)

Get the cvi data specified by the data_file path.
"""
function get_cvi_data(data_file::String)
    # Parse the data
    data = readdlm(data_file, ',')
    data = permutedims(data)
    train_x = data[1:2, :]
    train_y = convert(Array{Int64}, data[3, :])

    return train_x, train_y
end # get_cvi_data(data_file::String)

"""
    get_bernoulli_subset(data::Array{N, 2}, labels::Array{M, 1}) where {N<:Real, M<:Int}
"""
function get_bernoulli_subset(data::Array{N, 2}, labels::Array{M, 1}, p::Float64) where {N<:Real, M<:Int}
    # Get the dimensions of the data
    dim, n_samples = size(data)

    # Get a random subsamplin of the data
    subset = randsubseq(1:n_samples, p)

    # Return the subset
    return data[:, subset], labels[subset]
end # get_bernoulli_subset(data::Array{N, 2}, labels::Array{M, 1}) where {N<:Real, M<:Int}


text:
# default_freqs(t, nw=1) = LinRange(0,0.5/mean(diff(t))-1/length(t)/2,length(t)÷2÷nw)

function default_freqs(n::Int,fs=1)
    f = FFTW.rfftfreq(n,fs)
    0:f.multiplier:f[end]
end
default_freqs(t::AbstractVector,fs=1/mean(diff(t))) = default_freqs(length(t),fs)

default_freqs(t::AbstractVector,n::Int) = default_freqs(t[1:n])

# function default_freqs(t, nw=1)
#     fs = 1/mean(diff(t))
#     N = length(t)
#     l = N/2/nw
#     range(0,stop=0.5fs-1/N/2, length=floor(Int,l))
# end

function check_freq(f)
    zerofreq = findfirst(iszero, f)
    zerofreq !== nothing && zerofreq != 1 && throw(ArgumentError("If zero frequency is included it must be the first frequency"))
    zerofreq
end

function get_fourier_regressor(t::AbstractArray{T},f::AbstractArray{T}) where T
    zerofreq = check_freq(f)
    N  = length(t)
    Nf = length(f)
    Nreg = zerofreq === nothing ? 2Nf : 2Nf-1
    # N >= Nreg || throw(ArgumentError("Too many frequency components $Nreg > $N"))
    A  = zeros(T,N,Nreg)
    sinoffset = Nf
    π2 = T(2π)
    dd = 1/(sqrt(2*length(f))) # This is to give the spectrum the same power as using DSP.periodogram
    for fn=1:Nf
        if fn == zerofreq
            sinoffset = Nf-1
        end
        @inbounds for n = 1:N
            phi        = π2*f[fn]*t[n]
            A[n,fn]    = cos(phi)*dd
            if fn != zerofreq
                A[n,fn+sinoffset] = -sin(phi)*dd
            end
        end
    end
    A, zerofreq
end

"""`x,f = ls_spectral(y,t,f=(0:((length(y)-1)/2))/length(y); λ=0)`

perform spectral estimation using the least-squares method
`y` is the signal to be analyzed
`t` is the sampling points
`f` is a vector of frequencies

The difference between `ls_spectral` and `rfft` is `abs.(rfft) = √(2π)abs.(x)`

See also `ls_sparse_spectral` `tls_spectral`
"""
function ls_spectral(y,t,f=default_freqs(t); λ=1e-10, verbose=false)
    A, zerofreq = get_fourier_regressor(t,f)
    x = fourier_solve(A,y,zerofreq,λ)
    verbose && @info("Condition number: $(round(cond(A'A), digits=2))\n")
    x, f
end

"""`x,f = ls_spectral(y,t,f,W::AbstractVector)`
`W` is a vector of weights, same length as `y`, for weighted least-squares
"""
function ls_spectral(y,t,f,W::AbstractVector; verbose=false, λ=1e-10)
    A, zerofreq = get_fourier_regressor(t,f)
    Wd = Diagonal(W)
    x = (A'Wd*A + λ*I)\(A'Wd)*y
    verbose && @info("Condition number: $(round(cond(A'*Wd*A), digits=2))\n")
    fourier2complex(x,zerofreq), f
end

"""`x,f = tls_spectral(y,t,f)`
Perform total least-squares spectral estimation using the SVD-method. See `ls_spectral` for additional help
"""
function tls_spectral(y,t,f=default_freqs(t)[1:end-1])
    zerofreq = check_freq(f)
    A, zerofreq = get_fourier_regressor(t,f)
    AA    = [A y]
    # s     = svd(AA, full=true)
    _,_,Vt = LAPACK.gesvd!('S','S',AA)
    n     = size(A,2)
    V21   = Vt[n+1,1:n]
    V22   = Vt[n+1,n+1]
    x     = -V21/V22

    fourier2complex(x,zerofreq), f
end

"""`S,f = ls_windowpsd(y,t,freqs; nw = 8, noverlap = -1, window_func=rect, estimator=ls_spectral, kwargs...)`

perform widowed spectral estimation using the least-squares method.
`window_func` defaults to `Windows.rect`

`estimator` is the spectral estimatio function to use, default is `ls_spectral`. For sparse estimation, try
`estimator = ls_sparse_spectral` See `ls_sparse_spectral` for more help. `kwargs` are passed to `estimator`.

See `ls_spectral` for additional help.
"""
function ls_windowpsd(y,t,freqs=nothing; nw = 8, noverlap = -1, window_func=rect, estimator=ls_spectral, kwargs...)
    n = length(y)÷nw
    freqs === nothing && (freqs = default_freqs(t,n))
    windows = Windows2(y,t,n,noverlap,window_func)
    nw = length(windows)
    S       = zeros(eltype(y), length(freqs))
    noverlap = windows.ys.noverlap
    poverlap = noverlap/n
    for (yi,ti) in windows
        x  = estimator(yi,ti,freqs,windows.W; kwargs...)[1]
        S .+= abs2.(x)
    end
    # @show (1+poverlap),nw
    return S./nw^2, freqs
end

"""`ls_windowcsd(y,u,t,freqs; nw = 10, noverlap = -1, window_func=rect, estimator=ls_spectral, kwargs...)`

Perform windowed cross spectral density estimation using the least-squares method.

`y` and `u` are the two signals to be analyzed and `t::AbstractVector` are their sampling points
`window_func` defaults to `Windows.rect`

`estimator` is the spectral estimatio function to use, default is `ls_spectral`. For sparse estimation, try
`estimator = ls_sparse_spectral` See `ls_sparse_spectral` for more help.

See `ls_spectral` for additional help.
"""
function ls_windowcsd(y,u,t,freqs=nothing; nw = 10, noverlap = -1, window_func=rect, estimator=ls_spectral, kwargs...)
    n = length(y)÷nw
    freqs === nothing && (freqs = default_freqs(t,n))
    S       = zeros(ComplexF64,length(freqs))
    windowsy = Windows2(y,t,n,noverlap,window_func)
    windowsu = Windows2(u,t,n,noverlap,window_func)
    nw = length(windowsy)
    noverlap = windowsy.ys.noverlap
    poverlap = max(noverlap,1)/n
    for ((y,t), (u,_)) in zip(windowsy, windowsu)
        xy = estimator(y,t,freqs,windowsy.W; kwargs...)[1]
        xu = estimator(u,t,freqs,windowsu.W; kwargs...)[1]
        # Cross spectrum
        S += xy.*conj.(xu)
    end
    return S./nw, freqs
end

# function lscohere(y,u,t,freqs, nw = 10, noverlap = -1)
#         Syu     = lswindowcsd(y,u,t,freqs, nw, noverlap)
#         Syy     = lswindowpsd(y,  t,freqs, nw, noverlap)
#         Suu     = lswindowpsd(u,  t,freqs, nw, noverlap)
#         Sch     = (abs(Syu).^2)./(Suu.*Syy);
# end

"""`ls_cohere(y,u,t,freqs; nw = 10, noverlap = -1, estimator=ls_spectral, kwargs...)`

Perform spectral coherence estimation using the least-squares method.

`estimator` is the spectral estimatio function to use, default is `ls_spectral`. For sparse estimation, try
`estimator = ls_sparse_spectral` See `ls_sparse_spectral` for more help.
See also `ls_windowcsd` and `ls_spectral` for additional help.
"""
function ls_cohere(y,u,t,freqs=nothing; nw = 10, noverlap = -1, estimator=ls_spectral, kwargs...)
    n = length(y)÷nw
    freqs === nothing && (freqs = default_freqs(t,n))
    Syy     = zeros(length(freqs))
    Suu     = zeros(length(freqs))
    Syu     = zeros(ComplexF64,length(freqs))
    windows = Windows3(y,t,u,n,noverlap,hanning)
    for (y,t,u) in windows
        xy      = estimator(y,t,freqs,windows.W; kwargs...)[1]
        xu      = estimator(u,t,freqs,windows.W; kwargs...)[1]
        # Cross spectrum
        Syu .+= xy.*conj.(xu)
        Syy .+= abs2.(xy)
        Suu .+= abs2.(xu)
    end
    Sch = abs2.(Syu)./(Suu.*Syy)
    return Sch, freqs
end

@inline _K(V,vc,gamma) = exp.(-gamma*(V.-vc).^2)

@inline function _K_norm(V,vc,gamma)
    r = _K(V,vc,gamma)
    r ./=sum(r)
end

@inline _Kcoulomb(V,vc,gamma) = _K(V,vc,gamma).*(sign.(V) .== sign.(vc))

@inline function _Kcoulomb_norm(V,vc,gamma)
    r = _Kcoulomb(V,vc,gamma)
    r ./=sum(r)
end

"""psd(se::SpectralExt)
Compute the power spectral density For a SpectralExt object

See also `ls_windowpsd_lpv`
"""
function psd(se::SpectralExt)
    rp = LPVSpectral.reshape_params(copy(se.x),length(se.w))
    return abs2.(sum(rp,dims=2))
end

"""
`ls_spectral_lpv(Y,X,V,w,Nv::Int; λ = 1e-8, coulomb = false, normalize=true)`

Perform LPV spectral estimation using the method presented in
Bagge Carlson et al. "Linear Parameter-Varying Spectral Decomposition."
See the paper For additional details.

`Y` output\n
`X` sample locations\n
`V` scheduling signal\n
`w` frequency vector\n
`Nv` number of basis functions\n
`λ` Regularization parameter\n
`coulomb` Assume discontinuity at `v=0` (useful for signals where, e.g., Coulomb friction might cause issues.)\n
`normalize` Use normalized basis functions (See paper for details).

The method will issue a warning If less than 90% of the variance in `Y` is described by the estimated model. If this is the case, try increasing either the number of frequencies or the number of basis functions per frequency. Alternatively, try lowering the regularization parameter `λ`.

See also `psd`, `ls_sparse_spectral_lpv` and `ls_windowpsd_lpv`
"""
function ls_spectral_lpv(Y::AbstractVector,X::AbstractVector,V::AbstractVector,w,Nv::Integer;  λ = 1e-8, coulomb = false, normalize=true)
    w        = w[:]
    N        = length(Y)
    Nf       = length(w)
    K        = basis_activation_func(V,Nv,normalize,coulomb)
    M(w,X,V) = vec(vec(exp.(im*w.*X))*K(V)')'
    A        = zeros(ComplexF64,N, ifelse(coulomb,2,1)*Nf*Nv)
    for n = 1:N
        A[n,:] = M(w,X[n],V[n])
    end

    params      = real_complex_bs(A,Y,λ)
    real_params = [real.(params); imag.(params)]
    AA          = [real.(A) imag.(A)]
    e           = AA*real_params-Y
    Σ           = var(e)*inv(AA'AA + λ*I)
    fva         = 1-var(e)/var(Y)
    fva < 0.9 && @warn("Fraction of variance explained = $(fva)")
    SpectralExt(Y, X, V, w, Nv, λ, coulomb, normalize, params, Σ)

end

"""ls_windowpsd_lpv(Y::AbstractVector,X::AbstractVector,V::AbstractVector,w,Nv::Integer, nw::Int=10, noverlap=0;  kwargs...)

Perform windowed psd estimation using the LPV method. A rectangular window is always used.

See `?ls_spectral_lpv` for additional help.
"""
function ls_windowpsd_lpv(Y::AbstractVector,X::AbstractVector,V::AbstractVector,w,Nv::Integer, nw::Int=10, noverlap=0;  kwargs...)
    S       = zeros(length(w))
    windows = Windows3(Y,X,V,length(Y)÷nw,noverlap,rect) # ones produces a rectangular window
    for (y,x,v) in windows
        x  = ls_spectral_lpv(y,x,v,w,Nv; kwargs...)
        rp = reshape_params(x.x,length(w))
        S += abs2.(sum(rp,dims=2))
    end
    return S

end


text:
using Test: @test_nowarn, @test
import Trixi

"""
    test_trixi_run(parameters_file; l2=nothing, linf=nothing, atol=10*eps(), rtol=0.001, parameters...)

Test Trixi by calling `Trixi.run(parameters_file; parameters...)`.
By default, only the absence of error output is checked.
If `l2` or `linf` are specified, in addition the resulting L2/Linf errors
are compared approximately against these reference values, using `atol, rtol`
as absolute/relative tolerance.
"""
function test_trixi_run(parameters_file; l2=nothing, linf=nothing, atol=200*eps(), rtol=0.001, parameters...)
  # Run basic test to ensure that there is no output to STDERR
  l2_measured, linf_measured, _ = @test_nowarn Trixi.run(parameters_file; parameters...)

  # If present, compare L2 and Linf errors against reference values
  if !isnothing(l2)
    for (l2_expected, l2_actual) in zip(l2, l2_measured)
      @test isapprox(l2_expected, l2_actual, atol=atol, rtol=rtol)
    end
  end
  if !isnothing(linf)
    for (linf_expected, linf_actual) in zip(linf, linf_measured)
      @test isapprox(linf_expected, linf_actual, atol=atol, rtol=rtol)
    end
  end
end


text:
function dispersion_function(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        dim = 2, tol::T = 1e-4, symmetry = :sphere,
        kws...) where T<:Number

    low_tol = max(1e-4, tol) # a tolerance used for a first pass with time_limit

    MM = if dim == 2
        wavematrix2D(ω, medium, species; tol=tol, kws...)
    elseif dim == 3 && symmetry == :plane
        wavematrix3DPlane(ω, medium, species; tol=tol, kws...)
    elseif dim == 3
        wavematrix3D(ω, medium, species; tol=tol, kws...)
    else error("the dimension dim = $dim has no dispersion function implemented.")
    end

    # the constraint uses keff_vec[2] < -low_tol to better specify solutions where imag(k_effs)<0
    constraint(keff_vec::Array{T}) = (keff_vec[2] < -low_tol) ? (-one(T) + exp(-T(100.0)*keff_vec[2])) : zero(T)

    detMM(keff) = det(MM(keff))
    detMM2 = if dim == 3 && symmetry != :plane
        function(keff_vec::Array{T})
            constraint(keff_vec) + sqrt(abs(detMM(keff_vec[1]+im*keff_vec[2])))
        end
    else
        function detMM2(keff_vec::Array{T})
            constraint(keff_vec) + abs2(detMM(keff_vec[1]+im*keff_vec[2]))
        end
    end

    return detMM2
end

function wavematrix2D(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        dim = 2, tol::T = 1e-4,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier::T = T(1.005),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, species; hankel_order = hankel_order, dim = dim),
        kws...) where T<:Number

    k = real(ω/medium.c)
    S = length(species)
    ho = hankel_order

    len = (2ho+1) * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    function M_component(keff,j,l,m,n)
        (n == m ? 1.0 : 0.0)*(j == l ? 1.0 : 0.0) - 2.0pi*species[l].num_density*t_vecs[l][m+ho+1]*
            kernelN(n-m,k*as[j,l],keff*as[j,l])/(k^2.0-keff^2.0)
    end

    # this matrix is needed to calculate the eigenvectors
    function MM(keff::Complex{T})
        ind2 = 1
        for l = 1:S for n = -ho:ho
            ind1 = 1
            for s = 1:S for m = -ho:ho
                MM_mat[ind1,ind2] = M_component(keff,s,l,m,n)
                ind1 += 1
        end end
            ind2 += 1
        end end
        return MM_mat
    end

    return MM
end

function wavematrix3D(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        dim = 3, tol::T = 1e-4,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier::T = T(1.005),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, species; hankel_order = hankel_order, dim = dim),
        kws...) where T<:Number

    k = real(ω/medium.c)
    S = length(species)
    ho = hankel_order
    len = (ho+1)^4 * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    function M_component(keff,Ns,l,m,l2,m2,s1,dl,dm,l1,m1,s2)::Complex{T}
        (m == dm && l == dl && m1 == m2 && l1 == l2 && s1 == s2 ? 1.0 : 0.0) +
        as[s1,s2] * species[s2].num_density * t_vecs[s1][l+ho+1] *
        sum(l3 ->
            if abs(m1-m2) <= l3
                gaunt_coefficients(l,m,dl,dm,l3,m1-m2) *
                gaunt_coefficients(l1,m1,l2,m2,l3,m1-m2) * Ns[l3+1,s1,s2]
            else
                zero(Complex{T})
            end
        , 0:ho) / (keff^2.0 - k^2.0)
    end

    # The order of the indices below is important
    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [kernelN(l3,k*as[s1,s2],keff*as[s1,s2]; dim = dim) for l3 = 0:ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S for dl = 0:ho for dm = -dl:dl for l1 = 0:ho for m1 = -l1:l1
            ind1 = 1
            for s1 = 1:S for l = 0:ho for m = -l:l for l2 = 0:ho for m2 = -l2:l2
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,m,l2,m2,s1,dl,dm,l1,m1,s2)
                ind1 += 1
            end end end end end
            ind2 += 1
        end end end end end
        return MM_mat
    end

    return MM
end

function wavematrix3D_allocate(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        dim = 3, tol::T = 1e-4,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier::T = T(1.005),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, species; hankel_order = hankel_order, dim = dim),
        kws...) where T<:Number

    k = real(ω/medium.c)
    S = length(species)
    ho = hankel_order
    len = (ho+1)^4 * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    cs = reshape(
        [
            gaunt_coefficients(l1,m1,l2,m2,l3,m3)
        for l1 = 0:ho for m1 = -l1:l1 for l2 = 0:ho for m2 = -l2:l2 for l3 = 0:ho for m3 = -l3:l3]
    , ((ho+1)^2,(ho+1)^2,(ho+1)^2));

    # NOTE that cs[n3,n2,n1] == gaunt_coefficients(l1,m1,l2,m2,l3,m3)
    ls, ms = spherical_harmonics_indices(ho);
    n_max = length(ls)

    function M_component2(keff,Ns,n,dn,n1,n2,s1,s2)::Complex{T}
        (n == dn && n1 == n2 && s1 == s2 ? 1.0 : 0.0) +
        as[s1,s2] * species[s2].num_density * t_vecs[s1][ls[n]+ho+1] *
        sum(cs[:,n,dn] .* cs[:,n1,n2] .* Ns[ls[:].+1,s1,s2]) / (keff^2.0 - k^2.0)
    end
    function MM2(keff::Complex{T})
        Ns = [kernelN(l3,k*as[s1,s2],keff*as[s1,s2]; dim = dim) for l3 = 0:ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S for dn = 1:n_max for n1 = 1:n_max
            ind1 = 1
            for s1 = 1:S for n = 1:n_max for n2 = 1:n_max
                MM_mat[ind1, ind2] = M_component2(keff,Ns,n,dn,n1,n2,s1,s2)
                ind1 += 1
            end end end
            ind2 += 1
        end end end
        return MM_mat
    end
    # keff = 1.0 + rand()*im
    # s1 = s2 = 1

    # norm([
    #     M_component2(keff,Ns,n,dn,n1,n2,s1,s2) - M_component(keff,Ns,ls[n],ms[n],ls[n2],ms[n2],s1,ls[dn],ms[dn],ls[n1],ms[n1],s2)
    # for n = 1:n_max for n2 = 1:n_max for dn = 1:n_max for n1 = 1:n_max ])
    # #
    # @time [
    #     M_component2(keff,Ns,n,dn,n1,n2,s1,s2)
    # for n = 1:n_max for dn = 1:n_max for n1 = 1:n_max for n2 = 1:n_max ];
    # #
    # @time [
    #      M_component(keff,Ns,ls[n],ms[n],ls[n2],ms[n2],s1,ls[dn],ms[dn],ls[n1],ms[n1],s2)
    # for n = 1:n_max for n2 = 1:n_max for dn = 1:n_max for n1 = 1:n_max ];
    return MM2
end

function wavematrix3DPlane(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        θ_inc::T = T(0), # incident on the halfspace z >0, with kp_vec being in the y-z plane
        φ_inc::T = T(0),
        dim = 3, tol::T = 1e-4,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier::T = T(1.005),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, species; hankel_order = hankel_order, dim = dim),
        kws...) where T<:Number

    k = real(ω/medium.c)
    S = length(species)
    ho = hankel_order
    len = (ho+1)^2 * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    Ys = spherical_harmonics(ho, θ_inc, φ_inc);
    lm_to_n = lm_to_spherical_harmonic_index

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    function M_component(keff,Ns,l,m,s1,dl,dm,s2)::Complex{T}
        (m == dm && l == dl && s1 == s2 ? one(Complex{T}) : zero(Complex{T})) +
        Complex{T}(4pi) * as[s1,s2] * species[s2].num_density * t_vecs[s1][l+ho+1] *
        sum(
            if abs(dm-m) <= l1
                Complex{T}(im)^(-l1) * Ys[lm_to_n(l1,dm-m)] * Ns[l1+1,s1,s2] *
                gaunt_coefficients(dl,dm,l,m,l1,dm-m)
            else
                zero(Complex{T})
            end
        for l1 in 0:ho) / (keff^2.0 - k^2.0)
    end

    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [kernelN(l1,k*as[s1,s2],keff*as[s1,s2]; dim = dim) for l1 = 0:ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S, dl = 0:ho for dm = -dl:dl
            ind1 = 1
            for s1 = 1:S, l = 0:ho for m = -l:l
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,m,s1,dl,dm,s2)
                ind1 += 1
            end end
            ind2 += 1
        end end
        return MM_mat
    end

    return MM
end


text:
include("default_params.jl")
include("reflecting_model/model_def.jl")

sim_data = CSV.read(
    (@__DIR__)*
    "/../hitting_times_model/reflecting_model/transient_distribution/"*
    ic*"/data/sim.csv",DataFrame)
x_vals=range(0,model.b;length=n_err_evals)

function cdf(x::AbstractVector,data::DataFrame)
    sort!(data,:X)

    x1 = data.X[data.φ.==1]
    x2 = data.X[data.φ.==2]

    c1 = 1
    c2 = 1

    F1 = zeros(length(x))
    F2 = zeros(length(x))

    for (c,xi) in enumerate(x)
        tmp1 = findfirst(y->y>xi,x1[c1:end])
        if !(tmp1===nothing)
            c1 += tmp1-1
        else
            c1 = length(x1)+1
        end
        F1[c] = c1-1

        tmp2 = findfirst(y->y>xi,x2[c2:end])
        if !(tmp2===nothing)
            c2 += tmp2-1
        else
            c2 = length(x2)+1
        end
        F2[c] = c2-1
    end
    return F1./length(data.X), F2./length(data.X)
end

F1,F2 = cdf(x_vals,sim_data)

CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/sims_cdf_evaluated.csv",DataFrame(phase_1=F1,phase_2=F2))

function bootstrap(sim_data,x_vals,n_boot)
    rng = StableRNGs.StableRNG(14129)
    n_sims = length(sim_data.X)
    for n in 1:n_boot
        ind = rand(rng,1:n_sims,n_sims)
        boot_data = sim_data[ind,:]
        F1,F2 = cdf(x_vals,boot_data)
        CSV.write(
            (@__DIR__)*"/reflecting_model/data/"*ic*"/data/bootstrap/sample_"*string(n)*".csv",
            DataFrame(phase_1=F1,phase_2=F2),
        )
    end
    return nothing
end

bootstrap(sim_data,x_vals,n_boot)


text:
import Markdown

struct TransducerLister
    m::Module
end

TransducerLister() = TransducerLister(@__MODULE__)

is_transducer_type(t) = t isa Type && t <: Transducer
is_transducer_type(::typeof(Zip)) = true

Transducer(tl::TransducerLister) =
    Filter(x -> is_transducer_type(getproperty(tl.m, x))) |>
    Map(x -> Docs.Binding(tl.m, x))

(tl::TransducerLister)() = eduction(Transducer(tl), names(tl.m))

header_code(m::Markdown.MD) =
    if !isempty(m.content)
        header_code(m.content[1])
    end
header_code(c::Markdown.Code) = c
header_code(::Any) = nothing

first_paragraph(m::Markdown.MD) =
    if length(m.content) == 1
        first_paragraph(m.content[1])
    elseif length(m.content) > 1
        first_paragraph(m.content[2])
    end
first_paragraph(p::Markdown.Paragraph) = p
first_paragraph(::Any) = nothing

hasdoc(m::Module, b::Docs.Binding) = haskey(Docs.meta(m), b)

function Base.show(io::IO, ::MIME"text/markdown", tl::TransducerLister)
    println(io, "| **Transducer** | **Summary** |")
    println(io, "|:-- |:-- |")
    foreach(tl()) do binding
        hasdoc(tl.m, binding) || return
        d = Docs.doc(binding)
        h = header_code(d)
        if h === nothing
            shown = binding.var
        else
            shown, = split(h.code, "\n", limit=2)
        end

        p = first_paragraph(d)
        if p === nothing
            gist = ""
        else
            gist, = split(sprint(show, "text/markdown", Markdown.MD(p)),
                          "\n", limit=2)
        end

        println(io,
                "| [`", shown, "`](@ref ", binding, ")",
                " | ", gist,
                " |")
    end
end

function Markdown.MD(tl::TransducerLister)
    io = IOBuffer()
    show(io, "text/markdown", tl)
    seek(io, 0)
    return Markdown.parse(io)
end

Base.show(io::IO, mime::MIME"text/plain", tl::TransducerLister) =
    _show(io, mime, tl)
Base.show(io::IO, mime::MIME"text/html", tl::TransducerLister) =
    _show(io, mime, tl)

_show(io, mime, tl) = show(io, mime, Markdown.MD(tl))


text:
module DynaPlanning

using LightGraphs, LightGraphsFlows
using GraphPlot, Plots
using Colors
using Compose, Cairo
using SparseArrays
using DataStructures
using Random
using Setfield
using LinearAlgebra
using JuMP, GLPK
using Logging
using .Threads
using Actors
import Actors: spawn

abstract type AbstractStrategy end

abstract type AbstractInstance end

abstract type AbstractAgent end

################################################################################

mutable struct DynaTuple
    e::LightGraphs.SimpleGraphs.SimpleEdge{Int64}
    strategy::AbstractStrategy
    s::LightGraphs.SimpleGraphs.SimpleEdge{Int64}
end

mutable struct DictSrv{L}
    lk::L
end

(ds::DictSrv)(f::Function, args...) = call(ds.lk, f, args...)

# Indexing interface
Base.getindex(d::DictSrv, key) = call(d.lk, getindex, key)
Base.setindex!(d::DictSrv, value, key) = call(d.lk, setindex!, value, key)

# Behavior
ds(d::Dict, f::Function, args...) = f(d, args...)

# Start dictsrv
dictsrv(d::Dict; remote = false) = DictSrv(spawn(ds, d; remote))

export DictSrv, dictsrv

mutable struct EnvironmentModel
    model::DictSrv
    EnvironmentModel() = (x = new(); x.model = dictsrv(Dict{DynaTuple,Int}()))

    function Base.isequal(A::DynaTuple, B::DynaTuple)
        A.e == B.e && A.strategy == B.strategy && A.s == B.s
    end

    function Base.hash(A::DynaTuple)
        hash(A.e + A.strategy + A.s)
    end
end

mutable struct StackelbergStrategy <: AbstractStrategy
    strategy::Array{Float64,1}
    StackelbergStrategy() = (x = new())
end

mutable struct SelfishStrategy <: AbstractStrategy
    strategy::Array{Float64,1}
    SelfishStrategy() = (x = new())
end

mutable struct Incentiver
    strategy_followed::AbstractStrategy
    Incentiver() = (x = new())
    Incentiver(AbstractStrategy) = new(AbstractStrategy)
end

mutable struct Q
    q::Array{
        Pair{LightGraphs.SimpleGraphs.SimpleEdge{Int64},AbstractStrategy},
        Pair{LightGraphs.SimpleGraphs.SimpleEdge{Int64},Float64},
    }
    Q() = (x = new())
end

mutable struct DynaAgentPlus <: AbstractAgent
    path::Queue{Edge}
    state::Edge
    actions::Array{Edge,1}              # history
    γ::Float64                          # discount rate
    α::Float64                          # step-size
    ϵ::Float64                          # exploration parameter
    episodes::Int64                     # episodes
    τ::Int64                            # revisit time steps (Dyna-Q+)
    κ::Float64                          # weight
    θ::Float64                          # priority sweeping
    queue::PriorityQueue{Tuple{LightGraphs.SimpleGraphs.SimpleEdge{Int64},AbstractStrategy}}
    predecessors::Array{
        Pair{
            LightGraphs.SimpleGraphs.SimpleEdge{Int64},
            Array{Tuple{Int64,LightGraphs.SimpleGraphs.SimpleEdge{Int64},AbstractStrategy}},
        },
    }
    model::EnvironmentModel
    q::Q
    strategy::AbstractStrategy
    time::Int64
    DynaAgentPlus() = (x = new())

    function chooseAction(a)

        max_reward = -9999
        # ϵ-greedy
        if rand(Float64, 1) <= ϵ
            # randomly choose action
            action = a[rand(1:end)]
        else
            # greedy action
            current_state = state
            if length(Set(q)) == 1
                action = rand(actions)
            else
                for a in actions
                    this_reward = q[current_state][a]
                    if this_reward >= max_reward
                        action = a
                        max_reward = this_reward
                    end
                end
            end
        end

        return action

    end

    function learn()

    end
end

mutable struct StackelbergInstance <: AbstractInstance
    network::SimpleDiGraph
    costfunct::Array{Function,1}
    β::Float32                              # fraction
    agents::Array{DynaAgentPlus,1}
    StackelbergInstance() = (x = new())
end

end # module

max_stars_repo_name: mgoar/dynaplanning

max_issues_repo_name: mgoar/dynaplanning

max_forks_repo_name: mgoar/dynaplanning


text:
#=
  Subtle rework
=#
function bin_length_long(s::AbstractString)
  binNum = parse(UInt, s)

  finNum = 0
  seq = 1

  for i in 1:binNum
    if (i == seq)
      finNum += 1
      seq *= 2
    end
  end

  return string( finNum )
end

function bin_length_short(s::AbstractString)
  binNum = parse(UInt, s)

  finNum = 0
  i = 1

  while i <= binNum
    i *= 2
     finNum += 1
  end

  return string( finNum )
end

function bin_length_shift(s::AbstractString)
  binNum = parse(UInt, s)
  i = 1
  finNum = 1

  while(i < binNum)
    i = i << 1
    finNum += 1
  end

  return string( finNum )
end


text:
using Bridge
using StaticArrays
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

struct SIR <: ContinuousTimeProcess{SVector{3,Float64}}
    β::Float64
    c::Float64
    γ::Float64
end

function Bridge.b(t, u, P::SIR)
    (S,I,R) = u
    N = S + I + R
    dS = -P.β*P.c*S*I/N
    dI = P.β*P.c*S*I/N - P.γ*I
    dR = P.γ*I
    return @SVector [dS,dI,dR]
end

function Bridge.σ(t, u, P::SIR)
    (S,I,R) = u
    N = S + I + R
    ifrac = abs(P.β*P.c*I/N*S)
    rfrac = abs(P.γ*I)
    return @SMatrix Float64[
     sqrt(ifrac)      0.0
    -sqrt(ifrac)  -sqrt(rfrac)
     0.0   sqrt(rfrac)
    ]
end

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
ts = 0.0:δt:tmax;

u0 = @SVector [990.0,10.0,0.0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob = SIR(p...);

W = sample(ts, Wiener{SVector{2,Float64}}());

sol = solve(Bridge.EulerMaruyama(), u0, W, prob);

df_sde = DataFrame(Bridge.mat(sol.yy)')
df_sde[!,:t] = ts;

@df df_sde plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark begin
    W = sample(ts, Wiener{SVector{2,Float64}}());
    solve(Bridge.EulerMaruyama(), u0, W, prob);
end


text:
import MultipleShootingPEM
ms = MultipleShootingPEM
using Plots
pyplot()
using LaTeXStrings

# Generate Data
k0 = 1
N = 200
θ_opt = 3.78
function logistic_f(x_next, x, k, θ)
    x_next .=  θ*x - θ*x.^2
    return
end
function logistic_g(y, x, k, θ)
    y .= x
    return
end
# Define initial conditions
x0 = [0.5]
y = Vector{Vector{Float64}}(N)
for i = 1:length(y)
    y[i] = zeros(1)
end
x_final = ms.simulate_space_state!(y, logistic_f, logistic_g,
                                   x0, (k0, k0+N-1), (θ_opt,))

# Multiple shooting error estimation
M = 40
k0_list = collect(k0:Int(N//M):k0+N-1)
x0_list = y[k0_list]
list_procs = ones(Int, M)

# Define polynomial model
# f(x) = θ1*x - θ1*x^2
@everywhere function f(y, dx, dθ, x, k, θ)
    y .= θ[1]*x - θ[1]*x.^2
    dx .= θ[1] - 2*θ[1]*x
    dθ .=  x-x.^2
    return
end
# g(x) = x
@everywhere function g(y, dx, dθ, x, k, θ)
    y .= x
    dx .= 1
    dθ .= 0
    return
end

θ = 0
multiple_shoot = ms.MultipleShooting(f, g, x0_list,
                                     y, k0_list, θ,
                                     list_procs)
ms.cost_function(multiple_shoot)

# Get cost function
θ_min = 2.0
θ_delta = 0.01
θ_max = 4.0
θ_range = θ_min:θ_delta:θ_max
npoints = length(θ_range)
cost = Vector{Float64}(npoints)
i = 1
for θ = θ_range
    ms.new_simulation!(multiple_shoot, x0_list, [θ])
    cost[i] = ms.cost_function(multiple_shoot)
    i += 1
end

nrealiz = 20
cost_matrix = Matrix{Float64}(npoints, nrealiz)
srand(1)
j = 1
for j = 1:nrealiz
    σ = 0.005
    x0_list_modified = copy(x0_list)
    for i = 1:M
        x0_list_modified .+= σ*randn()
    end
    i = 1
    for θ = θ_range
        ms.new_simulation!(multiple_shoot, x0_list_modified, [θ])
        cost_matrix[i, j] = ms.cost_function(multiple_shoot)
        i += 1
    end
    j += 1
end

# Get solver solution
srand(4)
θ_list = Float64[]
cost_list = Float64[]
res_list = []
for θp = 3.2:0.05:3.9
    σ = 0.002
    x0_list_modified = copy(x0_list)
    for i = 1:M
        x0_list_modified .+= σ*randn()
    end
    θ0 = [θp]
    opt = ms.OptimizationProblem(f, g, x0_list_modified, y,
                                 k0_list, θ0, list_procs)

    res = ms.solve(opt, options=Dict("gtol" => 1e-10,
                                     "xtol" => 1e-10,
                                     "maxiter" => 2000))
    θ_est = res["x"][1]
    cost_est = res["fun"]

    push!(θ_list, θ_est)
    push!(cost_list, cost_est)
    push!(res_list, res)
end

plot(θ_range, cost, color=:black, lw=2, legend=false, grid=false,
     xlims=[θ_min, θ_max], xguide=L"\theta", yguide=L"V(\theta)")
plot!(θ_range, cost_matrix, color=:blue, linealpha=0.2)
scatter!(θ_list, cost_list, marker=:c, color=:green, markersize=5)
vline!([θ_opt], ls=:dot, color=:red)

# savefig("ratio4_ms_logistic.tex")

# Evaluate solver performance
niter_list = [res["niter"] for res in res_list]
median(niter_list)
maximum(niter_list)
minimum(niter_list)

nfun_list = [res["nfev"] for res in res_list]
median(nfun_list)
maximum(nfun_list)
minimum(nfun_list)

exect_list = [res["execution_time"] for res in res_list]
median(exect_list)
maximum(exect_list)
minimum(exect_list)


text:
using QuantumInfo, RandomQuantum
using LinearAlgebra: svdvals, norm, ishermitian, eigvals
using StatsBase: var

using Test

trnorm(A) = sum(svdvals(A))

@testset "Randomized GinUE tests" begin
for _ in 1:10
    rows = rand(2:10)
    cols = rand(1:10)
    X = rand(GinibreEnsemble(rows,cols))
    @test var(vec(real(X))) > 0
    @test var(vec(imag(X))) > 0
end
end

@testset "Randomized FubiniStudyEnsemble tests" begin
for _ in 1:10
    d = rand(2:10)
    ψ = rand(FubiniStudyPureState(d))
    @test var(real(ψ)) > 0
    @test var(imag(ψ)) > 0
    @test isapprox(norm(ψ), 1.0)
end
end

@testset "Randomized FSOE tests" begin
for _ in 1:10
    d = rand(2:10)
    db = rand(2:10)
    ρ = rand(FubiniStudyMixedState(d,db))
    @test isapprox(trnorm(ρ), 1)
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end
end

@testset "Randomized HSOE tests" begin
for _ in 1:10
    d = rand(2:10)
    ρ = rand(HilbertSchmidtMixedState(d))
    @test isapprox(trnorm(ρ), 1)
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end
end

@testset "BuresMixedState" begin
for _ in 1:10
    d = rand(2:10)
    ρ = rand(BuresMixedState(d))
    @test isapprox(trnorm(ρ), 1)
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end
end

@testset "ClosedHaarEnsemble" begin
for _ in 1:10
    d = rand(2:10)
    U = rand(ClosedHaarEnsemble(d))
    @test isapprox(U*U', QuantumInfo.eye(d), atol=1e-12)
    @test isapprox(U'*U, QuantumInfo.eye(d), atol=1e-12)
    @test isapprox(abs.(eigvals(U)),  ones(d), atol=1e-12)
end
end

@testset "OpenHaarEnsemble" begin
for _ in 1:10
    d = 2 # rand(2:10)
    db = rand(2:10)
    E = rand(OpenHaarEnsemble(d,db))
    @test istp(E,tol=1e-14)
    @test iscp(E,tol=1e-14)
    @test ischannel(E,tol=1e-14)
end
end

@testset "ClosedEvolution" begin
for _ in 1:10
    α = rand()
    d = rand(2:10)
    U = rand(RandomClosedEvolution(d,α))
    @test isapprox(U*U', QuantumInfo.eye(d), atol=1e-12)
    @test isapprox(U'*U, QuantumInfo.eye(d), atol=1e-12)
    @test isapprox(abs.(eigvals(U)), ones(d), atol=1e-12)
end
end

@testset "OpenEvolution" begin
for _ in 1:10
    α = rand()
    d = rand(2:10)
    db = rand(2:10)
    E = rand(RandomOpenEvolution(d,db,α))
    @test ischannel(E,tol=1e-12)
end
end


text:
__precompile__()

module Media

include("system.jl")
include("compat.jl")

end


text:
# # Adaptive Cruise Controller (ACC)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/ACC.ipynb)
#
# The Adaptive Cruise Control (ACC) benchmark is a system that tracks a set
# velocity and maintains a safe distance from a lead vehicle by adjusting
# the longitudinal acceleration of an ego vehicle.
# The neural network computes optimal control actions while satisfying safe
# distance, velocity, and acceleration constraints.

module ACC  #jl

using NeuralNetworkAnalysis, MAT
using NeuralNetworkAnalysis: FunctionPreprocessing

# ## Model
#
# For this case study, the ego car is set to travel at a set speed $v_{set} = 30$
# and maintains a safe distance $D_{safe}$ from the lead car.  The car’s dynamics
# are described as follows:
#
# ```math
# \left\{ \begin{array}{lcl}
# \dot{x}_{lead}(t) &=& v_{lead}(t) \\
# \dot{v}_{lead}(t) &=& \gamma_{lead}(t) \\
# \dot{\gamma}_{lead}(t) &=& -2\gamma_{lead}(t) + 2a_{lead}(t) - uv_{lead}^2(t)  \\
# \dot{x}_{ego}(t) &=& v_{ego}(t) \\
# \dot{v}_{ego}(t) &=& \gamma_{ego}(t) \\
# \dot{\gamma}_{ego}(t) &=& -2\gamma_{ego}(t) + 2a_{ego}(t) - uv_{ego}^2(t)
# \end{array} \right.
# ```
# where $x_i$ is the position, $v_i$ is the velocity, $γ_i$ is the
# acceleration of the car, $a_i$ is the acceleration control input applied
# to the car, and $u = 0.0001$ is the friction parameter, where
# $i ∈ \{ego, lead\}$. For this benchmark we are given four neural network
# controllers with 3, 5, 7, and 10 hidden layers of 20 neurons each, but only evaluate
# the controller with 5 hidden layers. All of them have the same number of inputs
# $(v_{set}, T_{gap}, v_{ego}, D_{rel}, v_{rel})$ and one output ($a_{ego}$).

const u = 0.0001  # friction parameter
const a_lead = -2.0  # acceleration control input applied to the lead vehicle

@taylorize function ACC!(dx, x, p, t)
    v_lead = x[2]  # lead car velocity
    γ_lead = x[3]  # lead car acceleration
    v_ego = x[5]  # ego car velocity
    γ_ego = x[6]  # ego car acceleration
    a_ego = x[7]  # ego car acceleration control input

    ## lead car dynamics
    dx[1] = v_lead
    dx[2] = γ_lead
    dx[3] = 2 * (a_lead - γ_lead) - u * v_lead^2

    ## ego car dynamics
    dx[4] = v_ego
    dx[5] = γ_ego
    dx[6] = 2 * (a_ego - γ_ego) - u * v_ego^2
    dx[7] = zero(a_ego)
    return dx
end

# ## Specification

# The uncertain initial condition is:
#
# $x_{lead}(0) ∈ [90, 110], v_{lead}(0) ∈ [32, 32.2]$,
# $γ_{lead}(0) = γ_{ego}(0) = 0, v_{ego}(0) ∈ [30, 30.2], x_{ego} ∈ [10, 11]$
#
# The controller input is
# $(v_{set}, T_{gap}, v_{ego}, x_{lead} - x_{ego}, v_{lead} - v_{ego})$.
# To extract from the current state the last three inputs to the network
# we define a projection matrix $M$.
#
# The verification objective of this system is that given a scenario where
# both cars are driving safely, the lead car suddenly slows down with
# $a_{lead} = -2$.
# We want to check whether there is a collision in the following 5 seconds.
# A control period of 0.1 seconds is used.
#
# Formally, the safety specification can be expressed as:
# ```math
#     D_{rel} = x_{lead} - x_{ego} ≥ D_{safe},
# ```
# where $D_{safe} = D_{default} + T_{gap} * v_{ego}$,
# $T_{gap} = 1.4$ sec, and $D_{default} = 10$.
# After substitution, the specification reduces to:
# $x_{lead} - x_{ego} - T_{gap} * v_{ego} ≥ D_{default}$.

## We choose the controller with 5 hidden layers.
controller = read_nnet_mat(@modelpath("ACC", "controller_5_20.mat");
                           act_key="act_fcns");

## The initial states according to the specification are:
X₀ = Hyperrectangle(low=[90, 32, 0, 10, 30, 0],
                    high=[110, 32.2, 0, 11, 30.2, 0])
U₀ = ZeroSet(1);

## The system has 6 state variables and 1 control variable:
vars_idx = Dict(:state_vars=>1:6, :control_vars=>7)
ivp = @ivp(x' = ACC!(x), dim: 7, x(0) ∈ X₀ × U₀)
period = 0.1;  # control period

## Preprocessing function for the network input:
v_set = 30.0  # ego car's set speed
T_gap = 1.4
M = zeros(3, 6)
M[1, 5] = 1.0
M[2, 1] = 1.0
M[2, 4] = -1.0
M[3, 2] = 1.0
M[3, 5] = -1.0
function preprocess(X::LazySet)  # version for set computations
    Y1 = Singleton([v_set, T_gap])
    Y2 = linear_map(M, X)
    return cartesian_product(Y1, Y2)
end
function preprocess(X::AbstractVector)  # version for simulations
    Y1 = [v_set, T_gap]
    Y2 = M * X
    return vcat(Y1, Y2)
end
control_preprocessing = FunctionPreprocessing(preprocess)

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       preprocessing=control_preprocessing);

## Safety specification
T = 5.0  # time horizon

D_default = 10.0
d_rel = [1.0, 0, 0, -1, 0, 0, 0]
d_safe = [0, 0, 0, 0, T_gap, 0, 0]

d_prop = d_rel - d_safe
safe_states = HalfSpace(-d_prop, -D_default)
predicate = X -> X ⊆ safe_states;

# ## Results

# To integrate the ODE, we use the Taylor-model-based algorithm:
alg = TMJets(abstol=1e-6, orderT=6, orderQ=1);

# To propagate sets through the neural network, we use the `Ai2` algorithm:
alg_nn = Ai2()

function benchmark(; silent::Bool=false)
    ## We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res_sol.value
    silent || print_timed(res_sol)

    ## Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    solz = overapproximate(sol, Zonotope)
    res_pred = @timed predicate(solz)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is satisfied.")
    else
        silent || println("The property may be violated.")
    end
    return solz
end

benchmark(silent=true)  # warm-up
res = @timed benchmark()  # benchmark
sol = res.value
println("total analysis time")
print_timed(res);

# We also compute some simulations:

import DifferentialEquations

println("simulation")
res = @timed simulate(prob, T=T, trajectories=10, include_vertices=true)
sim = res.value
print_timed(res);

# Finally we plot the results:

using Plots
import DisplayAs

fig = plot(leg=(0.4, 0.3))
xlabel!(fig, "time")
F = flowpipe(sol)

fp_rel = linear_map(Matrix(d_rel'), F)
output_map_rel = d_rel

fp_safe = affine_map(Matrix(d_safe'), [D_default], F)
output_map_safe = vcat([D_default], d_safe)

plot!(fig, fp_rel, vars=(0, 1), c=:red, alpha=.4)
plot!(fig, fp_safe, vars=(0, 1), c=:blue, alpha=.4)

plot_simulation!(fig, sim; output_map=output_map_rel, color=:red, lab="Drel")
plot_simulation!(fig, sim; output_map=output_map_safe, color=:blue, lab="Dsafe")

fig = DisplayAs.Text(DisplayAs.PNG(fig))
## savefig("ACC.png")
fig

#-

end  #jl
nothing  #jl


text:
# XML nodes

abstract AbstractXMLNode

#### Types of attributes

const XML_ATTRIBUTE_CDATA = 1
const XML_ATTRIBUTE_ID = 2
const XML_ATTRIBUTE_IDREF = 3
const XML_ATTRIBUTE_IDREFS = 4
const XML_ATTRIBUTE_ENTITY = 5
const XML_ATTRIBUTE_ENTITIES = 6
const XML_ATTRIBUTE_NMTOKEN = 7
const XML_ATTRIBUTE_NMTOKENS = 8
const XML_ATTRIBUTE_ENUMERATION = 9
const XML_ATTRIBUTE_NOTATION = 10

#### Types of nodes

const XML_ELEMENT_NODE = 1
const XML_ATTRIBUTE_NODE = 2
const XML_TEXT_NODE = 3
const XML_CDATA_SECTION_NODE = 4
const XML_ENTITY_REF_NODE = 5
const XML_ENTITY_NODE = 6
const XML_PI_NODE = 7
const XML_COMMENT_NODE = 8
const XML_DOCUMENT_NODE = 9
const XML_DOCUMENT_TYPE_NODE = 10
const XML_DOCUMENT_FRAG_NODE = 11
const XML_NOTATION_NODE = 12
const XML_HTML_DOCUMENT_NODE = 13
const XML_DTD_NODE = 14
const XML_ELEMENT_DECL = 15
const XML_ATTRIBUTE_DECL = 16
const XML_ENTITY_DECL = 17
const XML_NAMESPACE_DECL = 18
const XML_XINCLUDE_START = 19
const XML_XINCLUDE_END = 20
const XML_DOCB_DOCUMENT_NODE = 21

##### Generic methods

is_elementnode(nd::AbstractXMLNode) = (nodetype(nd) == XML_ELEMENT_NODE)
is_textnode(nd::AbstractXMLNode) = (nodetype(nd) == XML_TEXT_NODE)
is_commentnode(nd::AbstractXMLNode) = (nodetype(nd) == XML_COMMENT_NODE)
is_cdatanode(nd::AbstractXMLNode) = (nodetype(nd) == XML_CDATA_SECTION_NODE)

#######################################
#
#  XML Attributes
#
#######################################

immutable _XMLAttrStruct
    # common part
    _private::Ptr{Void}
    nodetype::Cint
    name::Xstr
    children::Xptr
    last::Xptr
    parent::Xptr
    next::Xptr
    prev::Xptr
    doc::Xptr

    # specific part
    ns::Xptr
    atype::Cint
    psvi::Ptr{Void}
end

type XMLAttr
    ptr::Xptr
    _struct::_XMLAttrStruct

    function XMLAttr(ptr::Xptr)
        s = unsafe_load(convert(Ptr{_XMLAttrStruct}, ptr))
        @assert s.nodetype == XML_ATTRIBUTE_NODE
        new(ptr, s)
    end
end

name(a::XMLAttr) = bytestring(a._struct.name)

function value(a::XMLAttr)
    pct = ccall(xmlNodeGetContent, Xstr, (Xptr,), a._struct.children)
    (pct != nullptr ? _xcopystr(pct) : "")::String
end

# iterations

immutable XMLAttrIter
    p::Xptr
end

Base.start(it::XMLAttrIter) = it.p
Base.done(it::XMLAttrIter, p::Xptr) = (p == nullptr)
Base.next(it::XMLAttrIter, p::Xptr) = (a = XMLAttr(p); (a, a._struct.next))

#######################################
#
#  Base XML Nodes
#
#######################################

immutable _XMLNodeStruct
    # common part
    _private::Ptr{Void}
    nodetype::Cint
    name::Ptr{Cchar}
    children::Xptr
    last::Xptr
    parent::Xptr
    next::Xptr
    prev::Xptr
    doc::Xptr

    # specific part
    ns::Xptr
    content::Xstr
    attrs::Xptr
    nsdef::Xptr
    psvi::Ptr{Void}
    line::Cushort
    extra::Cushort
end

type XMLNode <: AbstractXMLNode
    ptr::Xptr
    _struct::_XMLNodeStruct

    function XMLNode(ptr::Xptr)
        s = unsafe_load(convert(Ptr{_XMLNodeStruct}, ptr))
        new(ptr, s)
    end
end

name(nd::XMLNode) = bytestring(nd._struct.name)
nodetype(nd::XMLNode) = nd._struct.nodetype
has_children(nd::XMLNode) = (nd._struct.children != nullptr)

# whether it is a white-space only text node
is_blanknode(nd::XMLNode) = bool(ccall(xmlIsBlankNode, Cint, (Xptr,), nd.ptr))

# iteration over children

immutable XMLNodeIter
    p::Xptr
end

Base.start(it::XMLNodeIter) = it.p
Base.done(it::XMLNodeIter, p::Xptr) = (p == nullptr)
Base.next(it::XMLNodeIter, p::Xptr) = (nd = XMLNode(p); (nd, nd._struct.next))

child_nodes(nd::XMLNode) = XMLNodeIter(nd._struct.children)

function content(nd::XMLNode)
    pct = ccall(xmlNodeGetContent, Xstr, (Xptr,), nd.ptr)
    (pct != nullptr ? _xcopystr(pct) : "")::String
end

# dumping

const DEFAULT_DUMPBUFFER_SIZE = 4096

function Base.string(nd::XMLNode)
    buf = XBuffer(DEFAULT_DUMPBUFFER_SIZE)
    ccall(xmlNodeDump, Cint, (Xptr, Xptr, Xptr, Cint, Cint),
        buf.ptr, nd._struct.doc, nd.ptr, 0, 1)
    r = content(buf)
    free(buf)
    return r
end

Base.show(io::IO, nd::XMLNode) = println(io, string(nd))

#######################################
#
#  XML Elements
#
#######################################

type XMLElement <: AbstractXMLNode
    node::XMLNode

    function XMLElement(node::XMLNode)
        if !is_elementnode(node)
            throw(ArgumentError("The input node is not an element."))
        end
        new(node)
    end

    XMLElement(ptr::Xptr) = XMLElement(XMLNode(ptr))
end

name(x::XMLElement) = name(x.node)
nodetype(x::XMLElement) = XML_ELEMENT_NODE
has_children(x::XMLElement) = has_children(x.node)
child_nodes(x::XMLElement) = child_nodes(x.node)
content(x::XMLElement) = content(x.node)

Base.string(x::XMLElement) = string(x.node)
Base.show(io::IO, x::XMLElement) = show(io, x.node)

# attribute access

function attribute(x::XMLElement, name::String; required::Bool=false)
    pv = ccall(xmlGetProp, Xstr, (Xptr, Xstr), x.node.ptr, name)
    if pv != nullptr
        return _xcopystr(pv)
    else
        if required
            throw(XMLAttributeNotFound())
        else
            return nothing
        end
    end
end

function has_attribute(x::XMLElement, name::String)
    p = ccall(xmlHasProp, Xptr, (Xptr, Xstr), x.node.ptr, name)
    return p != nullptr
end

has_attributes(x::XMLElement) = (x.node._struct.attrs != nullptr)
attributes(x::XMLElement) = XMLAttrIter(x.node._struct.attrs)

function attributes_dict(x::XMLElement)
    # make an dictionary based on attributes

    dct = Dict{String,String}()
    if has_attributes(x)
        for a in attributes(x)
            dct[name(a)] = value(a)
        end
    end
    return dct
end

# element access

immutable XMLElementIter
    parent_ptr::Xptr
end

Base.start(it::XMLElementIter) = ccall(xmlFirstElementChild, Xptr, (Xptr,), it.parent_ptr)
Base.done(it::XMLElementIter, p::Xptr) = (p == nullptr)
Base.next(it::XMLElementIter, p::Xptr) = (XMLElement(p), ccall(xmlNextElementSibling, Xptr, (Xptr,), p))

child_elements(x::XMLElement) = XMLElementIter(x.node.ptr)

# elements by tag name

function find_element(x::XMLElement, n::String)
    for c in child_elements(x)
        if name(c) == n
            return c
        end
    end
    return nothing
end

function get_elements_by_tagname(x::XMLElement, n::String)
    lst = Array(XMLElement, 0)
    for c in child_elements(x)
        if name(c) == n
            push!(lst, c)
        end
    end
    return lst
end

#######################################
#
#  XML Tree Construction
#
#######################################

function new_element(name::String)
    p = ccall(xmlNewNode, Xptr, (Xptr, Xstr), nullptr, name)
    XMLElement(p)
end

function add_child(xparent::XMLElement, xchild::XMLNode)
    p = ccall(xmlAddChild, Xptr, (Xptr, Xptr), xparent.node.ptr, xchild.ptr)
    p != nullptr || throw(XMLTreeError("Failed to add a child node."))
end

add_child(xparent::XMLElement, xchild::XMLElement) = add_child(xparent, xchild.node)

function new_child(xparent::XMLElement, name::String)
    xc = new_element(name)
    add_child(xparent, xc)
    return xc
end

function new_textnode(txt::String)
    p = ccall(xmlNewText, Xptr, (Xstr,), txt)
    XMLNode(p)
end

add_text(x::XMLElement, txt::String) = add_child(x, new_textnode(txt))

function set_attribute(x::XMLElement, name::String, val::String)
    a = ccall(xmlSetProp, Xptr, (Xptr, Xstr, Xstr), x.node.ptr, name, val)
    return XMLAttr(a)
end

set_attribute(x::XMLElement, name::String, val) = set_attribute(x, name, string(val))

function set_attributes{P<:NTuple{2}}(x::XMLElement, attrs::AbstractArray{P})
    for (nam, val) in attrs
        set_attribute(x, string(nam), string(val))
    end
end

set_attributes(x::XMLElement, attrs::Associative) = set_attributes(x, collect(attrs))

function set_attributes(x::XMLElement; attrs...)
    for (nam, val) in attrs
        set_attribute(x, string(nam), string(val))
    end
end


text:
############################################################
# jo Types - miscellaneous constractors ####################
############################################################

############################################################
# includes
include("joMiscTypesMethods/joNumber.jl")
include("joMiscTypesMethods/joPAsetup.jl")
include("joMiscTypesMethods/joHelpers.jl")


text:
function add(obj::JVector, arg0::JVector)
    return jcall(obj, "add", JVector, (JVector,), arg0)
end

function add(obj::JVector, arg0::jdouble, arg1::JVector)
    return jcall(obj, "add", JVector, (jdouble, JVector), arg0, arg1)
end

function distance(obj::Point, arg0::Point)
    return jcall(obj, "distance", jdouble, (Point,), arg0)
end

function distance1(obj::JVector, arg0::JVector)
    return jcall(obj, "distance1", jdouble, (JVector,), arg0)
end

function distance_inf(obj::JVector, arg0::JVector)
    return jcall(obj, "distanceInf", jdouble, (JVector,), arg0)
end

function distance_sq(obj::JVector, arg0::JVector)
    return jcall(obj, "distanceSq", jdouble, (JVector,), arg0)
end

function dot_product(obj::JVector, arg0::JVector)
    return jcall(obj, "dotProduct", jdouble, (JVector,), arg0)
end

function get_norm(obj::JVector)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_norm1(obj::JVector)
    return jcall(obj, "getNorm1", jdouble, ())
end

function get_norm_inf(obj::JVector)
    return jcall(obj, "getNormInf", jdouble, ())
end

function get_norm_sq(obj::JVector)
    return jcall(obj, "getNormSq", jdouble, ())
end

function get_space(obj::Point)
    return jcall(obj, "getSpace", Space, ())
end

function get_zero(obj::JVector)
    return jcall(obj, "getZero", JVector, ())
end

function is_infinite(obj::JVector)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::Point)
    return jcall(obj, "isNaN", jboolean, ())
end

function negate(obj::JVector)
    return jcall(obj, "negate", JVector, ())
end

function normalize(obj::JVector)
    return jcall(obj, "normalize", JVector, ())
end

function scalar_multiply(obj::JVector, arg0::jdouble)
    return jcall(obj, "scalarMultiply", JVector, (jdouble,), arg0)
end

function subtract(obj::JVector, arg0::JVector)
    return jcall(obj, "subtract", JVector, (JVector,), arg0)
end

function subtract(obj::JVector, arg0::jdouble, arg1::JVector)
    return jcall(obj, "subtract", JVector, (jdouble, JVector), arg0, arg1)
end

function to_string(obj::JVector, arg0::NumberFormat)
    return jcall(obj, "toString", JString, (NumberFormat,), arg0)
end


text:
module SignalAnalysis

using Requires
using DocStringExtensions

using SignalBase
using SignalBase.Units
export nframes, nchannels, sampletype, framerate, duration

include("signals.jl")

include("basic.jl")
include("dsp.jl")
include("generate.jl")
include("array.jl")
include("rand.jl")

function __init__()
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    include("plot.jl")
    include("iplot.jl")
  end
end

end # module


text:
module TestTreeSet

using Phylo
using DataFrames
#using JuliaDB

using Compat.Test

species = ["Dog", "Cat", "Human"]
ntips = 10
df = DataFrame(species = species, count=[10, 20, 3])
observations = ["Dog", "Cat", "Dog", "Dog"]
#jdb = table(@NT(species = observations, count = [1, 2, 3, 4]))
jdb = DataFrame(species = observations, count = [1, 2, 3, 4])

@testset "TestSet" begin
    @test length(rand(Ultrametric(ntips), 10)) ==
        length(rand(Nonultrametric(ntips), 10))
    ts = rand(Ultrametric(ntips), 1:10)
    @test length(treeiter(ts)) == length(treenameiter(ts)) == ntrees(ts) == 10
    @test length(getbranchnames(ts)) == ntips * 2 - 2
    for name in treenameiter(ts)
        @test name ∈ 1:10
        @test length(getleafnames(ts[name])) == ntips
    end
    for tree in treeiter(ts)
        @test length(getnodenames(tree)) == ntips * 2 - 1
    end

    for tree in ts
        @test length(getbranchnames(tree)) == ntips * 2 - 2
        @test nleaves(tree) == ntips
    end
    @test nleaves(ts) == ntips
    @test getleafnames(TreeSet(Dict{String, NamedTree}(), Dict{String, Any}())) == String[]
    @test getbranchnames(TreeSet(Dict{String, NamedTree}(), Dict{String, Any}())) == Int[]
    @test nleaves(rand(Ultrametric(df), 20)) == length(species)
    @test nleaves(rand(Ultrametric(jdb), 10)) ==
        length(unique(observations))
    @test getleafinfo(rand(Ultrametric(df))) == df
    @test getleafinfo(rand(Ultrametric(jdb), 2)) == jdb
    @test nleaves(rand(Nonultrametric(df), 20)) == length(species)
    @test nleaves(rand(Nonultrametric(jdb), 10)) ==
        length(unique(observations))
    @test getleafinfo(rand(Nonultrametric(df))) == df
    @test getleafinfo(rand(Nonultrametric(jdb), 2)) == jdb

    @test nodetype(ts) == nodetype(ts[1])
    @test branchtype(ts) == branchtype(ts[1])
    ts0 = rand(Nonultrametric(ntips), 0)
    @test getleafnames(ts0) == String[]
    @test getnodenames(ts0) == String[]
    @test isempty(getbranchnames(ts0))
    @test isempty(getleafinfo(ts0))
    @test nleaves(ts0) == 0
end
end


text:
"""
Resistor Abstrato
"""
abstract type AbstractResistor end

doc"""
    R = Thermistor(R₀, B, T₀)

Modela um Termistor NTC (e também um resistor constante);
Neste construtor, `R₀` é a resistência na temperatura de referência `T₀`.
No construtor a temperatura de referência deve estar em °C, por conveniência,
mas é armazenado em K. `B` é um coeficiente que caracteriza como o termistor
varia de resistência com a temperatura. No modelo, utilizado, a resistência
varia com a temperatura de acordo com o seguinte modelo:

```math
R = R₀ \exp\left(\frac{1}{T₀} - \frac{1}{T}\right)
```

onde
 * R₀ é a resistência na temperatura de referência T₀
 * T₀ é a temperatura de referência em K, geralmente vale  293.15 K
 * T é a temperatura, em K, na qual se quer calcular a resistência
 * B é um coeficiente empírico que caracateriza a variação da resistência com a  temperatura (unidade 1/K)

    Construtor de objetos `Thermistor`. Cuidado que no construtor a unidade de `T₀` é °C e não K.

# Examples
```jldoctest
julia> R = Thermistor(5e3, 3200, 20) # Criar um objeto Thermistor
THW.Thermistor(5000.0, 3200.0, 293.15)

julia> resistance(R) # Resistência na temperatura de referência
5000.0

julia> resistance(R, 25) # Reistência a 25°C
4163.587774917559

julia> R() # Igual a resistance(R)
5000.0

julia> R(25) # Igual a resistance(R, 25)
4163.587774917559

julia> temperature(R)
20.0

julia> temperature(R, 4163.588)
24.999998498264233

```

"""
struct Thermistor <: AbstractResistor
    "Resistência na temperatura  T₀"
    R₀::Float64
    "Coeficiente empírico que modela o comportamento do termistor unidade 1/K"
    B::Float64
    "Temperatura de referência em K"
    T₀::Float64
    Thermistor(R₀=5e3, B=0.0, T₀=25.0) = new(R₀, B, T₀+273.15)
end

"""
    resistance(th::Thermistor)
    resistance(th::Thermistor, Tc)

Calcula a resistência a qualquer temperatura.
Se a temperatura não for fornecida retorna a resistência de referência `R₀`

Ver [`Thermistor`](@ref) para ver mais detalhes da classe `Thermistor`.

"""
resistance(th::Thermistor) = th.R₀
resistance(th::Thermistor, T) = th.R₀ * exp( -th.B * (1/th.T₀ - 1/(T+273.15) ) )

(th::Thermistor)(T) = resistance(th, T)
(th::Thermistor)() = th.R₀

"""
    temperature(th::Thermistor, R)
    temperature(th::Thermistor)

Dada uma resitência, calcular a temperatura correspondente. Caso não seja fornecida a resistência, retornar a temperatura de referência.

"""
temperature(th::Thermistor) = th.T₀-273.15
temperature(th::Thermistor, R) = 1/( 1/th.T₀ + 1/th.B * log(R/th.R₀) ) - 273.15

"""
    r = Resistor(R0=1e3, α=0.0, T₀=20.0)

Um resistor cuja resistência varia linearmente com a temperatura segundo a seguinte equação:

```math
R = R₀ (1 + α (T - T₀))
```

**IMPORTANTE** No construtor o coeficiente `α` entra como um *porcentual*
 linear.

# Exemples

```jldoctest
julia> R = Resistor(1e3, 0.01)
THW.Resistor(1000.0, 0.0001, 20.0)

julia> resistance(R)
1000.0

julia> resistance(R, 120)
1010.0

julia> R()
1000.0

julia> R(120)
1010.0

julia> temperature(R)
20.0

julia> temperature(R, 1010)
120.00000000000009

```
"""
struct Resistor <: AbstractResistor
    "Resistência na temperatura de referência"
    R₀::Float64
    "Coeficiente linear da resistência"
    α::Float64
    "Temperatura de referência"
    T₀::Float64
    Resistor(R0=1e3, α=0.0, T₀=20.0) = new(R0, α/100, T₀)
end

resistance(r::Resistor) = r.R₀
resistance(r::Resistor, T) = r.R₀ * (1 +  r.α*(T - r.T₀))
(r::Resistor)(T) = resistance(r, T)
(r::Resistor)() = resistance(r)

temperature(r::Resistor) = r.T₀
temperature(r::Resistor, R) = 1/r.α * (R/r.R₀ - 1) + r.T₀


text:
"""
    se_kernel_base(λ, δ)

Squared exonential GP kernel (~Gaussian). Infinitely mean square differentiable
(a.k.a. very smooth).

# Arguments
- `λ::Number`: The kernel lengthscale
- `δ::Number`: The difference between the inputs (e.g. `t1 - t2`)
"""
function se_kernel_base(λ::Number, δ::Number)
    return exp(-δ * δ / (2 * λ * λ))
end

# """
# Periodic kernel (for random cyclic functions)
# SE kernel where δ^2 -> 4 sin(π δ/P)^2
# """
# function periodic_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
#     P, λ = hyperparameters
#
#     sin_τ = sin(π * δ / P)
#     return exp(-2 * sin_τ * sin_τ / (λ * λ))
# end

# "Quasi-periodic kernel"
# function quasi_periodic_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 3 "incompatible amount of hyperparameters passed"
#     SE_λ, P_P, P_λ = hyperparameters
#
#     return se_kernel_base(SE_λ, δ) * periodic_kernel_base([P_P, P_λ], δ)
# end

# "Ornstein–Uhlenbeck (Exponential) kernel"
# function exp_kernel_base(λ::Number, δ::Number)
#     return exp(-abs(δ) / λ)
# end

# "Exponential-periodic kernel"
# function exp_periodic_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 3 "incompatible amount of hyperparameters passed"
#     OU_λ, P_P, P_λ = hyperparameters
#
#     return ou_kernel_base(OU_λ, δ) * periodic_kernel_base([P_P, P_λ], δ)
# end

# "general Matern kernel"
# function matern_kernel_base(λ::Number, δ::Number, nu::Number)
#
#     #limit of the function as it apporaches 0 (see https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function)
#     if δ == 0
#         return kernel_amplitude * kernel_amplitude
#     else
#         x = (sqrt(2 * nu) * abs(δ)) / λ
#         return ((2 ^ (1 - nu)) / (gamma(nu))) * x ^ nu * besselk(nu, x)
#     end
# end

# "Matern 3/2 kernel"
# function matern32_kernel_base(λ::Number, δ::Number)
#     x = sqrt(3) * abs(δ) / λ
#     return (1 + x) * exp(-x)
# end

"""
    matern52_kernel_base(λ, δ)

Matern 5/2 kernel. Twice mean square differentiable

# Arguments
- `λ::Number`: The kernel lengthscale
- `δ::Number`: The difference between the inputs (e.g. `t1 - t2`)
"""
function matern52_kernel_base(λ::Number, δ::Number)
    x = sqrt(5) * abs(δ) / λ
    return (1 + x * (1 + x / 3)) * exp(-x)
end

# "Matern 7/2 kernel"
# function matern72_kernel_base(λ::Number, δ::Number)
#     x = sqrt(7) * abs(δ) / λ
#     return (1 + x * (1 + x * (2 / 5 + x / 15))) * exp(-x)
# end

# "Matern 9/2 kernel"
# function matern92_kernel_base(λ::Number, δ::Number)
#     x = 3 * abs(δ) / λ
#     return (1 + x * (1 + x * (3 / 7 + x * (2 / 21 + x / 105)))) * exp(-x)
# end

"""
    pp_kernel_base(λ, δ)

Piecewise polynomial kernel. Twice mean square differentiable. Equation 4.21 in
Rasmussen and Williams.

# Arguments
- `λ::Number`: The kernel lengthscale, and cutoff variable
- `δ::Number`: The difference between the inputs (e.g. `t1 - t2`)
"""
function pp_kernel_base(λ::Number, δ::Number)
    # D = 1  # 1 dimension
    # q = 2  # twice MS differentiable
    # j = q + 1 + floor(D / 2) = 3
    r = abs(δ) / λ
    return (1 - r) ^ 5 * (8 * r * r + 5 * r + 1)
end

# """
# Rational Quadratic kernel
# Equivalent to adding together SE kernels with the inverse square of the
# lengthscales (τ = SE_λ^-2) are distributed as a Gamma distribution of p(τ|k,θ)
# where k (sometimes written as α) is the shape parameter and θ is the scale
# parameter. When α→∞, the RQ is identical to the SE with λ = (k θ)^-1/2.
# """
# function rq_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
#     k, θ = hyperparameters
#
#     return (1 + δ * δ * θ / 2) ^ -k
#     # return (1 / θ + δ * δ / 2) ^ -k
# end
# """
# Rational Quadratic kernel
# Equivalent to adding together SE kernels with the inverse square of the
# lengthscales (τ = SE_λ^-2) are distributed as a Gamma distribution of p(τ|α,β)
# where α (sometimes written as k) is the shape parameter and β is the rate
# parameter. When α→∞, the RQ is identical to the SE with λ = (α / β)^-1/2.
# """
# function rq_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
#     α, β = hyperparameters
#
#     # return (1 + δ * δ / (2 * β)) ^ -α
#     return (β + δ * δ / 2) ^ -α
# end
"""
    rq_kernel_base(hyperparameters, δ)

Rational Quadratic kernel. Equivalent to adding together SE kernels with the
inverse square of the lengthscales (τ = SE_λ^-2) are distributed as a Gamma
distribution of p(τ|α,μ) where α (sometimes written as k) is the shape parameter
and μ is the mean of the distribution. When α→∞, the RQ is identical to the SE
with λ = μ^-1/2.

# Arguments
- `hyperparameters::Vector`: The kernel shape parameter and mean (i.e. `[α, μ]`)
- `δ::Number`: The difference between the inputs (e.g. `t1 - t2`)
"""
function rq_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)

    @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
    α, μ = hyperparameters

    return (1 + δ * δ * μ / (2 * α)) ^ -α
    # return (α / μ + δ * δ / 2) ^ -α
end

"""
    rq_kernel_base(hyperparameters, δ)

Rational Matern 5/2 kernel. Equivalent to adding together Matern 5/2 kernels
with the inverse of the lengthscale (τ = M52_λ^-1) are distributed as a Gamma
distribution of p(τ|α,μ) where α (sometimes written as k) is the shape parameter
and μ is the mean of the distribution.

# Arguments
- `hyperparameters::Vector`: The kernel shape parameter and mean (i.e. `[α, μ]`)
- `δ::Number`: The difference between the inputs (e.g. `t1 - t2`)
"""
function rm52_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
    @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
    α, μ = hyperparameters

    x = sqrt(5) * abs(δ)
    y = x + α / μ
    return y ^ -α * (α * (α + x * (2 + α) * μ) + x * x * (1 + α) * (3 + α) * μ * μ / 3) / (y * y * μ * μ)
end

# """
# Periodic kernel (for random cyclic functions)
# RQ kernel where δ^2 -> 4 sin(π δ/P)^2
# """
# function periodic_rq_kernel_base(hyperparameters::Vector{<:Number}, δ::Number)
#
#     @assert length(hyperparameters) == 3 "incompatible amount of hyperparameters passed"
#     P, λ, α = hyperparameters
#
#     sin_τ = sin(π * δ / P)
#     return (1 + 2 * sin_τ * sin_τ / (α * λ * λ)) ^ -α
# end

# """
# Bessel (function of the first kind) kernel
# Bessel functions of the first kind, denoted as Jα(x), are solutions of Bessel's
# differential equation that are finite at the origin (x = 0) for integer or positive α
# http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/#bessel
# """
# function bessel_kernel_base(hyperparameters::Vector{<:Number}, δ::Number; nu=0)
#
#     @assert length(hyperparameters) == 2 "incompatible amount of hyperparameters passed"
#     λ, n = hyperparameters
#
#     @assert nu >= 0 "nu must be >= 0"
#
#     return besselj(nu + 1, λ * δ) / (δ ^ (-n * (nu + 1)))
# end

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
abstract type TypeLevel{T} end
struct TVal{T, Val} <: TypeLevel{T} end
struct TApp{Ret, Fn, Args} <: TypeLevel{Ret} end
struct TCons{T, Hd, Tl} <: TypeLevel{Cons{T}} end
struct TNil{T} <: TypeLevel{Nil{T}} end

function interpret(t::Type{TNil{T}}) where T
    nil(T)
end

function interpret(t::Type{TVal{T, V}}) where {T, V}
    V
end

function interpret(t::Type{TCons{T, Hd, Tl}}) where {T, Hd, Tl}
    tl :: List{T} = from_type(Tl)
    cons(from_type(Hd), tl)
end

function interpret(t::Type{TApp{Ret, Fn, Args}}) where {Fn, Args, Ret}
    args = from_type(Args)
    Fn(args...) :: Ret
end

Base.show(io::IO, t::Type{<:TypeLevel}) = show_t(io, t)
show_t(io::IO, @nospecialize(t)) = begin
    @match t begin
        ::Type{TypeLevel{L}} where L => print(io, "TypeLevel{",  L, "}")
        ::Type{TypeLevel} => print(io, "TypeLevel")
         _ => begin
            print(io, "TypeEncoding(")
            show_repr(io, t)
            print(io, ")")
         end
    end
end

@trait Typeable{T} begin
    to_type    :: T => Type{<:TypeLevel{T}}
    to_type(@nospecialize(x)) = TVal{T, x}
    from_type  :: Type{<:TypeLevel{T}} => T
    from_type(@nospecialize(t)) = interpret(t)

    show_repr :: [IO, Type{<:TypeLevel{T}}] => Nothing
    show_repr(io::IO, @nospecialize(t)) = begin
        print(io, from_type(t))
    end
end

to_typelist(@nospecialize(many)) =
    let T = eltype(many)
        foldr(many, init=TNil{T}) do each, prev
            TCons{T, to_type(each), prev}
        end
    end

types_to_typelist(@nospecialize(many)) =
    let T = eltype(many)
        foldr(many, init=TNil{T}) do each, prev
            TCons{T, each, prev}
        end
    end

@implement Typeable{L} where {T, L <: List{T}} begin
    to_type(@nospecialize(x)) = to_typelist(T[x...])
end

@implement Typeable{L} where {T, L <: Vector{T}} begin
    to_type(@nospecialize(x)) =
        let args = to_typelist(x),
            apply(@nospecialize(args...)) = T[args...]
            TApp{L, apply, args}
        end
end

@implement Typeable{Expr} begin
    function to_type(x::Expr)
        @when Expr(args...) = x begin
            args = to_typelist(args)
            f  = Expr
            TApp{Expr, f, args}
        @otherwise
            error("impossible")
        end
    end
end

@implement Typeable{LineNumberNode} begin
    function to_type(@nospecialize(ln))
        f = LineNumberNode
        args = Any[ln.line, ln.file] |> to_typelist
        TApp{LineNumberNode, f, args}
    end
end

@implement Typeable{QuoteNode} begin
    function to_type(@nospecialize(x))
        f = QuoteNode
        args = [x.value] |> to_typelist
        TApp{QuoteNode, f, args}
    end
end

@implement Typeable{Tp} where Tp <: Tuple  begin
    function to_type(@nospecialize(x))
        args = collect(x) |> to_typelist
        TApp{Tp, tuple, args}
    end
end

const named_tuple_maker(p...) = (;p...)

@implement Typeable{NamedTuple{Ks, Ts}} where {Ks, Ts} begin
    function to_type(@nospecialize(x))
        f = named_tuple_maker
        args = [kv for kv in zip(Ks, values(x))] |> to_typelist
        TApp{NamedTuple{Ks, Ts}, f, args}
    end
end

@implement Typeable{Symbol}
@implement Typeable{T} where T <: Number
@implement Typeable{T} where T <: Type
@implement Typeable{T} where T <: Function
@implement Typeable{Nothing}

@implement Typeable{String} begin
    function to_type(x::String)
        wrapped = Symbol(x) |> to_type
        TVal{String, wrapped}
    end
    function from_type(::Type{TVal{String, V}}) where V
        string(from_type(V))
    end
end

using Base.Threads: lock, unlock, SpinLock
const _modules = Module[]
const _lock = SpinLock()
function module_index(m::Module)
    lock(_lock)
    try
        i = findfirst(==(m), _modules)
        if i === nothing
            # TODO: thread safe
            push!(_modules, m)
            i = length(_modules)
        end
        i
    finally
        unlock(_lock)
    end
end

@implement Typeable{Module} begin
    function to_type(x::Module)
        TVal{Module, module_index(x)}
    end
    function from_type(:: Type{TVal{Module, V}}) where V
        _modules[V]
    end
end

@implement Typeable{GlobalRef} begin
    function to_type(x::GlobalRef)
        TVal{GlobalRef, to_type((x.mod, x.name))}
    end
    function from_type(::Type{TVal{GlobalRef, V}}) where V
        mod, v = from_type(V)
        GlobalRef(mod, v)
    end
end


text:
function simpson_integration(x_vec,y_vec)

    sum = 0.
    dx = x_vec[2]-x_vec[1]

    #First point
    sum += 3*y_vec[1]/8 + 3*y_vec[end]/8
    sum += 7*y_vec[2]/6 + 7*y_vec[end-1]/6
    sum += 23*y_vec[3]/24 + 23*y_vec[end-2]/24

    for idx in 4:(length(x_vec)-3)
        sum += y_vec[idx]
    end
    return sum*dx
end

function simpson_integration_2D(x_vec,y_vec,z_vec)

    sum_vector = zeros(length(y_vec))

    for idx in 1:length(sum_vector)
        sum_vector[idx] = simpson_integration(x_vec,z_vec[:,idx])
    end

    sum = simpson_integration(y_vec , sum_vector)

    return sum
end

function get_field_energy( fname::String, every::Int , tlimit::Float64)

    #Get Variables
    fid         = h5open( fname , "r")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Discretize space
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    #Create Differential Operators
    Dx  =   Diff_Operator_2D(1 ,1, [dx,dy] , [Nx,Ny])
    Dy  =   Diff_Operator_2D(1 ,2, [dx,dy] , [Nx,Ny])

    #Get max time iteration
    _ , tmax , max_sim_iter  = get_time_variables(fname)
    max_plot_iter = max_sim_iter

    if(tlimit <=  tmax)
        max_plot_iter = floor(Int , max_sim_iter * tlimit/tmax )
    else
        println("ATTENTION: max time input is larger than simulation time, using tmax")
    end

    #Potential Matrix
   # Vpot = potential_matrix(p)

    # Total volume of space
    Volume = ((ymax-ymin)*(xmax-xmin))
    dV = dx*dy
    I_factor = dV / Volume

    #Create storage matrix
    Energy_matrix = zeros( Nx , Ny )

    #Storage vectors
    E_vector = Float64[]
    t_vector = Float64[]

    t = 0
    E = 0

    #Loop to get energy
    for i in 0:every:max_plot_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(fname,i)

        #Gradient term (possibly wrong)
         vx = reshape( Dx * reshape(ψ, Nx*Ny) , (Nx,Ny))
         vy = reshape( Dy * reshape(ψ, Nx*Ny) , (Nx,Ny))

        #Potential term
        # v2 = Vpot * reshape(ψ, Nx*Ny)

        Energy_matrix = vx.^2 .+ vy.^2 .+ dψ.^2

        E = simpson_integration_2D(x , y ,Energy_matrix)

        #Add to storage vectors  and normalize
        push!(t_vector, t)
        push!(E_vector, I_factor*E  )

    end

    return t_vector , E_vector

end


text:
import CompactBases: nel, complex_rotate, real_locs

@testset "FE-DVR" begin
    include("basics.jl")
    include("complex_scaling.jl")
    include("block_structure.jl")
    include("scalar_operators.jl")
    include("inner_products.jl")
    include("function_interpolation.jl")
    include("derivatives.jl")
end


text:
using DiffEqBase, DiffEqJump, DataStructures
using Test, Statistics

using StableRNGs
rng = StableRNG(12345)

# using Plots; plotlyjs()
doplot = false

# using BenchmarkTools
# dobenchmark = false

dotestmean   = true
doprintmeans = true

# SSAs to test
SSAalgs = (RDirect(), RSSACR(),Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR(), MNRM())

Nsims        = 32000
tf           = .01
u0           = [200, 100, 150]
expected_avg = 84.876015624999994
reltol       = .01

# MODEL SETUP

# using DiffEqBiological
# rs = @reaction_network dtype begin
#     k1, 2A --> B
#     k2, B --> 2A
#     k3, A + B --> C
#     k4, C --> A + B
#     k5, 3C --> 3A
# end k1 k2 k3 k4

# model using mass action jumps
# ids: A = 1, B = 2, C = 3
reactstoch = [
    [1 => 2],
    [2 => 1],
    [1 => 1, 2 => 1],
    [3 => 1],
    [3 => 3]
]
netstoch = [
    [1 => -2, 2 => 1],
    [1 => 2, 2 => -1],
    [1 => -1, 2 => -1, 3 => 1],
    [1 => 1, 2 => 1, 3 => -1],
    [1 => 3, 3 => -3]
]
rates = [1., 2., .5, .75, .25]
spec_to_dep_jumps = [[1,3],[2,3],[4,5]]
jump_to_dep_specs = [[1,2],[1,2],[1,2,3],[1,2,3],[1,3]]
majumps = MassActionJump(rates, reactstoch, netstoch)

# average number of proteins in a simulation
function runSSAs(jump_prob)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[1,end]
    end
    mean(Psamp)
end

# TESTING:
prob = DiscreteProblem(u0, (0.0, tf), rates)

# plotting one full trajectory
if doplot
    for alg in SSAalgs
        local jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
        local sol = solve(jump_prob, SSAStepper())
        local plothand = plot(sol, seriestype=:steppost, reuse=false)
        display(plothand)
    end
end

# test the means
if dotestmean
    means = zeros(Float64,length(SSAalgs))
    for (i,alg) in enumerate(SSAalgs)
        local jump_prob = JumpProblem(prob, alg, majumps, save_positions=(false,false), vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
        means[i]  = runSSAs(jump_prob)
        relerr = abs(means[i] - expected_avg) / expected_avg
        if doprintmeans
            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ", relerr)
        end

        # if dobenchmark
        #      @btime (runSSAs($jump_prob);)
        # end

        @test abs(means[i] - expected_avg) < reltol*expected_avg
    end
end

# benchmark performance
# if dobenchmark
#     # exact methods
#     for alg in SSAalgs
#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
#         @btime solve($jump_prob, SSAStepper())
#     end
#     println()
# end

# add a test for passing MassActionJumps individually (tests combining)
if dotestmean
    majump_vec = Vector{MassActionJump{Float64,Vector{Pair{Int,Int}}}}()
    for i = 1:length(rates)
        push!(majump_vec, MassActionJump(rates[i], reactstoch[i], netstoch[i]))
    end
    jset = JumpSet((),(),nothing,majump_vec)
    jump_prob = JumpProblem(prob, Direct(), jset, save_positions=(false,false), vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
    meanval = runSSAs(jump_prob)
    relerr = abs(meanval - expected_avg) / expected_avg
    if doprintmeans
        println("Using individual MassActionJumps; Mean from method: ", typeof(Direct()), " is = ", meanval, ", rel err = ", relerr)
    end
    @test abs(meanval - expected_avg) < reltol*expected_avg
end


text:
@testset "Show" begin
    Mod.@polyvar x y z
    @test sprint(show, (x*y^2 + x + 1 + y)) == "xy² + x + y + 1"
    @test sprint(show, (x + 1 + y) / x^2) == "(x + y + 1) / (x²)"
    @test sprint(show, (x - y - x + y) / (x^2 - x)) == "(0) / (x² - x)"
    cp = CustomPoly(1 + x)
    @test sprint(show, cp) == "$(typeof(cp))(x + 1)"
    # Test taken from TypedPolynomials
    @test sprint(show, x) == "x"
    @test sprint(show, x^0) == "1"
    @test sprint(show, x^1) == "x"
    @test sprint(show, x^2) == "x²"
    @test sprint(show, 1x^2) == "x²"
    @test sprint(show, 5x) == "5x"
    @test sprint(show, x * y) == "xy"
    @test sprint(show, y * x) == "xy"
    @test sprint(show, 5 + y + x) == "x + y + 5"
    @test sprint(show, y + 5 + x) == "x + y + 5"
    @test sprint(show, x + x^2) == "x² + x"
    @test sprint(show, x^2 + x) == "x² + x"
    @test sprint(show, x^2 - 3.0x) == "x² - 3.0x"
    @test sprint(show, -x^3) == "-x³"
    @test sprint(show, -x^3 + y^2 - x) == "-x³ + y² - x"
    @test sprint(show, -2.0x^2) == "-2.0x²"
    @test sprint(show, (1.0 + 3.1im) * x*z) == "(1.0 + 3.1im)xz"
    @test sprint(show, -(1.0 + 3.1im) * z*x) == "(-1.0 - 3.1im)xz"
    @test sprint(show, x^2 + (1.0 + 3.1im) * x) == "x² + (1.0 + 3.1im)x"
    @test sprint(show, x^2 - (1.0 + 3.1im) * x) == "x² + (-1.0 - 3.1im)x"

    Mod.@polyvar x[0:9]
    @test sprint(show, sum(i*x[i]^i for i=1:10)) == "10x₉¹⁰ + 9x₈⁹ + 8x₇⁸ + 7x₆⁷ + 6x₅⁶ + 5x₄⁵ + 4x₃⁴ + 3x₂³ + 2x₁² + x₀"
    @test sprint(show, "text/latex", sum(i*x[i]^i for i=1:10)) == "\$\$ 10x_{9}^{10} + 9x_{8}^{9} + 8x_{7}^{8} + 7x_{6}^{7} + 6x_{5}^{6} + 5x_{4}^{5} + 4x_{3}^{4} + 3x_{2}^{3} + 2x_{1}^{2} + x_{0} \$\$"
    @test sprint(show, "text/latex", (x[2] + 1) / x[3]^2) == "\$\$ \\frac{x_{1} + 1}{x_{2}^{2}} \$\$"

    Mod.@polyvar x[1:11]
    @test sprint(show, "text/latex", x[10]) == "\$\$ x_{10} \$\$"
    @test sprint(show, x[10]) == "x₁₀"

    @test sprint(print, 2x[1]^2+3x[3]+1+x[4]) == "2*x[1]^2 + 3*x[3] + x[4] + 1"
end