text:
# In this tutorial, we will learn:
#
#   * How to formulate the weak form for a scalar time-harmonic electromagnetic problem
#   * How to implement a perfectly matched layer (PML) to absorb outgoing waves
#   * How to impose periodic boundary conditions in Gridap
#   * How to discretize PDEs with complex-valued solutions
#
# ## Problem statement
#
# We are going to solve a scalar electromagnetic wave scattering problem: a plane wave (Hz-polarized $H_{inc}$) scattering of a dielectric cylinder (of radius $R$ and permittivity $\varepsilon$), as illustrated below. The computational cell is of height $H$ and length $L$, and we employ a perfectly matched layer (PML) thickness of $d_{pml}$ to implement outgoing (radiation) boundary conditions for this finite domain.
#
# ![](../assets/emscatter/Illustration.png)
#
# From Maxwell's equations, considering a time-harmonic electromagnetic field, we can derive the governing equation of this problem in 2D (Helmholtz equation) [1]:
#
# ```math
# \left[-\nabla\cdot\frac{1}{\varepsilon(x)}\nabla -k^2\mu(x)\right] H = f(x),
# ```
#
# where $k=\omega/c$ is the wave number in free space and $f(x)$ is the source term (which corresponds to a magnetic current density in Maxwell's equations).
#
# In order to simulate this scattering problem in a finite computation cell, we need outgoing (radiation) boundary conditions such that all waves at the boundary would not be reflected back since we are simulating an infinite space.
# One commonly used technique to simulate such infinite space is through the so called "perfectly matched layers" (PML) [2]. Actually, PML is not a boundary condition but an artificial absorbing "layer" that absorbs waves with minimal reflections (going to zero as the resolution increases).
# There are many formulations of PML. Here, we use one of the most flexible formulations, the "stretched-coordinate" formulation, which takes the following replace in the PDE [3]:
#
# ```math
# \frac{\partial}{\partial x}\rightarrow \frac{1}{1+\mathrm{i}\sigma(u_x)/\omega}\frac{\partial}{\partial x},
# ```
#
# ```math
# \frac{\partial}{\partial y}\rightarrow \frac{1}{1+\mathrm{i}\sigma(u_y)/\omega}\frac{\partial}{\partial y},
# ```
#
# where $u_{x/y}$ is the depth into the PML, $\sigma$ is a profile function (here we chose $\sigma(u)=\sigma_0(u/d_{pml})^2$) and different derivative corresponds to different absorption directions.
# Note that at a finite mesh resolution, PML reflects some waves, and the standard technique to mitigate this is to "turn on" the PML absorption gradually—in this case we use a quadratic profile. The amplitude $\sigma_0$ is chosen so that in the limit of infinite resolution the "round-trip" normal-incidence is some small number.
#
# Since PML absorbs all waves in $x/y$ direction, the associated boundary condition is then usually the zero Dirichlet boundary condition. Here, the boundary conditions are zero Dirichlet boundary on the top and bottom side $\Gamma_D$ but periodic boundary condition on the left ($\Gamma_L$) and right side ($\Gamma_R$).
# The reason that we use a periodic boundary condition for the left and right side instead of zero Dirichlet boundary condition is that we want to simulate a plane wave exicitation, which then requires a periodic boundary condition.
#
# Consider $\mu(x)=1$ (which is mostly the case in electromagnetic problems) and denote $\Lambda=\operatorname{diagm}(\Lambda_x,\Lambda_y)$ where $\Lambda_{x/y}=\frac{1}{1+\mathrm{i}\sigma(u_{x/y})/\omega}$, we can formulate the problem as
#
# ```math
# \left\{ \begin{aligned}
# \left[-\Lambda\nabla\cdot\frac{1}{\varepsilon(x)}\Lambda\nabla -k^2\right] H &= f(x) & \text{ in } \Omega,\\
# H&=0 & \text{ on } \Gamma_D,\\
# H|_{\Gamma_L}&=H|_{\Gamma_R},&\\
# \end{aligned}\right.
# ```
# For convenience, in the weak form and Julia implementation below we represent $\Lambda$ as a vector instead of a diagonal $2 \times 2$ matrix, in which case $\Lambda\nabla$ becomes the elementwise product.
#
# ## Numerical scheme
#
# Similar to the previous tutorials, we need to construct the weak form for the above PDEs. After integral by part and removing the zero boundary integral term, we get:
#
# ```math
# a(u,v) = \int_\Omega \left[\nabla(\Lambda v)\cdot\frac{1}{\varepsilon(x)}\Lambda\nabla u-k^2uv\right]\mathrm{d}\Omega
# ```
#
# ```math
# b(v) = \int_\Omega vf\mathrm{d}\Omega
# ```
# Notice that the $\nabla(\Lambda v)$ is also a element-wise "product" of two vectors $\nabla$ and $\Lambda v$.
#
# ## Setup
#
# We import the packages that will be used, define the geometry and physics parameters.
#

using Gridap
using GridapGmsh
using Gridap.Fields
using Gridap.Geometry

λ = 1.0          # Wavelength (arbitrary unit)
L = 4.0          # Width of the area
H = 6.0          # Height of the area
xc = [0 -1.0]    # Center of the cylinder
r = 1.0          # Radius of the cylinder
d_pml = 0.8      # Thickness of the PML
k = 2*π/λ        # Wave number
const ϵ₁ = 3.0   # Relative electric permittivity for cylinder

# ## Discrete Model
#
# We import the model from the `geometry.msh` mesh file using the `GmshDiscreteModel` function defined in `GridapGmsh`. The mesh file is created with GMSH in Julia (see the file ../assets/emscatter/MeshGenerator.jl). Note that this mesh file already contains periodic boundary information for the left and right side, and that is enough for gridap to realize a periodic boundary condition should be implemented.
#

model = GmshDiscreteModel("../models/geometry.msh")

# ## FE spaces
#
# We use the first-order lagrangian as the finite element function space basis. The dirihlet edges are labeld with `DirichletEdges` in the mesh file. Since our problem involves complex numbers (because of PML), we need to assign the `vector_type` to be `Vector{ComplexF64}`.
#

# ### Test and trial finite element function space
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe,dirichlet_tags="DirichletEdges",vector_type=Vector{ComplexF64})
U = V # mathematically equivalent to TrialFESpace(V,0)

# ## Numerical integration
#
# We generate the triangulation and a second-order Gaussian quadrature for the numerial integration. Note that we create a boundary triangulation from a `Source` tag for the line excitation. Generally, we do not need such additional mesh tags for the source, we can use a delta function to approximate such line source excitation. However, by generating a line mesh, we can increase the accuracy of this source excitation.
#

# ### Generate triangulation and quadrature from model
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
# ### Source triangulation
Γ = BoundaryTriangulation(model;tags="Source")
dΓ = Measure(Γ,degree)

# ## PML formulation
#
# Here we first define a `s_PML` function: $s(x)=1+\mathrm{i}\sigma(u)/\omega,$ and its derivative `ds_PML`. The parameter `LH` indicates the size of the inner boundary of the PML regions. Finally, we create a function-like object `Λ` that returns the PML factors and define its derivative in gridap.
# Note that here we are defining a "callable object" of type `Λ` that encapsulates all of the PML parameters. This is convenient, both because we can pass lots of parameters around easily and also because we can define additional methods on `Λ`, e.g. to express the `∇(Λv)` operation.
#

# ### PML parameters
Rpml = 1e-12      # Tolerence for PML reflection
σ = -3/4*log(Rpml)/d_pml # σ_0
LH = [L,H] # Size of the PML inner boundary (a rectangular centere at (0,0))

# ### PML coordinate streching functions
function s_PML(x,σ,k,LH,d_pml)
    u = abs.(Tuple(x)).-LH./2  # get the depth into PML
    return @. ifelse(u > 0,  1+(1im*σ/k)*(u/d_pml)^2, $(1.0+0im))
end

function ds_PML(x,σ,k,LH,d_pml)
    u = abs.(Tuple(x)).-LH./2 # get the depth into PML
    ds = @. ifelse(u > 0, (2im*σ/k)*(1/d_pml)^2*u, $(0.0+0im))
    return ds.*sign.(Tuple(x))
end

struct Λ<:Function
    σ::Float64
    k::Float64
    LH::Vector{Float64}
    d_pml::Float64
end

function (Λf::Λ)(x)
    s_x,s_y = s_PML(x,Λf.σ,Λf.k,Λf.LH,Λf.d_pml)
    return VectorValue(1/s_x,1/s_y)
end

Fields.∇(Λf::Λ) = x->TensorValue{2,2,ComplexF64}(-(Λf(x)[1])^2*ds_PML(x,Λf.σ,Λf.k,Λf.LH,Λf.d_pml)[1],0,0,-(Λf(x)[2])^2*ds_PML(x,Λf.σ,Λf.k,Λf.LH,Λf.d_pml)[2])

# ## Weak form
#
# In the mesh file, we labeled the cylinder region with `Cylinder` to distinguish it from other regions. Using this tag, we can assign material properties correspondingly (basically a function with different value in different regions). The weak form is very similar to its mathematical form in gridap.
#

# ### Intermediate variables
labels = get_face_labeling(model)
dimension = num_cell_dims(model)
tags = get_face_tag(labels,dimension)
const cylinder_tag = get_tag_from_name(labels,"Cylinder")

function ξ(tag)
    if tag == cylinder_tag
        return 1/ϵ₁
    else
        return 1.0
    end
end

τ = CellField(tags,Ω)
Λf = Λ(σ,k,LH,d_pml)

# ### Bi-linear term (from weak form)
# Note that we use a element-wise product .* here for the vector-vector product $\Lambda \nabla$
a(u,v) = ∫(  (∇.*(Λf*v))⊙((ξ∘τ)*(Λf.*∇(u))) - (k^2*(v*u))  )dΩ

# ### Source term (uniform line source)
b(v) = ∫(v)*dΓ

# ## Solver phase
#
# We can assemble the finite element operator in Gridap with the bi-linear and linear form, then solve for the field.
#

op = AffineFEOperator(a,b,U,V)
uh = solve(op)

# ## Analytical solution
# ### Theoretical analysis
# In this section, we construct the semi-analytical solution to this scattering problem, for comparison to the numerical solution. This is possible because of the symmetry of the cylinder, which allows us to expand the solutions of the Helmoltz equation in Bessel functions and match boundary conditions at the cylinder interface. (In 3d, the analogous process with spherical harmonics is known as "Mie scattering".) For more information on this technique, see Ref [4].
# In 2D cylinder coordinates, we can expand the plane wave in terms of Bessel functions (this is the Jacobi–Anger identity [5]):
#
# ```math
# H_0=\sum_m i^mJ_m(kr)e^{im\theta},
# ```
#
# where $m=0,\pm 1,\pm 2,\dots$ and $J_m(z)$ is the Bessel function of the fisrt kind.
#
# For simplicity, we start with only the $m$-th component and take it as the incident part:
#
# ```math
# H_{inc}=J_m(kr).
# ```
#
# For the scattered field, since the scattered wave should be going out, we can then expand it in terms of the Hankel function of the first kind (outgoing and incoming cylindrical waves are Hankel functions of the first and second kind [6]):
#
# ```math
# H_1=\alpha_mH_m^1(kr).
# ```
#
# For the fields inside the cylinder, we require the field to be finite at $r=0$, which then constrains the field to be only the expansion of the Bessel functions:
#
# ```math
# H_2=\beta_mJ_m(nkr),
# ```
#
# where $n=\sqrt{\varepsilon}$ is the refractive index.
#
# Applying the boundary conditions (tangential part of the electric and magnetic field to be continuous):
#
# ```math
# H_{inc}+H_1=H_2|_{r=R},
# ```
#
# ```math
# \frac{\partial H_{inc}}{\partial r}+\frac{\partial H_1}{\partial r}=\frac{1}{\epsilon}\frac{\partial H_2}{\partial r}|_{r=R}.
# ```
#
# After some math, we get:
#
# ```math
# \alpha_m=\frac{J_m(nkR)J_m(kR)^\prime-\frac{1}{n}J_m(kR)J_m(nkR)^\prime}{\frac{1}{n}H_m^1(kR)J_m(nkr)^\prime-J_m(nkr)H_m^1(kr)^\prime},
# ```
#
# ```math
# \beta_m = \frac{H_m^1(kR)J_m(kR)^\prime-J_m(kR)H_m^1(kR)^\prime}{\frac{1}{n}J_m(nkR)^\prime H_m^1(kR)-J_m(nkR)H_m^1(kR)^\prime},
# ```
#
# where $^\prime$ denotes the derivative, and the derivatives of the Bessel functions are obtained with the recurrent relations:
#
# ```math
# Y_m(z)^\prime=\frac{Y_{m-1}(z)-Y_{m+1}(z)}{2}
# ```
#
# where $Y_m$ denotes any Bessel functions (Hankel functions).
#
#
# Finally, the analytical field is ($1/2k$ is the amplitude that comes from the unit line source excitation):
# ```math
# H(r>R)=\frac{1}{2k}\sum_m\left[\alpha_mi^mH_m^1(kr)+J_m(kr)\right]e^{im\theta}
# ```
#
# ```math
# H(r\leq R)=\frac{1}{2k}\sum_m\beta_mi^mJ_m(nkr)e^{im\theta}
# ```
#

# ### Define the analytical functions
using SpecialFunctions
dbesselj(m,z) = (besselj(m-1,z)-besselj(m+1,z))/2
dhankelh1(m,z)= (hankelh1(m-1,z)-hankelh1(m+1,z))/2
α(m,n,z) = (besselj(m,n*z)*dbesselj(m,z)-1/n*besselj(m,z)*dbesselj(m,n*z))/(1/n*hankelh1(m,z)*dbesselj(m,n*z)-besselj(m,n*z)*dhankelh1(m,z))
β(m,n,z) = (hankelh1(m,z)*dbesselj(m,z)-besselj(m,z)*dhankelh1(m,z))/(1/n*dbesselj(m,n*z)*hankelh1(m,z)-besselj(m,n*z)*dhankelh1(m,z))

function H_t(x,xc,r,ϵ,λ)
    n = √ϵ
    k = 2*π/λ
    θ = angle(x[1]-xc[1]+1im*(x[2]-xc[2]))+π
    M = 40 # Number of Bessel function basis used
    H0 = 0
    if norm([x[1]-xc[1],x[2]-xc[2]])<=r
        for m=-M:M
            H0 += β(m,n,k*r)*cis(m*θ)*besselj(m,n*k*norm([x[1]-xc[1],x[2]-xc[2]]))
        end
    else
        for m=-M:M
            H0 += α(m,n,k*r)*cis(m*θ)*hankelh1(m,k*norm([x[1]-xc[1],x[2]-xc[2]]))+cis(m*θ)*besselj(m,k*norm([x[1]-xc[1],x[2]-xc[2]]))
        end
    end
    return 1im/(2*k)*H0
end
#
# ### Construct the analytical solution in finite element basis
uh_t = CellField(x->H_t(x,xc,r,ϵ₁,λ),Ω)

# ## Output and compare results
#
# The simulated field is shown below. We can see that the simulated fields and the analytical solution matched closed except for the top and PML regions. This is because the simulated source generate plane waves in two directions but we only consider the downward propagating wave in the analytical solution and the PML effect is also not considered in the analytical solution. Therefore, we just need to focus on the "center" regions which excludes the PML and top region above the source, the difference is within 6% of the field amplitude integral. As we increase the resolution, this difference should decrease (until it becomes limited by the PML reflection coefficient from $\sigma_0$, the number of Bessel function basis $M$ or by floating-point error.)
# ![](../assets/emscatter/Results.png)

# ### Save to file and view
writevtk(Ω,"demo",cellfields=["Real"=>real(uh),
        "Imag"=>imag(uh),
        "Norm"=>abs2(uh),
        "Real_t"=>real(uh_t),
        "Imag_t"=>imag(uh_t),
        "Norm_t"=>abs2(uh_t),
        "Difference"=>abs(uh_t-uh)])
# ### Compare the difference in the "center" region
function AnalyticalBox(x) # Get the "center" region
    if abs(x[1])<L/2 && abs(x[2]+0.5)<2.5
        return 1
    else
        return 0
    end
end

Difference=sqrt(sum(∫(abs2(uh_t-uh)*AnalyticalBox)*dΩ)/sum(∫(abs2(uh_t)*AnalyticalBox)*dΩ))

@assert Difference < 0.1
#
# ## References
# [1] [Wikipedia: Electromagnetic wave equation](https://en.wikipedia.org/wiki/Electromagnetic_wave_equation)
#
# [2] [Wikipedia: Perfectly matched layer](https://en.wikipedia.org/wiki/Perfectly_matched_layer)
#
# [3] A. Oskooi and S. G. Johnson, “[Distinguishing correct from incorrect PML proposals and a corrected unsplit PML for anisotropic, dispersive media](http://math.mit.edu/~stevenj/papers/OskooiJo11.pdf),” Journal of Computational Physics, vol. 230, pp. 2369–2377, April 2011.
#
# [4] Stratton, J. A. (1941). Electromagnetic Theory. New York: McGraw-Hill.
#
# [5] [Wikipedia: Jacobi–Anger expansion](https://en.wikipedia.org/wiki/Jacobi%E2%80%93Anger_expansion)
#
# [6] [Wikipedia: Bessel function](https://en.wikipedia.org/wiki/Bessel_function)


text:
using LinearAlgebra, PyPlot
import HarmonicBalance: transform_solutions, plot
export transform_solutions, plot

"""

    ODEProblem(
            eom::HarmonicEquation;
            fixed_parameters,
            x0::Vector,
            sweep::ParameterSweep,
            timespan::Tuple
            )

Creates an ODEProblem object used by DifferentialEquations.jl from the equations in `eom` to simulate time-evolution within `timespan`.
`fixed_parameters` must be a dictionary mapping parameters+variables to numbers (possible to use a solution index, e.g. solutions[x][y] for branch y of solution x).
If `x0` is specified, it is used as an initial condition; otherwise the values from `fixed_parameters` are used.
"""
function ODEProblem(eom::HarmonicEquation, fixed_parameters; sweep::ParameterSweep=ParameterSweep(), x0::Vector=[], timespan::Tuple, perturb_initial=0.0, kwargs...)

    if !is_rearranged(eom) # check if time-derivatives of the variable are on the right hand side
        eom = HarmonicBalance.rearrange_standard(eom)
    end

    fixed = Dict(fixed_parameters)
    # substitute fixed parameters
    fixed = HarmonicBalance.filter_duplicate_parameters(sweep, fixed)
    p_values = [fixed[p] for p in keys(fixed)]
    subeqs = substitute_all(Num.(getfield.(eom.equations, :lhs)), Dict(zip(keys(fixed), p_values)))
    vars = get_variables(eom)

    # substitute the harmonic variables
    eqs(v) = [substitute(eq, Dict(zip(vars, v))) for eq in subeqs]
     # substitute  sweep parameters
    eqs(v,T) = [substitute(eq,Dict(zip(keys(sweep), [sweep[p](T) for p in keys(sweep)]))) for eq in eqs(v)]

    function f!(du,u,p,T)
        state = eqs(u,T)
        for j in 1:length(vars)
            du[j] = state[j]
        end
    end

    # the initial condition is x0 if specified, taken from fixed_parameters otherwise
    initial = isempty(x0) ? real.(collect(values(fixed_parameters))[1:length(vars)]) * (1-perturb_initial) : x0

    return DifferentialEquations.ODEProblem(f!, initial, timespan; kwargs...)
end

"""
$(TYPEDSIGNATURES)

Numerically investigate the stability of a solution `soln` of `eom` within `timespan`.
The initial condition is displaced by `perturb_initial`.

Return `true` the solution evolves within `tol` of the initial value (interpreted as stable).

"""
function is_stable(soln::StateDict, eom::HarmonicEquation; timespan, tol=1E-1, perturb_initial=1E-3)
    problem = ODEProblem(eom, steady_solution=soln, timespan=timespan)
    solution = solve(problem)
    dist = norm(solution[end] - solution[1]) / (norm(solution[end]) + norm(solution[1]))
    !is_real(solution[end]) || !is_real(solution[1]) ? error("the solution is complex!") : dist < tol
end

transform_solutions(soln::OrdinaryDiffEq.ODECompositeSolution, f::String, harm_eq::HarmonicEquation) = transform_solutions(soln.u, f, harm_eq)
transform_solutions(s::OrdinaryDiffEq.ODECompositeSolution, funcs::Vector{String}, he::HarmonicEquation) = [transform_solutions(s, f, he) for f in funcs]

function plot(soln::OrdinaryDiffEq.ODECompositeSolution, funcs, harm_eq::HarmonicEquation)
    if funcs isa String || length(funcs) == 1
        plot(soln.t, transform_solutions(soln, funcs, harm_eq))
    elseif length(funcs) == 2
        plot(transform_solutions(soln, funcs, harm_eq)...)
    else
        error("Invalid plotting argument: ", funcs)
    end
end

ext: jl


text:
"""
  τ2λ(τ::Vector{Real}, λ::Vector{Real})

Suplement the vector of λ patrameters of Marshall-Olkin copula, given some of those
parameters and a vector fo Kendall's τ correlations. Wroks fo 2,3 variate MO copula
"""
function τ2λ(τ::Vector{T}, λ::Vector{T}) where T <: Real
  if length(τ) == 1
    return [λ[1],λ[2],(λ[1]+λ[2])*τ[1]/(1-τ[1])]
  else
    t = τ./(1 .-τ)
    M = Matrix(1.0I, 3, 3) .-(ones(3,3) .-Matrix(1.0I, 3, 3)) .*t
    fm = hcat([1. 1. 0.; 0. 1. 1.; 1. 0. 1.].*t, .-[1.; 1.; 1.])
    ret = map(x -> (x > 0) ? (x) : (0.000000001), inv(M)*fm*λ)
    [λ[1:3]..., ret..., λ[end]]
  end
end

function moρ2τ(ρ::T) where T <: Real
  f(τ) = 1/2 .*sin.(τ[1]*pi/2)+τ[1]/2 - ρ[1]
  fzero(f, -0.999, 0.999)
end


text:
using PyCall
using Random

export create_qft_circuit, create_simple_preparation_circuit
export create_ghz_preparation_circuit
export create_RQC

"""
    create_qft_circuit(n::Int, m::Int=-1)

Generate QFT circuit acting on n qubits with approximation parameter m. To no approximation
m is set to -1 which is the default value.
"""
function create_qft_circuit(n::Int, m::Int=-1)
    qiskit = pyimport("qiskit")
    circ = qiskit.QuantumCircuit(n)
    m = (m == -1) ? n : m # if m is -1 then do not truncate
    for i in 1:n
        circ.h(i-1)
        for j in i:n-1
            if j - i < m
                circ.cu1(π/(2^(j-i+1)), j, i-1)
            end
        end
        circ.barrier()
    end
    # TODO: do not add swap gates but return ordering information instead
    for i = 1:convert(Int, floor(n//2))
        circ.swap(i - 1, n - i)
    end
    circ
end

"""
    create_simple_preparation_circuit(qubits::Int,
                                               depth::Int,
                                               seed::Int=nothing)

Create a simple preparation circuit which mixes
"""
function create_simple_preparation_circuit(qubits::Int,
                                           depth::Int,
                                           seed::Union{Int, Nothing}=nothing)
    if seed === nothing
        rng = MersenneTwister()
    else
        rng = MersenneTwister(seed)
    end

    qiskit = pyimport("qiskit")
    circ = qiskit.QuantumCircuit(qubits)

    # add layer of Hadamard gates
    for q in 1:qubits
        circ.h(q - 1)
    end
    circ.barrier()

    # add layers of random unitaries and contolled not gates
    for d in 1:depth
        params = rand(rng, qubits, 3)
        for q in 1:qubits
            circ.u3(params[q,:]..., q - 1)
        end
        start_qubit = ((d + 1) % 2) + 1
        for q in start_qubit:2:qubits - 1
            circ.cx(q - 1, q)
        end
        circ.barrier()
    end
    circ
end

"""
    create_ghz_preparation_circuit(qubits::Int)

Create a ghz preparation circuit
"""
function create_ghz_preparation_circuit(qubits::Int)
    qiskit = pyimport("qiskit")
    circ = qiskit.QuantumCircuit(qubits)

    circ.h(0)
    for q in 1:qubits-1
        circ.cx(q-1, q)
    end

    circ
end

# *************************************************************************** #
#                              RQC functions
# *************************************************************************** #

"Struct to represent a random quantum circuit"
struct RQC
    # The dimensions of the 2-dimensional array of qubits.
    n::Int
    m::Int

    # A qiskit circuit object for the R.Q.C.
    circ

    # A 2-dimensional array of integers to keep track of which single qubit
    # gate should next be applied to a given qubit in the circuit.
    # The integers map to gates as follows:
    # 1 => T, 2 => sqrt(X), 3 => sqrt(Y).
    # The sign of the integer is used to indicate if the corresponding qubit
    # has been hit by a two qubit gate and can be hit by a single qubit gate.
    # This should be initialised as an array of ones as the first single qubit
    # gate is always a T gate.
    next_gate::Array{Int, 2}

    # A dictionary to convert integers from the above array into gate names.
    single_qubit_gates::Dict{Int, String}
end

function RQC(n::Int, m::Int)
    qiskit = pyimport("qiskit")
    RQC(n, m, qiskit.QuantumCircuit(n*m), -ones(Int, n,m),
        Dict(1=>"t", 2=>"x", 3=>"y"))
end

"""
    random_gate!(rqc::RQC, i::Int, j::Int)

Returns the gate name for the next single qubit gate to be applied to the qubit
at (i, j).
"""
function random_gate!(rqc::RQC, i::Int, j::Int, rng::MersenneTwister)
    gate = rqc.next_gate[i, j]
    if gate > 0
        rqc.next_gate[i, j] = (rqc.next_gate[i, j] + rand(rng, [0,1]))%3 + 1
        rqc.next_gate[i, j] *= -1
        return rqc.single_qubit_gates[gate]
    else
        return ""
    end
end

"""
    patterns(rqc::RQC)

Generate a dictionary of patterns of qubit pairs to apply two qubit gates to.
The patterns are numbers 1 to 8 like in Boxio_2018 but in a different order.
(Note, the order used in Boxio_2018 is 3, 1, 6, 8, 5, 7, 2, 4)
"""
function patterns(rqc::RQC)
    key = 1
    gp = Dict()
    N = max(rqc.m, rqc.n)
    for v in [1, -1]
        for shift in [0, 1]
            for flip in [false, true]
                gp[key] = [[[i, j], [i, j+v]] for j in 1:2:N
                                              for i in 1+((j+shift)÷2)%2:2:N]
                if flip; gp[key] = map(x->reverse.(x), gp[key]); end
                filter!(gp[key]) do p
                    out_of_bounds(rqc, p)
                end
                key += 1
            end
        end
    end
    gp
end

"""
    out_of_bounds(rqc::RQC, targets::Array{Array{Int, 1}, 1})

Check if any of the qubits, in the given qubit pair, are out of bounds.
"""
function out_of_bounds(rqc::RQC, targets::Array{Array{Int, 1}, 1})
    ((i, j), (u, v)) = targets
    (0 < i <= rqc.n) && (0 < u <= rqc.n) &&
    (0 < j <= rqc.m) && (0 < v <= rqc.m)
end

"""
    create_RQC(rows::Int, cols::Int, depth::Int,
                        seed::Union{Int, Nothing}=nothing;
                        use_iswap::Bool=false,
                        final_Hadamard_layer::Bool=false)

Generate a random quantum circuit with a (rows x cols) grid of qubits. The
parameter 'depth' is the number of one and two qubit gate layers to be used in
the circuit.

The single qubit gates used are chossen randomly from the set
{T, sqrt(X), sqrt(Y)}. The seed for the random selection can be set with the
seed parameter.

CZ gates are used by default as the two qubit gates. Setting 'use_iswap' to
true will use iSWAP gates inplace of CZ gates.

Setting 'final_Hadamard_layer' to true will include a layer of Hadamard gates
at the end of the circuit.
"""
function create_RQC(rows::Int, cols::Int, depth::Int,
                    seed::Union{Int, Nothing}=nothing;
                    use_iswap::Bool=false,
                    final_Hadamard_layer::Bool=false)

    if seed === nothing
        rng = MersenneTwister()
    else
        rng = MersenneTwister(seed)
    end

    # Create an empty rqc.
    rqc = RQC(rows, cols)

    # A dictionary to map gate names to qiskit circuit methods that add the
    # corresponding gate to the rqc.
    # TODO: x and y should be replaced by sqrt of x and y
    gates = Dict("h" => rqc.circ.h, "cz" => rqc.circ.cz, "t" => rqc.circ.t,
                 "x" => target -> rqc.circ.rx(pi/2, target),
                 "y" => target -> rqc.circ.ry(pi/2, target),
                 "iswap" => rqc.circ.iswap)

    # A function to add a gate to the rqc with specific target qubits.
    function add_gate!(rqc::RQC, gate::String, target_qubits)
        if !(length(gate) == 0)
            target_qubits = [i + (j-1)*rqc.n - 1 for (i,j) in target_qubits]
            gates[gate](target_qubits...)
        end
    end

    # The rqc starts with a layer of hadamard gates applied to all qubits.
    for i in 1:rqc.n, j in 1:rqc.m
        add_gate!(rqc, "h", [[i, j]])
    end

    # Get the patterns to apply two qubit gates in.
    gate_patterns = patterns(rqc)

    # Loop through the patterns applying the two qubit gates and apply a random
    # single qubit gate to the appropriate qubits.
    pattern = [3, 1, 6, 8, 5, 7, 2, 4]
    two_qubit_gate = use_iswap ? "iswap" : "cz"
    for i in 0:depth-1
        i = pattern[i%8 + 1]

        if !(length(gate_patterns[i]) == 0)
            # Apply the two qubit gates.
            for targets in gate_patterns[i]
                add_gate!(rqc, two_qubit_gate, targets)
            end

            # Apply the single qubit gates
            qubits_hit = reduce(vcat, gate_patterns[i])
            for i in 1:rqc.n, j in 1:rqc.m
                if !([i, j] in qubits_hit)
                    gate = random_gate!(rqc, i, j, rng)
                    add_gate!(rqc, gate, [[i, j]])
                end
            end

            # Update rqc with which qubits were hit by a two qubit gate
            # and so may be hit by a single qubit gate next round.
            for ((i, j), (u, v)) in gate_patterns[i]
                rqc.next_gate[i, j] = abs(rqc.next_gate[i, j])
                rqc.next_gate[u, v] = abs(rqc.next_gate[u, v])
            end
        end

        rqc.circ.barrier()
    end

    if final_Hadamard_layer
        for i in 1:rqc.n, j in 1:rqc.m
            add_gate!(rqc, "h", [[i, j]])
        end
    end

    rqc.circ
end


text:
# # Classification of MNIST digits with a convolutional network.
# This program writes out saved model to the file "mnist_conv.bson".
# This file demonstrates basic implimentation of ~
# - data structure
# - construction
# - pooling
# - training
# - saving
# - conditional early-exit
# - learning rate scheduling.
#
#
# **This model, while simple, should hit around 99% test accuracy after training for approximately 20 epochs.**
#
# For core concepts of ML and Convolution check out [towardsdatascience's guide](https://towardsdatascience.com/a-comprehensive-guide-to-convolutional-neural-networks-the-eli5-way-3bd2b1164a53) for a detailed explanation of everything
#
# ## Importing required libraries
# Can't do machine learning in flux whithout `Flux` :)
using Flux
# Import the mnist dataset
using Flux.Data.MNIST
# `Statistics` to calculate the mean which is required for finding the accuracy
using Statistics
# The uses of these will be explained later
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
# `Printf` For outputing the accuracy and other information and `BSON` forsaving the model
using Printf, BSON

# ## Load labels and images from Flux.Data.MNIST
@info("Loading data set")
train_labels = MNIST.labels()
train_imgs = MNIST.images()

# ## Reshaping and batching
# The `train_imgs` is **not** of the right shape. So, it can't be directly divided into smaller parts for batching.
# For this, we are creating another array of the correct shape which will be our batch, and we'll fill this array with the data from our `train_imgs`.
# The correct shape, which Flux expects the batch to be in, is **(IMG_DATA, IMG_DATA, COLOR_CHANNEL, SAMPLES)**.
#
# A good way to think about this is that every image is a *2D* rectangle. Which has a hight and width of 28.
# Flux requires us to add another dim which represents the color. So a ***BLACK AND WHITE*** image would be like a cuboid of thickness 1, hight and width of 28.
# An ***RGB*** image would be like cuboid thickness 3, made up of 3 rectangles stackted on top of each other . each rectangle representing red, blue and green channels.
#
# This type of thinking is good for beginners as it helps visualize things.
#
# Note that we're adding image to the last (samples) dim
# So, a batch can be visualized as multiple cubiods, each of which is made up of *3* rectangles (R - rectangle, G - rectangle,)
function make_minibatch(X, Y, idxs)
    X_batch = Array{Float32}(undef, size(X[1])..., 1, length(idxs))
    for i in 1:length(idxs)

        X_batch[:, :, :, i] = Float32.(X[idxs[i]])
    end
    Y_batch = onehotbatch(Y[idxs], 0:9)
    return (X_batch, Y_batch)
end

# In the above, onehotencoding turns neumerical data into kind of truth tabels.
# example,`onehot(:b, [:a, :b, :c])` whill output
# ```julia
# 3-element Flux.OneHotVector:
# false
#  true
# false
# ```
# onehotbatch just one hot encodes all the data in a provided array
# and yields the encoding of every element in an output *OneHotMatrix*.

batch_size = 128
mb_idxs = partition(1:length(train_imgs), batch_size)
train_set = [make_minibatch(train_imgs, train_labels, i) for i in mb_idxs]

# We are dividing an array of numbers of the size of the imgs into the sizes of batches and saving it in `mb_idxs`.
# This **array of arrays** of numbers will be used as indexs.

# ## Prepare test set as one giant minibatch:
test_imgs = MNIST.images(:test)
test_labels = MNIST.labels(:test)
test_set = make_minibatch(test_imgs, test_labels, 1:length(test_imgs))

# ## Define our model.
# We will use a simple convolutional architecture with
# three iterations of Conv -> ReLU -> MaxPool, followed by a final Dense
# layer that feeds into a softmax probability output.
#
@info("Constructing model...")
model = Chain(

    Conv((3, 3), 1=>16, pad=(1,1), relu),
    MaxPool((2,2)),

    Conv((3, 3), 16=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    Conv((3, 3), 32=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10),

    softmax,
)

# A conv layer is defined as `Conv((3, 3), 1=>16, pad=(1,1), relu)`
#
# `(3,3)` is the size of the filter that will be convolving.
# `1=>16` follows the `input_size=>output_size` format.
# `relu` is the name of the activation function we're gonna be using
#
# First convolution, operating upon a 28x28 image
# Second convolution, operating upon a 14x14 image
# Third convolution, operating upon a 7x7 image
# Then reshape 3d tensor into a 2d one, at this point it should be (3, 3, 32, N)
# which is where we get the 288 in the `Dense`.
#
# Finally, softmax to get nice probabilities

# ## Generic Pre-processing
# Load model and datasets onto GPU, if enabled
train_set = gpu.(train_set)
test_set = gpu.(test_set)
model = gpu(model)

# Make sure our model is nicely precompiled before starting our training loop
model(train_set[1][1])

# `loss()` calculates the crossentropy loss between our prediction `y_hat`
# (calculated from `model(x)`) and the ground truth `y`.  We augment the data
# a bit, adding gaussian random noise to our image to make it more robust and avoid overfitting.
#
function loss(x, y)

    x_aug = x .+ 0.1f0*gpu(randn(eltype(x), size(x)))

    y_hat = model(x_aug)
    return crossentropy(y_hat, y)
end
accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
# Onecold does what is sounds like it might do, it turns one hotencoded vectors into the orignal data.
#
# There are two reasons for using onecold.
# One is that m(x) has its last layer as `softmax`. This means it returns a probability distribution.
# We can not equate PB distributions. So, we need to turn the probablity distribution into neumerical data.
# As it turns out, onecold works with probablity distributions too.
# This means it turns the PB into normal neumerical data which can be used with the `==` sign.

opt = ADAM(0.001)
# Train our model with the given training set using the ADAM optimizer and
# printing out performance against the test set as we go.

# ## Training
@info("Beginning training loop...")
best_acc = 0.0
last_improvement = 0
for epoch_idx in 1:100
    global best_acc, last_improvement
    # Train for a single epoch
    Flux.train!(loss, params(model), train_set, opt)

    # Calculate accuracy:
    acc = accuracy(test_set...)
    @info(@sprintf("[%d]: Test accuracy: %.4f", epoch_idx, acc))

    # If our accuracy is good enough, quit out.
    if acc >= 0.999
        @info(" -> Early-exiting: We reached our target accuracy of 99.9%")
        break
    end

    # If this is the best accuracy we've seen so far, save the model out
    if acc >= best_acc
        @info(" -> New best accuracy! Saving model out to mnist_conv.bson")
        BSON.@save joinpath(dirname(@__FILE__), "mnist_conv.bson") model epoch_idx acc
        best_acc = acc
        last_improvement = epoch_idx
    end

    # If we haven't seen improvement in 5 epochs, drop our learning rate:
    if epoch_idx - last_improvement >= 5 && opt.eta > 1e-6
        opt.eta /= 10.0
        @warn(" -> Haven't improved in a while, dropping learning rate to $(opt.eta)!")

        # After dropping learning rate, give it a few epochs to improve
        last_improvement = epoch_idx
    end

    if epoch_idx - last_improvement >= 10
        @warn(" -> We're calling this converged.")
        break
    end
end

max_stars_repo_name: Akshat-mehrotra/model-zoo

max_issues_repo_name: Akshat-mehrotra/model-zoo

max_forks_repo_name: Akshat-mehrotra/model-zoo


text:
@testset "update_y_l_parallel" begin
#test update_y_l_parallel.jl

TF=Float64
TI=Int64

#without BLAS
x=randn(100)
p=2
i=10
Blas_active = true

y=Vector{Vector{Float64}}(undef, p)
y[1]=randn(51)
y[2]=randn(100)
y_old=Vector{Vector{Float64}}(undef, p)
y_old[1]=randn(51)
y_old[2]=randn(100)
l_old=Vector{Vector{Float64}}(undef, p)
l_old[1]=randn(51)
l_old[2]=randn(100)
l=Vector{Vector{Float64}}(undef, p)
l[1]=randn(51)
l[2]=randn(100)

rho=[1.234;10.23432]
gamma=Vector{Float64}(undef, p)
gamma[1]=1.0
gamma[2]=1.345

prox=Vector{Any}(undef, p)
prox[1] = x -> 1.0.*x
m=randn(100)
prox[2] = x -> prox_l2s!(x,rho[2],m)
TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef, p)
TD_OP[1] = sparse(1.0*I, 51,100)*2.0
TD_OP[2] = sparse(1.0*I, 100, 100)
maxit=39
log_PSDMM = log_type_PARSDMM(zeros(maxit,p-1),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());
P_sub = Vector{Any}(undef, 1)
P_sub[1] = prox[1]
counter=12
x_hat=Vector{Vector{Float64}}(undef, p)
x_hat[1]=randn(51)
x_hat[2]=randn(100)
r_pri=Vector{Vector{Float64}}(undef, p)
r_pri[1]=randn(51)
r_pri[2]=randn(100)
s=Vector{Vector{Float64}}(undef, p)
s[1]=randn(51)
s[2]=randn(100)

set_feas=zeros(length(TD_OP))

#distributed quantities
s_d      = distribute(s)
r_pri_d  = distribute(r_pri)
x_hat_d  = distribute(x_hat)
prox_d   = distribute(prox)
P_sub_d  = distribute(P_sub)
TD_OP_d  = distribute(TD_OP)
rho_d    = distribute(rho)
gamma_d  = distribute(gamma)
y_d      = distribute(y)
y_old_d  = distribute(y_old)
l_d      = distribute(l)
l_old_d  = distribute(l_old)
set_feas_d=distribute(set_feas)

#distributed computation
[ @spawnat pid  update_y_l_parallel(x,i,Blas_active,
  y_d[:L],y_old_d[:L],l_d[:L],l_old_d[:L],rho_d[:L],gamma_d[:L],prox_d[:L],TD_OP_d[:L]
  ,P_sub_d[:L],x_hat_d[:L],r_pri_d[:L],s_d[:L],set_feas_d[:L]) for pid in y_d.pids]

  rhs_next_iter=zeros(length(x))
  rhs_next_iter=@distributed (+) for i=1:2
      #rhs_d[i]
      TD_OP_d[i]'*(l_d[i].+y_d[i])
  end

# (y,l,r_pri,s,log_PSDMM,counter,y_old,l_old)=update_y_l_parallel(x,p,i,Blas_active,
#   y_d[:L],y_old_d[:L],l_d[:L],l_old_d[:L],rho_d[:L],gamma_d[:L],prox_d[:L],TD_OP_d[:L],
#   log_PSDMM,P_sub_d[:L],counter,x_hat_d[:L],r_pri_d[:L],s_d[:L])

#construct rerefrence solution
x_hat2=deepcopy(x_hat)
y2=deepcopy(y)
l2=deepcopy(l)
x2=deepcopy(x)

for i=1:p
  x_hat2[i]=gamma[i]*TD_OP[i]*x2 + (1-gamma[i])*y2[i]
  y2[i]=prox[i](x_hat2[i]-l2[i]./rho[i])
  l2[i].=l2[i].+rho[i].*(y2[i].-x_hat2[i])
end
rhs_next_iter_ref=zeros(length(x))
for i  = 1:p
    #rhs_d[i]
    rhs_next_iter_ref=rhs_next_iter_ref+TD_OP[i]'*(l2[i].+y2[i])
end

#use serial update_y_l.jl
(y,l,r_pri,s,log_PSDMM,counter,y_old,l_old)=update_y_l(x,p,i,Blas_active,
  y,y_old,l,l_old,rho,gamma,prox,TD_OP,log_PSDMM,P_sub,counter,x_hat,r_pri,s)
  rhs_next_iter_serial=zeros(length(x))
  for i  = 1:p
      #rhs_d[i]
      rhs_next_iter_serial=rhs_next_iter_serial+TD_OP[i]'*(l[i].+y[i])
  end

@test isapprox(rhs_next_iter_serial,rhs_next_iter_ref,rtol=10*eps())
@test isapprox(rhs_next_iter,rhs_next_iter_ref,rtol=10*eps())
@test isapprox(s,s_d,rtol=10*eps())
@test isapprox(l,l_d,rtol=10*eps())
@test isapprox(r_pri,r_pri_d,rtol=10*eps())

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using MonotoneCubicInterpolation
using Base.Test

# interpolate a linear function
include("interpolate_knownfun.jl")


text:
@testset "Geometry" begin
    @testset "Boosts" begin
        β0 = Vec3(0, 0, 0)
        u1 = Vec4(1.2, 0.2, 0.2, 1.)
        @test boost(u1, β0) ≈ u1
        β1 = Vec3(0, 0, 0.85)
        v1 = boost(u1, β1)
        @test u1⋅u1 ≈ v1⋅v1
        β2 = u1 / u1.t
        v2 = boost(u1, β2)
        @test v2.t^2 ≈ u1⋅u1
        ε = 1e-16
        @test norm(Vec3(v2)) < ε
        β3 = Vec3(0.55, -0.37, 0.11)
        v3 = boost(u1, β3)
        u3 = boost(v3, -β3)
        @test u1 ≈ u3
    end

    @testset "Rotations" begin
        x̂ = Vec3(1, 0, 0)
        ŷ = Vec3(0, 1, 0)
        ẑ = Vec3(0, 0, 1)
        @test rotate(x̂, ẑ, π/2) ≈ ŷ
        @test rotate(ŷ, x̂, π/2) ≈ ẑ
        @test rotate(ẑ, ŷ, π/2) ≈ x̂
        @test rotate(x̂, ẑ, π) ≈ -x̂
        @test rotate(x̂, ŷ, π) ≈ -x̂
        @test rotate(ŷ, ẑ, π) ≈ -ŷ
        @test rotate(ŷ, x̂, π) ≈ -ŷ
        @test rotate(ẑ, x̂, π) ≈ -ẑ
        @test rotate(ẑ, ŷ, π) ≈ -ẑ
        u = Vec3(1.1, -0.3, 0.9)
        k̂ = normalize(Vec3(-0.6, 0.8, 1.3))
        @test rotate(u, k̂, -π/3) ≈ rotate(u, -k̂, π/3)
        @test norm(rotate(u, k̂, 0.6)) ≈ norm(u)
        @test rotate(u, normalize(u), 1.3) ≈ u
        v̂ = normalize(k̂ × u)
        @test rotate(u, v̂, π) ≈ -u
    end
end


text:
function download()
  # File Names
  DataDir = joinpath(dirname(@__FILE__),"data")
  CompressFile = joinpath(DataDir,"bindata.tar.gz")
  CifarOrigDir = joinpath(DataDir,"cifar-10-batches-bin")
  CifarFinalDir = joinpath(DataDir,"bin")

  # Download Compressed Files
  download("http://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz",CompressFile)

  # Decompress Dataset
  cd(DataDir)
  run(`tar -zxvf $CompressFile`)
  rm(CompressFile)

  # Rename decompressed directory (for ease)
  mv(CifarOrigDir,CifarFinalDir;remove_destination=true)
end


text:
# using EffectiveWaves, LinearAlgebra
# import StaticArrays: SVector

import MultipleScattering: outgoing_translation_matrix

function discrete_system(ω::T, source::AbstractSource{Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::WithoutSymmetry{Dim}; kws...) where {T,Dim}

    return discrete_system(ω, source, material, AzimuthalSymmetry{Dim}(); kws...)
end

"""
    discrete_system(k,r1_vec::Vector,r2,θ2; basis_order::Int = 3)

documentation
"""
function discrete_system(ω::T, source::AbstractSource{Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::AbstractAzimuthalSymmetry{Dim};
        basis_order::Int = 1,
        basis_field_order::Int = 2,
        legendre_order::Int = basis_field_order + 1,
        mesh_points::Int = Int(round(sqrt(1.1 * (basis_field_order) * legendre_order ))) + 2,
        rtol::T = 1e-2,
        maxevals::Int = Int(2e4),
        numdensity = (x1, s1) -> number_density(s1),
        pair_corr = hole_correction_pair_correlation
    ) where {T,Dim}

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity(x1,s1) = scale_number_density * numdensity(x1,s1)

    R = outer_radius(material.shape)

    gs = regular_spherical_coefficients(source)(basis_field_order,origin(material.shape),ω);

    v = regular_basis_function(source.medium,  ω)

    Uout = outgoing_translation_matrix(ω, source.medium, material;
        basis_order = basis_order, tol = rtol
    )

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rθφ2xyz = radial_to_cartesian_coordinates

    r1s = LinRange(0,R-outer_radius(s1), mesh_points)
    θ1s = LinRange(0,π, mesh_points)

    ls, ms = spherical_harmonics_indices(basis_order)
    azi_inds(m::Int) = lm_to_spherical_harmonic_index.(abs(m):basis_field_order,-m)

    len = basisorder_to_basislength(Acoustic{T,Dim}, basis_order)
    len_p = sum(legendre_order for nd in 1:len for i in azi_inds(ms[nd]))

    function incident_coefficients(r1s::AbstractVector{T},θ1s::AbstractVector{T})
        lm2n = lm_to_spherical_harmonic_index

        coefs = [
            begin
                vs = v(basis_order + basis_field_order, rθφ2xyz(SVector(r1,θ1,zero(T))))
                data = [
                    t_diags[1][lm2n(l,m)] *
                    sum(
                        gaunt_coefficient(dl,0,l,m,l1,-m) * vs[lm2n(l1,-m)] * gs[lm2n(dl,0)]
                    for dl in 0:basis_field_order for l1 in max(abs(m),abs(dl-l)):(dl+l))
                for l = 0:basis_order for m = -l:l]

                data
            end
        for r1 in r1s, θ1 in θ1s][:];

        return vcat(coefs...)
    end

    function field_basis(rθφ::AbstractVector{T})
        Ys = spherical_harmonics(basis_field_order, rθφ[2], rθφ[3]);
        P = Legendre{T}()

        # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
        Prs = P[2 * rθφ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]

        return [Pr * Yθ for Pr in Prs, Yθ in Ys]
    end
    # function field_basis(rθ::AbstractVector{T})
    #     P = Legendre{T}()
    #
    #     # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
    #     Pθs = P[cos(rθ[2]), 1:legendre_order]
    #     Prs = P[2 * rθ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]
    #
    #     # [Pr * Pθ for Pr in Prs, Pθ in Pθs][:]
    #     return (Prs * transpose(Pθs))[:]
    # end

    function kernel_function(rθ1::SVector{2,T})
        x1 = rθφ2xyz(SVector(rθ1[1],rθ1[2],zero(T)))
        Kzero = zeros(Complex{T},len,len_p)

        fun = function (rθφ::SVector{3,T})
            x2 = rθφ2xyz(rθφ)
            if pair_corr(x1,s1,x2,s1) ≈ zero(T)
                return Kzero
            end
            basis2 = field_basis(rθφ)

            U = Uout(x1 - x2)
            U = U .* (bar_numdensity(x2,s1) * pair_corr(x1,s1,x2,s1) * sin(rθφ[2]) * rθφ[1]^2)

            return reshape(
                [
                    t_diags[1][n] * U[nd,n] * b2
                for nd in 1:len for b2 in basis2[:,azi_inds(ms[nd])][:] for n in 1:len],
            (len, :))
        end

        return fun
    end

    function δφj(rθφ1::AbstractVector{T})
        basis1 = field_basis(rθφ1)
        return reshape(
            [
                (nd == n) ? b1 : zero(Complex{T})
            for nd in 1:len for b1 in basis1[:,azi_inds(ms[nd])][:] for n in 1:len],
        (len, len_p))
    end

    δφj(rθ1::SVector{2,T}) = δφj(SVector(rθ1[1],rθ1[2],zero(T)))

    test_ker = kernel_function(SVector(mean(r1s),θ1s[1]))
    (I,E) = hcubature(test_ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
        rtol=rtol, maxevals=maxevals
    );

    println("The estimated max coefficient of the integrated kernel is:")
    println("I = ", maximum(abs.(I)) )
    println("with an estimated error of: ")
    println("E = ", E)

    Ks = [
        begin
            rθ1 = SVector(r1,θ1)
            ker = kernel_function(rθ1)
            ker_integrated = hcubature(ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
                rtol=rtol, maxevals=maxevals
            )[1]

            δφj(rθ1) - ker_integrated
        end
    for r1 in r1s, θ1 in θ1s][:];

    bigK = vcat(Ks...);

    bs = incident_coefficients(r1s,θ1s);

    Fs = bigK \ bs;

    ## Alternative:
    # as = inv(transpose(conj.(bigK)) * bigK) * transpose(conj.(bigK)) * bs;

    if norm(bigK * Fs - bs) / norm(bs) > rtol
        @warn "Numerical solution has a relative residual error of $(norm(bigK * Fs - bs) / norm(bs)), where the requested relative tolernance was: $rtol. This residual error can be decreased by increasing the basis_field_order (current value: $basis_field_order) for the field."
    end

    function scattered_field(xs::AbstractVector{T})
        rθφ = cartesian_to_radial_coordinates(xs)
        return δφj(rθφ) * Fs
        # The factor exp(-im * m * φ) is due to azimuthal symmetry
        # azi_factor = exp.((-im*rθφ[3]) .* ms)
        # return azi_factor .* (as * field_basis(rθφ[1:2]))
    end

    return ScatteringCoefficientsField(ω, source.medium, material, scattered_field;
        symmetry = AzimuthalSymmetry{Dim}(),
        basis_order = basis_order,
        basis_field_order = basis_field_order
    )
end

function discrete_system(ω::T, source::AbstractSource{Acoustic{T,3}}, material::Material{3,Sphere{T,3}}, ::RadialSymmetry{3};
        basis_order::Int = 1,
        basis_field_order::Int = Int(round(T(2) * real(ω / source.medium.c) * outer_radius(material.shape))) + 1,
        legendre_order::Int = 2basis_field_order + 1,
        # mesh_points::Int = 3 * (legendre_order + 1),
        mesh_points::Int = 3 * (legendre_order + 1),
        rtol::T = 1e-2, ptol::T = 1e-2,
        maxevals::Int = Int(2e4),
        numdensity = (x1, s1) -> number_density(s1),
        pair_corr = hole_correction_pair_correlation
    ) where T

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    a12 = 2.0 * s1.exclusion_distance * outer_radius(s1);

    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity(x1,s1) = scale_number_density * numdensity(x1,s1)

    R = outer_radius(material.shape)

    g0 = regular_spherical_coefficients(source)(1,origin(material.shape),ω)[1];

    v = regular_basis_function(source.medium,  ω)

    Uout = outgoing_translation_matrix(ω, source.medium, material;
        basis_order = basis_order, tol = rtol
    )

    lm_to_n = lm_to_spherical_harmonic_index
    ls_to_ns = lm_to_n.(0:basis_order,0)

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rθφ2xyz = radial_to_cartesian_coordinates

    r1s = LinRange(0,R-outer_radius(s1), mesh_points)

    len = basis_order + 1
    len_p = legendre_order * len

    function incident_coefficients(r1s::AbstractVector{T})

        coefs = [
            begin
                vs = v(basis_order, rθφ2xyz([r1,zero(T),zero(T)]))
                t_diags[1][ls_to_ns] .* vs[ls_to_ns] .* (-T(1)) .^ (0:basis_order)
            end
        for r1 in r1s][:];

        return vcat((sqrt(4pi) * g0) .* coefs...)
    end

    function field_basis(rθφ::AbstractVector{T})
        Ys = conj.(spherical_harmonics(basis_order, rθφ[2], rθφ[3]));
        P = Legendre{T}()

        Prs = P[2 * rθφ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]

        return [Pr * Yθ for Pr in Prs, Yθ in Ys]
    end

    function kernel_function(r1::T)
        x1 = rθφ2xyz([r1,zero(T),zero(T)])
        Kzero = zeros(Complex{T},len,len_p)

        fun = function (rθφ::SVector{3,T})
            x2 = rθφ2xyz(rθφ)
            if r1 + rθφ[1] < a12
                return Kzero
            end
            if pair_corr(x1,s1,x2,s1) < ptol || norm(x1 - x2) < a12
                return Kzero
            end
            basis2 = field_basis(rθφ)

            U = Uout(x1 - x2)
            U = U .* (bar_numdensity(x2,s1) * pair_corr(x1,s1,x2,s1) * sin(rθφ[2]) * rθφ[1]^2)

            data = [
                sum([
                        t_diags[1][lm_to_n(l,0)] * U[lm_to_n(dl,dm),lm_to_n(l,0)] * b2
                    for b2 in basis2[:,lm_to_n(dl,dm)][:] for l in 0:basis_order]
                for dm = -dl:dl)
            for dl in 0:basis_order];

            return reshape(vcat(data...), (len, len_p))
        end

        return fun
    end

    function δφj(rθφ1::AbstractVector{T})
        basis1 = field_basis(rθφ1)
        return reshape(
            [
                (dl == l) ? b1 : zero(Complex{T})
            for dl in 0:basis_order for b1 in basis1[:,lm_to_n(dl,0)][:] for l in 0:basis_order],
        (len, len_p))
    end

    δφj(r1::T) = δφj([r1,zero(T),zero(T)])

    test_ker = kernel_function(mean(r1s))
    (maxvalue,estimate_error) = hcubature(test_ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
        rtol=rtol, maxevals=maxevals
    );

    println("The estimated max coefficient of the integrated kernel is:")
    println(maximum(abs.(maxvalue)) )
    println("with an estimated error of: ")
    println(estimate_error)

    Ks = [
        begin
            ker = kernel_function(r1)
            ker_integrated = hcubature(ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
                rtol=rtol, maxevals=maxevals
            )[1]

            δφj(r1) - ker_integrated
        end
    for r1 in r1s];

    bigK = vcat(Ks...);

    bs = incident_coefficients(r1s);
    Fs = bigK \ bs;

    ## Alternative:
    # as = inv(transpose(conj.(bigK)) * bigK) * transpose(conj.(bigK)) * bs;

    if norm(bigK * Fs - bs) / norm(bs) > rtol
        @warn "Numerical solution has a relative residual error of $(norm(bigK * Fs - bs) / norm(bs)), where the requested relative tolernance was: $rtol. This residual error can be decreased by increasing the basis_field_order (current value: $basis_field_order) for the field."
    end

    function scattered_field(x::AbstractVector{T})
        rθφ = cartesian_to_radial_coordinates(x)
        basis1 = field_basis(rθφ)

        Fs = reshape(Fs, (legendre_order,len))

        # F = [[k,dl] for dl in 0:basis_order for k = 1:3]
        # basis1[k,(l,m)] * F[(k,l)]
        return [sum(basis1[:, lm_to_n(l,m)] .* Fs[:,l+1]) for l = 0:basis_order for m = -l:l]
    end

    return ScatteringCoefficientsField(ω, source.medium, material, scattered_field;
        symmetry = RadialSymmetry{3}(),
        basis_order = basis_order,
        basis_field_order = basis_field_order
    )
end

function discrete_system_radial(ω::T, source::AbstractSource{Acoustic{T,3}}, material::Material{3,Sphere{T,3}}, ::RadialSymmetry{3};
        basis_order::Int = 1,
        basis_field_order::Int = Int(round(T(2) * real(ω / source.medium.c) * outer_radius(material.shape))) + 1,
        # basis_field_order = Int(round(T(2) * real(ω / source.medium.c) * outer_radius(material.shape))) + 1,
        polynomial_order::Int = 10,
        legendre_order::Int = 2basis_field_order + 1,
        mesh_points::Int = 3 * (legendre_order + 1),
        # mesh_points = Int(basis_field_order + 1)^2,
        numdensity::Function = (x1, s1) -> number_density(s1),
        # numdensity = (x1, s1) -> number_density(s1),
        scheme = :trapezoidal, pair_kws...
    ) where T

    # basis_field_order = 14

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity(x1,s1) = scale_number_density * numdensity(x1,s1)

    a12 = 2.0 * s1.exclusion_distance * outer_radius(s1)
    R = outer_radius(material.shape)
    k = ω / source.medium.c

    function calculate_gls_fun(;
            pair_corr::Function = hole_correction_pair_correlation,
            pair_corr_distance::Function = z ->  pair_corr(
                [0.0,0.0,0.0],material.species[1], [0.0,0.0,z],material.species[1]
            ),
            sigma_approximation = false,
            gls_pair_radial::Function = gls_pair_radial_fun(
                pair_corr_distance, a12;
                polynomial_order = polynomial_order,
                sigma_approximation = sigma_approximation
            ),
            h12 = R
        )
        gls_pair_radial, h12
    end

    gls_function, h12 = calculate_gls_fun(; pair_kws...)

    ls = 0:basis_order
    lm_to_n = lm_to_spherical_harmonic_index
    ls_to_ns = lm_to_n.(ls,0)

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rs = LinRange(0, R - outer_radius(s1), mesh_points)
    σs = integration_scheme(rs; scheme = scheme)

    # incident wave coefficients
    b0 = regular_spherical_coefficients(source)(1,origin(material.shape),ω)[1];
    # regular_radial_basis(source.medium, ω, basis_order, r1)[ls_to_ns]

    Bs = (sqrt(4pi) * b0) .* [
        t_diags[1][lm_to_n(l,0)] * (-T(1))^l * sbesselj(l, k*r1)
    for l in ls, r1 in rs][:]

    chi(l5::Int,l6::Int,r1::T,r2::T) = T(-1)^l6 *
        if(r1 < r2)
            shankelh1(l6, k*r2) * sbesselj(l5, k*r1)
        else
            shankelh1(l5, k*r1) * sbesselj(l6, k*r2)
        end

    function C_kernal(l2::Int,j2::Int)

        term2 = σs[j2] * rs[j2]^2 * bar_numdensity([T(0),T(0),rs[j2]],s1)
        a12 = 2*outer_radius(s1)*s1.exclusion_distance

        data = term2 .* [
        begin
            t_diags[1][lm_to_n(l,0)] *
            if r1 + rs[j2] < a12
                zero(Complex{T})
            elseif abs(r1 - rs[j2]) > h12
                chi(l,l2,r1,rs[j2]) * 4pi * T(-1)^l * T(2l2 + 1)
            else
                gls = gls_function(r1,rs[j2])
                L1 = polynomial_order
                sum(
                    gls[l1+1] *
                    sum(
                        (chi(l5,l6,r1,rs[j2]) * Complex{T}(im)^(l-l5-l6-l2) /
                        (4pi)^2) *
                        sum(
                            T(-1)^m2 *
                            gaunt_coefficient(l,0,l1,m1,l5,-m1) * gaunt_coefficient(l1,m1,l2,m2,l6,m1-m2) *
                            gaunt_coefficient(l2,m2,l,0,l4,m2) * gaunt_coefficient(l4,m2,l5,m1,l6,m2-m1)
                        for m1 = max(-l1,-l5):min(l1,l5)
                        for m2 = max(-l2,-l4,m1-l6):min(l2,l4,m1+l6))
                    for l4 = abs(l-l2):(l+l2)
                    for l5 = abs(l-l1):(l+l1)
                    for l6 = max(abs(l4-l5),abs(l1-l2)):min(l4+l5,l1+l2))
                for l1 = 0:L1)
            end
        end
        for l in ls, r1 in rs];

        return data[:]
    end

    Cs = [C_kernal(l2,j2) for l2 in ls, j2 in eachindex(rs)][:];

    bigC = hcat(Cs...);
    #NOTE bigC[:,M2] == Cs[M2]

    A = I - bigC

    Fs = A \ Bs;

    # NOTE: G * Fs = 0

    # CHECK:
    # norm(Fs - (Bs + bigC * Fs)) / norm(Bs)
    # norm(Fs - (Bs + bigC * Fs))

    # Fc = reshape(Fc,(basis_order+1,length(rs)))
    Fs = reshape(Fs,(basis_order+1,length(rs)))
    # reshape([(l,r) for l in ls, r in rs][:],(basis_order+1,length(rs)))

    # Approximate with a Legendre series
    P = Legendre{T}()

    r1_max = maximum(rs);
    rbars = T(2.0) .* rs ./ r1_max .- T(1.0);
    Pmat = P[rbars, 1:(legendre_order + 1)];

    # Pmat * pls_arr ~ tranpose(Fs)
    pls_arr = transpose(Pmat \ transpose(Fs))

    FsP = pls_arr * transpose(Pmat)
    # Fs ~ pls_arr * transpose(Pmat)

    if  norm(A * FsP[:] - Bs)/norm(Bs) > 0.001
        @warn "Numerical solution has a relative residual error of $(norm((I - bigC) * FsP[:] - Bs)/norm(Bs)), which is larger than expected. This residual error can be decreased by increasing the basis_field_order (current value: $basis_field_order) for the field."
    end

    function scattered_field(xs::AbstractVector{T})
        r,θ,φ = cartesian_to_radial_coordinates(xs)

        Ps = P[2r / r1_max - T(1), 1:(legendre_order + 1)]
        Fls = pls_arr * Ps

        # Fls_arr = [
        #     pls_arr * P[2r / r1_max - T(1), 1:(legendre_order + 1)]
        # for r in rs]
        # norm(Fs - hcat(Fls_arr...)) / norm(Fs)

        conj_Ys = conj.(spherical_harmonics(basis_order, θ, φ));
        fns = [
            Fls[l+1] * conj_Ys[lm_to_n(l,m)]
        for l = 0:basis_order for m = -l:l]

        return fns
    end

    # NOTE: need to think about this type, it is exactly the same as the discrete method above
    return ScatteringCoefficientsField(ω, source.medium, material, scattered_field;
        symmetry = RadialSymmetry{3}(),
        basis_order = basis_order,
        basis_field_order = basis_field_order
    )
end

"""
    outgoing_translation_matrix(ω, ::Acoustic, material::Material{Dim,Sphere{T,Dim}};
        basis_order = 2, tol = 1e-3)

    return a function U where U(X) for X ∈ material.shape gives the the outgoing translation matrix
"""
function outgoing_translation_matrix(ω::T, medium::Acoustic{T,Dim}, material::Material{Dim,Sphere{T,Dim}};
        basis_order::Int = 2, tol::T = 1e-3) where {T,Dim}

    L = basisorder_to_basislength(typeof(medium), basis_order);
    R = outer_radius(material.shape) - minimum(outer_radius.(material.species))
    a12 = T(2) * minimum(outer_radius(s) * s.exclusion_distance for s in material.species)

    rθφ2xyz = radial_to_cartesian_coordinates

    # Just to get the type Tinter.
        xs = LinRange(a12,2R,4);
        data = [zeros(Complex{T},L,L) for x in xs];
        inter = interpolate(
            [data[j][1] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        Tinterx = typeof(scale(inter, xs))

        data = [ zeros(Complex{T},L,L) for x in xs, y in xs, z in xs];
        inter = interpolate(
            [data[j][1] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        Tinterxyz = typeof(scale(inter, xs, xs, xs))

    # Determine the number of mesh points
        Ns = 10:4:66;
        Nθs = 10:4:62;

        xs = LinRange(a12,2R,Ns[end]);
        θs = LinRange(zero(T),T(pi),Nθs[end]);

        datax = [
            outgoing_translation_matrix(medium, basis_order, ω, [x,zero(T),zero(T)])
        for x in xs];

        dataθ = [
            outgoing_translation_matrix(medium, basis_order, ω, rθφ2xyz(SVector(a12,θ,zero(T))))
        for θ in θs];

        CI1s = CartesianIndices(datax[1]);
        inter_xs = Array{Tinterx}(undef, size(CI1s));
        inter_θs = Array{Tinterx}(undef, size(CI1s));

        for i in CI1s
            inter = interpolate(
                [datax[j][i] for j in CartesianIndices(datax)],
                BSpline(Cubic(Line(OnGrid())))
            );
            inter_xs[i] = scale(inter, xs)

            inter = interpolate(
                [dataθ[j][i] for j in CartesianIndices(dataθ)],
                BSpline(Cubic(Line(OnGrid())))
            );
            inter_θs[i] = scale(inter, θs)
        end

        errorsx = Array{T}(undef, length(Ns))
        for n in eachindex(Ns)
            xNs = LinRange(a12,2R,Ns[n]);
            inter_xNs = Array{Tinterx}(undef, size(CI1s));

            for i in CI1s
                inter = interpolate(inter_xs[i].(xNs),BSpline(Cubic(Line(OnGrid()))));
                inter_xNs[i] = scale(inter, xNs)
            end

            errorsx[n] = norm(norm(inter_xNs[i].(xs) - inter_xs[i].(xs)) for i in CI1s) / norm(norm(inter_xs[i].(xs)) for i in CI1s)
        end
        n = findfirst(errorsx .< tol)
        Nr = Ns[n]

        errorsθ = Array{T}(undef, length(Nθs))
        for n in eachindex(Nθs)
            θNs = LinRange(zero(T),T(pi),Nθs[n]);
            inter_θNs = Array{Tinterx}(undef, size(CI1s));

            for i in CI1s
                inter = interpolate(inter_θs[i].(θNs),BSpline(Cubic(Line(OnGrid()))));
                inter_θNs[i] = scale(inter, θNs)
            end

            errorsθ[n] = norm(norm(inter_θNs[i].(θs) - inter_θs[i].(θs)) for i in CI1s) / norm(norm(inter_θs[i].(θs)) for i in CI1s)
        end
        n = findfirst(errorsθ .< tol)
        Nθ = Nθs[n]

    # interpolate the whole spherical region based on the xs mesh size
    rs = LinRange(a12,2R,Nr);
    φs = LinRange(-T(pi),T(pi),Nθ);
    θs = LinRange(zero(T),T(pi),Nθ);

    data = [
            outgoing_translation_matrix(medium, basis_order, ω, rθφ2xyz(SVector(r,θ,φ)))
    for r in rs, θ in θs, φ in φs];

    # reorganise the data to interpolate each element of outgoing_translation_matrix in terms of x, y, z.
    CI1s = CartesianIndices(data[1]);
    inters = Array{Tinterxyz}(undef, size(CI1s));

    for i in CI1s
        inter = interpolate(
            [data[j][i] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        inters[i] = scale(inter, rs, θs, φs)
    end

    function Uout(X::AbstractVector{T})
        if norm(X) < a12
            zeros(Complex{T},L,L)
        else
            rθφ = cartesian_to_radial_coordinates(X)
            map(CI1s) do i
                inters[i](rθφ[1], rθφ[2], rθφ[3])
            end
        end
    end

    return Uout
end

max_stars_repo_name: arturgower/EffectiveWaves

max_issues_repo_name: arturgower/EffectiveWaves

max_forks_repo_name: arturgower/EffectiveWaves


text:
using InvertibleNetworks, NNlib, LinearAlgebra, Test

# Data
nx = 32
ny = 32
nz = 32
n_in = 3
batchsize = 4
k = 3   # kernel size
s = 1   # stride
p = 1   # padding

###################################################################################################
# Hyperbolic layer

# Data
X_prev = randn(Float32, nx, ny, nz, n_in, batchsize)
X_curr = randn(Float32, nx, ny, nz, n_in, batchsize)

# Network architecture ((decrease channel no?, no. of hidden units), ...)
architecture = ((0, 3), (-1, 3,), (0, 3), (1, 4), (0, 8))

# Hyperbolic network
HN = H = NetworkHyperbolic3D(n_in, architecture; k=3, s=1, p=1, logdet=true, α=1f0)

# Forward/inverse pass
Y_curr_, Y_new_, lgdet = HN.forward(X_prev, X_curr)
X_prev_, X_curr_ = HN.inverse(Y_curr_, Y_new_)

# Test invertibility
@test isapprox(norm(X_prev - X_prev_)/norm(X_prev_), 0f0; atol=1f-3)
@test isapprox(norm(X_curr - X_curr_)/norm(X_prev_), 0f0; atol=1f-3)


text:
struct ConsistentSampling <: BootstrapSampling
    nrun::Int
end

"""
    bootstrap(statistic, data, sampling::ConsistentSampling[;
              rng::AbstractRNG = Random.GLOBAL_RNG])

Obtain non-parametric bootstrap samples of `statistic` from the `data` set of predictions
and corresponding labels using consistency resampling, i.e., resampling of the labels.

Optionally you can specify a random number generator `rng` that is used for resampling
the data.
"""
function Bootstrap.bootstrap(statistic,
                             data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}},
                             sampling::ConsistentSampling;
                             rng::AbstractRNG = Random.GLOBAL_RNG)
    # check arguments
    predictions, labels = data
    nclasses, nsamples = size(predictions)
    nsamples == length(labels) ||
        throw(DimensionMismatch("number of predictions and labels must be equal"))

    # use same heuristic as StatsBase to decide whether to sample predictions
    # directly or to build an alias table
    # TODO: needs proper benchmarking
    if nsamples < 40
        bootstrap_direct(statistic, data, sampling; rng = rng)
    else
        t = nsamples < 500 ? 64 : 32
        if nclasses < t
            bootstrap_direct(statistic, data, sampling; rng = rng)
        else
            bootstrap_alias(statistic, data, sampling; rng = rng)
        end
    end
end

# sample labels directly (without alias table)
function bootstrap_direct(statistic,
                          data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}},
                          sampling::ConsistentSampling;
                          rng::AbstractRNG = Random.GLOBAL_RNG)
    # check arguments
    predictions, labels = data
    nclasses, nsamples = size(predictions)
    nsamples == length(labels) ||
        throw(DimensionMismatch("number of predictions and labels must be equal"))

    # evaluate statistic
    t0 = tx(statistic(data))

    # create caches
    resampled_predictions = similar(predictions)
    resampled_labels = similar(labels)
    resampled_data = (resampled_predictions, resampled_labels)

    # create sampler
    sp = Random.RangeGenerator(Base.OneTo(nsamples))

    # create output
    m = nrun(sampling)
    t1 = zeros_tuple(typeof(t0), m)

    # for each resampling step
    @inbounds for i in 1:m
        # resample data
        for j in 1:nsamples
            # resample predictions
            idx = rand(rng, sp)
            for k in axes(predictions, 1)
                resampled_predictions[k, j] = predictions[k, idx]
            end

            # resample labels
            p = rand(rng)
            cw = resampled_predictions[1, j]
            label = 1
            while cw < p && label < nclasses
                label += 1
                cw += resampled_predictions[label, j]
            end
            resampled_labels[j] = label
        end

        # evaluate statistic
        for (j, t) in enumerate(tx(statistic(resampled_data)))
            t1[j][i] = t
        end
    end

    NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

# sample labels with alias table
function bootstrap_alias(statistic,
                         data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}},
                         sampling::ConsistentSampling;
                         rng::AbstractRNG = Random.GLOBAL_RNG)
    # check arguments
    predictions, labels = data
    nclasses, nsamples = size(predictions)
    nsamples == length(labels) ||
        throw(DimensionMismatch("number of predictions and labels must be equal"))

    # evaluate statistic
    t0 = tx(statistic(data))

    # create alias table
    accept = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
    alias = [Vector{Int}(undef, nclasses) for _ in 1:nsamples]
    @inbounds for i in axes(predictions, 2)
        StatsBase.make_alias_table!(view(predictions, :, i), 1.0, accept[i], alias[i])
    end

    # create sampler of labels
    splabels = Random.RangeGenerator(Base.OneTo(nclasses))

    # create caches
    resampled_predictions = similar(predictions)
    resampled_labels = similar(labels)
    resampled_data = (resampled_predictions, resampled_labels)

    # create sampler
    sp = Random.RangeGenerator(Base.OneTo(nsamples))

    # create output
    m = nrun(sampling)
    t1 = zeros_tuple(typeof(t0), m)

    # for each resampling step
    @inbounds for i in 1:m
        # resample data
        for j in 1:nsamples
            # resample predictions
            idx = rand(rng, sp)
            for k in axes(predictions, 1)
                resampled_predictions[k, j] = predictions[k, idx]
            end

            # resample labels
            l = rand(rng, splabels)
            resampled_labels[j] = rand(rng) < accept[j][l] ? l : alias[j][l]
        end

        # evaluate statistic
        for (j, t) in enumerate(tx(statistic(resampled_data)))
            t1[j][i] = t
        end
    end

    NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

@deprecate bootstrap(statistic, predictions::AbstractMatrix{<:Real}, labels::AbstractVector{<:Integer}, sampling::ConsistentSampling) bootstrap(statistic, (predictions, labels), sampling)
@deprecate bootstrap(statistic, rng::AbstractRNG, predictions::AbstractMatrix{<:Real}, labels::AbstractVector{<:Integer}, sampling::ConsistentSampling) bootstrap(statistic, (predictions, labels), sampling; rng = rng)


text:
export get_config

function get_config(filename::String)
    cfg = YAML.load_file(filename)
    cfg["activation_functions"]=[]

    for k in keys(activ_dict)
        if cfg[k]
            push!(cfg["activation_functions"], activ_dict[k])
        end
        if cfg["hyperNEAT"] && cfg["hn_activation"]==k
            cfg["hn_activ_func"] = activ_dict[k]
        end
    end

    cfg["Species"]=Dict()

    cfg
end


text:
# a intercept-free ridge regressor for testing meta-alogorithms

export FooBarRegressor

import MLJBase
using LinearAlgebra
using MLJScientificTypes

mutable struct FooBarRegressor <: MLJBase.Deterministic
    lambda::Float64
end

function FooBarRegressor(; lambda=0.0)
    simpleridgemodel = FooBarRegressor(lambda)
    message = MLJBase.clean!(simpleridgemodel)
    isempty(message) || @warn message
    return simpleridgemodel
end

function MLJBase.clean!(model::FooBarRegressor)
    warning = ""
    if model.lambda < 0
        warning *= "Need lambda ≥ 0. Resetting lambda=0. "
        model.lambda = 0
    end
    return warning
end

function MLJBase.fitted_params(::FooBarRegressor, fitresult)
    return (coefficients=fitresult)
end

function MLJBase.fit(model::FooBarRegressor, verbosity::Int, X, y)
    x = MLJBase.matrix(X)
    fitresult = (x'x - model.lambda*I)\(x'y)
    cache = nothing
    report = NamedTuple()
    return fitresult, cache, report
end

function MLJBase.predict(model::FooBarRegressor, fitresult, Xnew)
    x = MLJBase.matrix(Xnew)
    return x*fitresult
end

# metadata:
MLJBase.load_path(::Type{<:FooBarRegressor}) = "MLJBase.FooBarRegressor"
MLJBase.package_name(::Type{<:FooBarRegressor}) = "MLJBase"
MLJBase.package_uuid(::Type{<:FooBarRegressor}) = ""
MLJBase.is_pure_julia(::Type{<:FooBarRegressor}) = true
MLJBase.input_scitype(::Type{<:FooBarRegressor}) = Table(Continuous)
MLJBase.target_scitype(::Type{<:FooBarRegressor}) = AbstractVector{Continuous}


text:
module SLATEC

export DRC, DRD, DRF, DRJ

const D1MACH1 = floatmin(Float64)
const D1MACH2 = floatmax(Float64)
const D1MACH3 = eps(Float64)/2
const D1MACH4 = eps(Float64)
const D1MACH5 = log10(2.)

#***BEGIN PROLOGUE  DRF
#***PURPOSE  Compute the incomplete or complete elliptic integral of the
#            1st kind.  For X, Y, and Z non-negative and at most one of
#            them zero, RF(X,Y,Z) = Integral from zero to infinity of
#                                -1/2     -1/2     -1/2
#                      (1/2)(t+X)    (t+Y)    (t+Z)    dt.
#            If X, Y or Z is zero, the integral is complete.
#***LIBRARY   SLATEC
#***CATEGORY  C14
#***TYPE      DOUBLE PRECISION (RF-S, DRF-D)
#***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
#             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE FIRST KIND,
#             TAYLOR SERIES
#***AUTHOR  Carlson, B. C.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Notis, E. M.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Pexton, R. L.
#             Lawrence Livermore National Laboratory
#             Livermore, CA  94550

function DRF(X::Float64, Y::Float64, Z::Float64)

    ERRTOL = (4.0*D1MACH3)^(1.0/6.0)
    LOLIM  = 5.0 * D1MACH1
    UPLIM  = D1MACH2/5.0
    C1 = 1.0/24.0
    C2 = 3.0/44.0
    C3 = 1.0/14.0

    ans = 0.0
    if min(X,Y,Z) < 0.0
        return ans, 1
    end

    if max(X,Y,Z) > UPLIM
        return ans, 3
    end

    if min(X+Y,X+Z,Y+Z) < LOLIM
        return ans, 2
    end

    XN = X
    YN = Y
    ZN = Z
    MU = 0.
    XNDEV = 0.
    YNDEV = 0.
    ZNDEV = 0.

    while true
        MU = (XN+YN+ZN)/3.0
        XNDEV = 2.0 - (MU+XN)/MU
        YNDEV = 2.0 - (MU+YN)/MU
        ZNDEV = 2.0 - (MU+ZN)/MU
        EPSLON = max(abs(XNDEV),abs(YNDEV),abs(ZNDEV))
        if (EPSLON < ERRTOL) break end
        XNROOT = sqrt(XN)
        YNROOT = sqrt(YN)
        ZNROOT = sqrt(ZN)
        LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
        XN = (XN+LAMDA)*0.250
        YN = (YN+LAMDA)*0.250
        ZN = (ZN+LAMDA)*0.250
    end

    E2 = XNDEV*YNDEV - ZNDEV*ZNDEV
    E3 = XNDEV*YNDEV*ZNDEV
    S  = 1.0 + (C1*E2-0.10-C2*E3)*E2 + C3*E3
    ans = S/sqrt(MU)

    return ans, 0
end

#***BEGIN PROLOGUE  DRD
#***PURPOSE  Compute the incomplete or complete elliptic integral of
#            the 2nd kind. For X and Y nonnegative, X+Y and Z positive,
#            DRD(X,Y,Z) = Integral from zero to infinity of
#                                -1/2     -1/2     -3/2
#                      (3/2)(t+X)    (t+Y)    (t+Z)    dt.
#            If X or Y is zero, the integral is complete.
#***LIBRARY   SLATEC
#***CATEGORY  C14
#***TYPE      DOUBLE PRECISION (RD-S, DRD-D)
#***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
#             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE SECOND KIND,
#             TAYLOR SERIES
#***AUTHOR  Carlson, B. C.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Notis, E. M.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Pexton, R. L.
#             Lawrence Livermore National Laboratory
#             Livermore, CA  94550

function DRD(X::Float64, Y::Float64, Z::Float64)

    ERRTOL = (D1MACH3/3.0)^(1.0/6.0)
    LOLIM  = 2.0/(D1MACH2)^(2.0/3.0)
    TUPLIM = D1MACH1^(1.0E0/3.0E0)
    TUPLIM = (0.10*ERRTOL)^(1.0E0/3.0E0)/TUPLIM
    UPLIM  = TUPLIM^2.0
    C1 = 3.0/14.0
    C2 = 1.0/6.0
    C3 = 9.0/22.0
    C4 = 3.0/26.0

    ans = 0.0
    if min(X,Y) < 0.0
        return ans, 1
    end

    if max(X,Y,Z) > UPLIM
        return ans, 3
    end

    if min(X+Y,Z) < LOLIM
        return ans, 2
    end

    XN = X
    YN = Y
    ZN = Z
    SIGMA = 0.0
    POWER4 = 1.0
    MU = 0.
    XNDEV = 0.
    YNDEV = 0.
    ZNDEV = 0.

    while true
        MU = (XN+YN+3.0*ZN)*0.20
        XNDEV = (MU-XN)/MU
        YNDEV = (MU-YN)/MU
        ZNDEV = (MU-ZN)/MU
        EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV))
        if (EPSLON < ERRTOL) break end
        XNROOT = sqrt(XN)
        YNROOT = sqrt(YN)
        ZNROOT = sqrt(ZN)
        LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
        SIGMA = SIGMA + POWER4/(ZNROOT*(ZN+LAMDA))
        POWER4 = POWER4*0.250
        XN = (XN+LAMDA)*0.250
        YN = (YN+LAMDA)*0.250
        ZN = (ZN+LAMDA)*0.250
    end

    EA = XNDEV*YNDEV
    EB = ZNDEV*ZNDEV
    EC = EA - EB
    ED = EA - 6.0*EB
    EF = ED + EC + EC
    S1 = ED*(-C1+0.250*C3*ED-1.50*C4*ZNDEV*EF)
    S2 = ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA))
    ans = 3.0*SIGMA + POWER4*(1.0+S1+S2)/(MU*sqrt(MU))

    return ans, 0
end

#***BEGIN PROLOGUE  DRC
#***PURPOSE  Calculate a double precision approximation to
#             DRC(X,Y) = Integral from zero to infinity of
#                              -1/2     -1
#                    (1/2)(t+X)    (t+Y)  dt,
#            where X is nonnegative and Y is positive.
#***LIBRARY   SLATEC
#***CATEGORY  C14
#***TYPE      DOUBLE PRECISION (RC-S, DRC-D)
#***KEYWORDS  DUPLICATION THEOREM, ELEMENTARY FUNCTIONS,
#             ELLIPTIC INTEGRAL, TAYLOR SERIES
#***AUTHOR  Carlson, B. C.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Notis, E. M.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Pexton, R. L.
#             Lawrence Livermore National Laboratory
#             Livermore, CA  94550

function DRC(X::Float64, Y::Float64)

    ERRTOL = (D1MACH3/16.0)^(1.0/6.0)
    LOLIM  = 5.0 * D1MACH1
    UPLIM  = D1MACH2 / 5.0
    C1 = 1.0/7.0
    C2 = 9.0/22.0

    ans = 0.0
    if X < 0.0 || Y <= 0.0
        return ans, 1
    end

    if max(X,Y) > UPLIM
        return ans, 3
    end

    if X+Y < LOLIM
        return ans, 2
    end

    XN = X
    YN = Y
    MU = 0.
    SN = 0.

    while true
        MU = (XN+YN+YN)/3.0
        SN = (YN+MU)/MU - 2.0
        if abs(SN) < ERRTOL break end
        LAMDA = 2.0*sqrt(XN)*sqrt(YN) + YN
        XN = (XN+LAMDA)*0.250
        YN = (YN+LAMDA)*0.250
    end

    S = SN*SN*(0.30+SN*(C1+SN*(0.3750+SN*C2)))
    ans = (1.0+S)/sqrt(MU)

    return ans,0
end

#***BEGIN PROLOGUE  DRJ
#***PURPOSE  Compute the incomplete or complete (X or Y or Z is zero)
#            elliptic integral of the 3rd kind.  For X, Y, and Z non-
#            negative, at most one of them zero, and P positive,
#             RJ(X,Y,Z,P) = Integral from zero to infinity of
#                              -1/2     -1/2     -1/2     -1
#                    (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt.
#***LIBRARY   SLATEC
#***CATEGORY  C14
#***TYPE      DOUBLE PRECISION (RJ-S, DRJ-D)
#***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
#             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE THIRD KIND,
#             TAYLOR SERIES
#***AUTHOR  Carlson, B. C.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Notis, E. M.
#             Ames Laboratory-DOE
#             Iowa State University
#             Ames, IA  50011
#           Pexton, R. L.
#             Lawrence Livermore National Laboratory
#             Livermore, CA  94550

function DRJ(X::Float64, Y::Float64, Z::Float64, P::Float64)

    ERRTOL = (D1MACH3/3.0)^(1.0/6.0)
    LOLIM  = (5.0 * D1MACH1)^(1.0/3.0)
    UPLIM  = 0.30*( D1MACH2 / 5.0)^(1.0/3.0)

    C1 = 3.0/14.0
    C2 = 1.0/3.0
    C3 = 3.0/22.0
    C4 = 3.0/26.0

    ans = 0.0
    if min(X,Y,Z) < 0.0
        return ans, 1
    end

    if max(X,Y,Z,P) > UPLIM
        return ans, 3
    end

    if min(X+Y,X+Z,Y+Z,P) < LOLIM
        return ans, 2
    end

    IER = 0
    XN = X
    YN = Y
    ZN = Z
    PN = P
    SIGMA = 0.0
    POWER4 = 1.0
    MU = 0.
    XNDEV = 0.
    YNDEV = 0.
    ZNDEV = 0.
    PNDEV = 0.

    while true
        MU = (XN+YN+ZN+PN+PN)*0.20
        XNDEV = (MU-XN)/MU
        YNDEV = (MU-YN)/MU
        ZNDEV = (MU-ZN)/MU
        PNDEV = (MU-PN)/MU
        EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV), abs(PNDEV))
        if EPSLON < ERRTOL break end
        XNROOT =  sqrt(XN)
        YNROOT =  sqrt(YN)
        ZNROOT =  sqrt(ZN)
        LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
        ALFA = PN*(XNROOT+YNROOT+ZNROOT) + XNROOT*YNROOT*ZNROOT
        ALFA = ALFA*ALFA
        BETA = PN*(PN+LAMDA)*(PN+LAMDA)
        drc,IER = DRC(ALFA,BETA)
        SIGMA = SIGMA + POWER4*drc
        POWER4 = POWER4*0.250
        XN = (XN+LAMDA)*0.250
        YN = (YN+LAMDA)*0.250
        ZN = (ZN+LAMDA)*0.250
        PN = (PN+LAMDA)*0.250
    end

    EA = XNDEV*(YNDEV+ZNDEV) + YNDEV*ZNDEV
    EB = XNDEV*YNDEV*ZNDEV
    EC = PNDEV*PNDEV
    E2 = EA - 3.0*EC
    E3 = EB + 2.0*PNDEV*(EA-EC)
    S1 = 1.0 + E2*(-C1+0.750*C3*E2-1.50*C4*E3)
    S2 = EB*(0.50*C2+PNDEV*(-C3-C3+PNDEV*C4))
    S3 = PNDEV*EA*(C2-PNDEV*C3) - C2*PNDEV*EC
    ans = 3.0*SIGMA + POWER4*(S1+S2+S3)/(MU* sqrt(MU))

    return ans, IER
end

end # module


text:
using Documenter
using TrajectoryOptimization

makedocs(
    sitename = "TrajectoryOptimization",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Introduction" => "index.md",
        "Getting Started" => [
            "models.md",
            "costfunctions.md",
            "constraints.md",
            "creating_problems.md",
            "solving.md"
        ],
        "Interfaces" => [
            "costfunction_interface.md",
            "constraint_interface.md",
            "solver_interface.md"
        ],
        "Documentation" => [
            "model_types.md",
            "discretization.md",
            "cost_api.md",
            "constraint_api.md",
            "problem.md",
            "solvers.md",
            "rotations.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/RoboticExplorationLab/TrajectoryOptimization.jl.git",
    deploy_config=Documenter.Travis(),
    devbranch = "v1.3",
)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

# Functions that can be called in getResidues! function

"""
    ModiaMath.getTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return actual **simulation time**.
"""
getTime(sim::SimulationState)                 = sim.time
getTime(m::ModiaMath.AbstractSimulationModel) = m.simulationState.time

"""
    ModiaMath.getStartTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return **startTime** of the actual simulation run.
"""
getStartTime(sim::SimulationState)                 = sim.startTime
getStartTime(m::ModiaMath.AbstractSimulationModel) = m.simulationState.startTime

"""
    ModiaMath.getStopTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return **stopTime** of the actual simulation run (when the simulation will be terminated).
"""
getStopTime(sim::SimulationState)                 = sim.stopTime
getStopTime(m::ModiaMath.AbstractSimulationModel) = m.simulationState.stopTime

"""
    ModiaMath.getTolerance(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return (relative) **tolerance** of the actual simulation run.
"""
getTolerance(sim::SimulationState)                 = sim.tolerance
getTolerance(m::ModiaMath.AbstractSimulationModel) = m.simulationState.tolerance

"""
    ModiaMath.isInitial(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if **initialization phase** of simulation.
"""
isInitial(sim::SimulationState)                 = sim.eventHandler.initial
isInitial(m::ModiaMath.AbstractSimulationModel) = m.simulationState.eventHandler.initial

"""
    ModiaMath.isTerminal(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if **terminal phase** of simulation.
"""
isTerminal(sim::SimulationState)                 = sim.eventHandler.terminal
isTerminal(m::ModiaMath.AbstractSimulationModel) = m.simulationState.eventHandler.terminal

"""
    ModiaMath.isEvent(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if **event phase** of simulation (including initialization)
"""
isEvent(sim::SimulationState)                 = sim.eventHandler.event
isEvent(m::ModiaMath.AbstractSimulationModel) = m.simulationState.eventHandler.event

"""
    ModiaMath.isZeroCrossing(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if simulation model shall compute zero-crossing functions (as required by the integrator).
"""
isZeroCrossing(sim::SimulationState)                 = sim.eventHandler.crossing
isZeroCrossing(m::ModiaMath.AbstractSimulationModel) = m.simulationState.eventHandler.crossing

"""
    ModiaMath.isAfterSimulationStart(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if after start of simulation and false if during initialization.
"""
isAfterSimulationStart(sim::SimulationState)                 = sim.eventHandler.afterSimulationStart
isAfterSimulationStart(m::ModiaMath.AbstractSimulationModel) = m.simulationState.eventHandler.afterSimulationStart

"""
    ModiaMath.isStoreResult(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if **communication point** and variables shall be stored in the result data structure.
"""
isStoreResult(sim::SimulationState)                 = sim.storeResult
isStoreResult(m::ModiaMath.AbstractSimulationModel) = m.simulationState.storeResult

ModiaMath.Logging.isLogInfos(sim::SimulationState) = ModiaMath.Logging.isLogInfos(sim.logger)
ModiaMath.Logging.isLogWarnings(sim::SimulationState) = ModiaMath.Logging.isLogWarnings(sim.logger)
ModiaMath.Logging.isLogEvents(sim::SimulationState) = ModiaMath.Logging.isLogEvents(sim.logger)

"""
    ModiaMath.setNominal!(sim::ModiaMath.SimulationState, x_nominal::Vector{Float64})

At initialization provide `x_nominal, the nominal values of the x-vector, and
store them in `sim`, the simulation state. These values are used to compute the
absolute tolerances of the x-vector for the integrator.

# Example

```julia
function getModelResidues!(m::Model, t, x, derx, r, w)
   sim = m.simulationState
   if ModiaMath.isInitial(sim)
      ModiaMath.setNominal!(sim, [1.0, 1e-5, 1e8])   # if length(x)=3
   end
   ...
end
```
"""
function setNominal!(sim::SimulationState, x_nominal::Vector{Float64})
    @assert(length(x_nominal) == sim.nx)
    @assert(x_nominal[ indmin(x_nominal) ] > 0.0)
    for i in eachindex(x_nominal)
        sim.x_nominal[i] = x_nominal[i]
    end
end

"""
    ModiaMath.setNextEvent!(sim, nextEventTime; integratoToEvent=true, restart = ModiaMath.Restart)

At an event instant (isEvent(sim) == true) trigger the next time event.

# Arguments

- `sim::ModiaMath.SimulationState`: (Internal) simulation state provided by simulation engine.
- `nextEventTime::Float64`: Time instant of the next time event.
- `integrateToEvent::Bool`: If true, the integrator integrates exactly to this event.
  If false, the integrator might integrate beyond `nextEventTime` (so the step size of the
  integrator is not influenced by this event). This option is only useful, if information is inquired
  at the event and restart=ModiaMath.NoRestart is used. The benefit is that the integrator is
  practically not influenced by this event.
- `restart::ModiaMath.`[`EventRestart`](@ref): Restart action after the event at which
  `setNextEvent!` was called.
"""
setNextEvent!(sim::SimulationState, nextEventTime::Float64;
              integrateToEvent::Bool=true, restart::EventRestart=Restart) =
              setNextEvent!(sim.eventHandler, nextEventTime; integrateToEvent=integrateToEvent, restart=restart)

"""
    ModiaMath.positive!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Return `crossing > 0` such that a state event is triggered whenever `crossing > 0` changes its value.

Note, a small hysteresis is added to the crossing function during continuous integration,
in order to reduce issues during event iterations due to small numerical errors.
However, at an event instant, `crossing > 0` is returned without any hysteresis.

# Arguments

- `sim::ModiaMath.SimulationState`: (Internal) simulation state provided by simulation engine.
- `nr::Int`: (> 0 required) Every call of `positive!(..)` must be identified by a unique value of `nr`
        This value is used as index in a vector that holds the internal memory for `crossing > 0`.
- `crossing::Float64`: Zero crossing function.
- `crossingAsString::String`: `crossing` as string representation. This string is used for log messages.
- `restart::ModiaMath.`[`EventRestart`](@ref): Restart action after the `crossing > 0` event occurred.

# Example

```julia
function computeVariables!(m::Model, sim::SimulationState)
   ...
   # f = s > 0 ? fMax : 0.0
   m.sPos.value = ModiaMath.positive!(sim, 1, m.s.value, "s")
   m.f.value    = m.sPos.value ? m.fMax.value : 0.0
   ...
end
```
"""
positive!(sim::SimulationState, nr::Int, crossing::Float64, crossingAsString::String; restart::EventRestart=Restart) =
              positive!(sim.eventHandler, nr, crossing, crossingAsString; restart=restart)

"""
    ModiaMath.negative!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Return `crossing < 0` such that a state event is triggered whenever `crossing < 0` changes its value.
"""
negative!(sim::SimulationState, nr::Int, crossing::Float64, crossingAsString::String; restart::EventRestart=Restart) =
              negative!(sim.eventHandler, nr, crossing, crossingAsString; restart=restart)

"""
    ModiaMath.change!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Trigger an event, whenever `crossing > 0` changes from `false` to `true` or from `true` to `false`.
The function returns always `false`.
"""
change!(sim::SimulationState, nr::Int, crossing::Float64, crossingAsString::String; restart::EventRestart=Restart) =
              change!(getDaeInfo(model).eventHandler, nr, crossing, crossingAsString; restart=restart)

"""
    ModiaMath.edge!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Trigger an event, whenever `crossing > 0` switches from `false` to `true`.
The function returns always `false`.
"""
edge!(sim::SimulationState, nr::Int, crossing::Float64, crossingAsString::String; restart::EventRestart=Restart) =
              edge!(sim.eventHandler, nr, crossing, crossingAsString; restart=restart)

ModiaMath.Logging.logOn!(sim::SimulationState)                     = ModiaMath.Logging.logOn!(sim.logger)
ModiaMath.Logging.logOn!(model::ModiaMath.AbstractSimulationModel) = ModiaMath.Logging.logOn!(model.simulationState.logger)

ModiaMath.Logging.logOff!(sim::SimulationState)                     = ModiaMath.Logging.logOff!(sim.logger)
ModiaMath.Logging.logOff!(model::ModiaMath.AbstractSimulationModel) = ModiaMath.Logging.logOff!(model.simulationState.logger)

"""
    ModiaMath.setLogCategories!(m::ModiaMath.[AbstractSimulationModel|SimulationState], categories; reinit=true)

Set log categories as vector of symbols, e.g. setLogCategories!(logger, [:LogProgess]).
Supported categories:
- `:LogStatistics`, print statistics information at end of simulation
- `:LogProgress`, print progress information during simulation
- `:LogInfos`, print information messages of the model
- `:LogWarnings`, print warning messages of the model
- `:LogEvents`, log events of the model

If option reinit=true, all previous set categories are reinitialized to be no longer present.
If reinit=false, previously set categories are not changed.
"""
ModiaMath.Logging.setLogCategories!(sim::SimulationState, categories::Vector{Symbol}; reinit::Bool=true) =
                  ModiaMath.Logging.setLogCategories!(sim.logger, categories; reinit=reinit)
ModiaMath.Logging.setLogCategories!(model::ModiaMath.AbstractSimulationModel, categories::Vector{Symbol}; reinit::Bool=true) =
                  ModiaMath.Logging.setLogCategories!(model.simulationState.logger, categories; reinit=reinit)

# Functions that can be called after a model is instantiated
getSimulationResult(model::ModiaMath.AbstractSimulationModel) = model.simulationState.result


text:
# Reconstruct a digit from an array
function nd(a)
  multiplicand,result = 1,0
  for i in reverse(1:length(a))
    result += multiplicand * a[i]
    multiplicand *= 10
  end
  result
end

# Check if one array is a permutation of another
function ispermutation(a, b)
  sort(a) == sort(b)
end

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler