text:
# Useful functions from Base

function Base.rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:AbstractQuaternion,S<:AbstractQuaternion}
    Base.rtoldefault(eltype(x), eltype(y), atol)
end
function Base.rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:AbstractQuaternion,S<:Number}
    Base.rtoldefault(eltype(x), y, atol)
end
function Base.rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:AbstractQuaternion}
    Base.rtoldefault(x, eltype(y), atol)
end

function Base.:(==)(q1::AbstractQuaternion{Symbolics.Num}, q2::AbstractQuaternion{Symbolics.Num})
    (
        iszero(Symbolics.simplify(q1.w-q2.w; expand=true)) &&
        iszero(Symbolics.simplify(q1.x-q2.x; expand=true)) &&
        iszero(Symbolics.simplify(q1.y-q2.y; expand=true)) &&
        iszero(Symbolics.simplify(q1.z-q2.z; expand=true))
    )
end
function Base.:(==)(q1::AbstractQuaternion{Symbolics.Num}, q2::Real)
    (
        iszero(Symbolics.simplify(q1.w-q2; expand=true)) &&
        iszero(Symbolics.simplify(q1.x; expand=true)) &&
        iszero(Symbolics.simplify(q1.y; expand=true)) &&
        iszero(Symbolics.simplify(q1.z; expand=true))
    )
end
function Base.:(==)(q1::Real, q2::AbstractQuaternion{Symbolics.Num})
    (
        iszero(Symbolics.simplify(q1-q2.w; expand=true)) &&
        iszero(Symbolics.simplify(q2.x; expand=true)) &&
        iszero(Symbolics.simplify(q2.y; expand=true)) &&
        iszero(Symbolics.simplify(q2.z; expand=true))
    )
end
function Base.:(==)(q1::AbstractQuaternion{Symbolics.Num}, q2::Symbolics.Num)
    (
        iszero(Symbolics.simplify(q1.w-q2; expand=true)) &&
        iszero(Symbolics.simplify(q1.x; expand=true)) &&
        iszero(Symbolics.simplify(q1.y; expand=true)) &&
        iszero(Symbolics.simplify(q1.z; expand=true))
    )
end
function Base.:(==)(q1::Symbolics.Num, q2::AbstractQuaternion{Symbolics.Num})
    (
        iszero(Symbolics.simplify(q1-q2.w; expand=true)) &&
        iszero(Symbolics.simplify(q2.x; expand=true)) &&
        iszero(Symbolics.simplify(q2.y; expand=true)) &&
        iszero(Symbolics.simplify(q2.z; expand=true))
    )
end
Base.:(==)(q1::AbstractQuaternion{<:Real}, q2::AbstractQuaternion{<:Real}) = (q1.w==q2.w) && (q1.x==q2.x) && (q1.y==q2.y) && (q1.z==q2.z)
Base.:(==)(q::AbstractQuaternion{<:Real}, x::Real) = isreal(q) && real(q) == x
Base.:(==)(x::Real, q::AbstractQuaternion) = isreal(q) && real(q) == x

Base.isequal(q1::AbstractQuaternion, q2::AbstractQuaternion) = isequal(q1.w,q2.w) && isequal(q1.x,q2.x) && isequal(q1.y,q2.y) && isequal(q1.z,q2.z)

Base.isreal(q::AbstractQuaternion{T}) where {T<:Real} = iszero(q.x) && iszero(q.y) && iszero(q.z)
Base.isinteger(q::AbstractQuaternion{T}) where {T<:Real} = isreal(q) && isinteger(real(q))
Base.isfinite(q::AbstractQuaternion{T}) where {T<:Real} = isfinite(q.w) && isfinite(q.x) && isfinite(q.y) && isfinite(q.z)
Base.isnan(q::AbstractQuaternion{T}) where {T<:Real} = isnan(q.w) || isnan(q.x) || isnan(q.y) || isnan(q.z)
Base.isinf(q::AbstractQuaternion{T}) where {T<:Real} = isinf(q.w) || isinf(q.x) || isinf(q.y) || isinf(q.z)
Base.iszero(q::AbstractQuaternion{T}) where {T<:Real} = iszero(q.w) && iszero(q.x) && iszero(q.y) && iszero(q.z)
Base.isone(q::AbstractQuaternion{T}) where {T<:Real} = isone(q.w) && iszero(q.x) && iszero(q.y) && iszero(q.z)

Base.in(q::AbstractQuaternion, r::AbstractRange{<:Real}) = isreal(q) && real(q) in r

Base.flipsign(q::AbstractQuaternion, x::Real) = ifelse(signbit(x), -q, q)

Base.bswap(q::Q) where {T<:Real, Q<:AbstractQuaternion{T}} = Q(bswap(q.w), bswap(q.x), bswap(q.y), bswap(q.z))

if UInt === UInt64
    const h_imagx = 0xdf13da9384000582
    const h_imagy = 0x437d0726f1028bcd
    const h_imagz = 0xcf13f7ab1f367e01
else
    const h_imagx = 0x27a4bf84
    const h_imagy = 0xccefdeeb
    const h_imagz = 0x1683854f
end
const hash_0_imagx = hash(0, h_imagx)
const hash_0_imagy = hash(0, h_imagy)
const hash_0_imagz = hash(0, h_imagz)

function Base.hash(q::AbstractQuaternion, h::UInt)
    # TODO: with default argument specialization, this would be better:
    # hash(q.w, h ⊻ hash(q.x, h ⊻ h_imagx) ⊻ hash(0, h ⊻ h_imagx) ⊻ hash(q.y, h ⊻ h_imagy) ⊻ hash(0, h ⊻ h_imagy) ⊻ hash(q.z, h ⊻ h_imagz) ⊻ hash(0, h ⊻ h_imagz))
    hash(q.w, h ⊻ hash(q.x, h_imagx) ⊻ hash_0_imagx ⊻ hash(q.y, h_imagy) ⊻ hash_0_imagy ⊻ hash(q.z, h_imagz) ⊻ hash_0_imagz)
end

function Base.show(io::IO, q::AbstractQuaternion)
    function pm(x)
        s = "$x"
        if s[1] ∉ "+-"
            s = "+" * s
        end
        if occursin(r"[+^/-]", s[2:end])
            s = " " * s[1] * " " * "(" * s[2:end] * ")"
        else
            s = " " * s[1] * " " * s[2:end]
        end
        s
    end
    print(
        io,
        q.w,
        pm(q.x), "𝐢",
        pm(q.y), "𝐣",
        pm(q.z), "𝐤"
    )
end

function Base.show(io::IO, ::MIME"text/latex", q::AbstractQuaternion)
    function pm(x)
        s = latexify(x, env=:raw, bracket=true)
        if s[1] ∉ "+-"
            s = "+" * s
        end
        if occursin(r"[+^/-]", s[2:end])
            s = " " * s[1] * " " * "\\left(" * s[2:end] * "\\right)"
        else
            s = " " * s[1] * " " * s[2:end]
        end
        s
    end
    s = latexstring(
        latexify(q.w, env=:raw, bracket=true),
        pm(q.x), "\\,\\mathbf{i}",
        pm(q.y), "\\,\\mathbf{j}",
        pm(q.z), "\\,\\mathbf{k}"
    )
    print(io, s)
end

function Base.read(s::IO, QT::Type{Q}) where {T<:Real, Q<:AbstractQuaternion{T}}
    w = read(s,T)
    x = read(s,T)
    y = read(s,T)
    z = read(s,T)
    QT(w,x,y,z)
end

function Base.write(s::IO, q::AbstractQuaternion)
    write(s,q.w,q.x,q.y,q.z)
end

#Broadcast.broadcasted(f, q::QT, args...) where {QT<:AbstractQuaternion{<:Real}} = wrapper(QT)(f.(q.components, args...))


text:
using FileIO

"""
    save_pcof(refFileName, pcof)

Save the parameter vector `pcof` on a JLD2 formatted file with handle `pcof`

# Arguments
- `refFileName`: String holding the name of the file.
- `pcof`: Vector of floats holding the parameters.
"""
function save_pcof(refFileName:: String, pcof:: Vector{Float64})
    save(refFileName, "pcof", pcof)
end

"""
    pcof = read_pcof(refFileName)

Read the parameter vector `pcof` from a JLD2 formatted file

# Arguments
- `refFileName`: String holding the name of the file.
"""
function read_pcof(refFileName:: String)
    dict = load(refFileName)
    pcof = dict["pcof"]
    return pcof
end

"""
	juq2qis(params, pcof, samplerate, q_ind, fileName="ctrl.dat", node_loc="c")

Evaluate control functions and export them into a format that is readable by Qiskit.

# Arguments
- `params:: objparams`: Struct with problem definition
- `pcof:: Array{Float64,1})`: Vector of parameter values
- `samplerate:: Float64`: Samplerate for quantum device (number of samples per ns for the IQ mixer)
- `q_ind:: Int64`: Index of the control function to output (`1 <= q_ind <= Nctrl*Nfreq`)
- `fileName:: String`: Name of output file containing controls to be handled by Qiskit
- `node_loc:: String`: Node location, "c" for cell centered, "n" for node-centered, default is cell-centered
"""
function juq2qis(params::objparams, pcof:: Array{Float64, 1}, samplerate:: Float64, q_ind:: Int64, fileName:: String="ctrl.dat", node_loc:: String="c")

    Nctrl = params.Ncoupled + params.Nunc
    Nfreq = params.Nfreq

    qiskit_dt = 1/samplerate
    @assert(q_ind >=1)
    @assert(q_ind <= Nctrl*Nfreq)

    if params.T < qiskit_dt
	throw(ArgumentError("Final simulation time shorter than qiskit_dt, please increase final simulation time.\n"))
    end

    @assert( params.Ncoupled == 0 || params.Nunc == 0) # can't have both

    # Check if we have an integer number of timesteps of size qiskit_dt in our signal
    # FG: This can be a problem if the width of a B-spline is small enough as
    # the qiskit_dt could ask for data outside of the control function definition
    # if we don't take value from left.
    nsteps = round(Int64,params.T/qiskit_dt)
    rmdr = params.T-qiskit_dt*nsteps
    println("Number of ctrl samples = ", nsteps)
    if rmdr > 1e-10
	throw(ArgumentError("Non-integral number of timesteps of size qiskit_dt, adjust final simulation time.\n"))
    end

    # Take midpoint value of signal in interval [(k-1)*qiskit_dt,k*qiskit_dt] ∀ k
    if(node_loc=="c")
        td = collect((0:nsteps-1).*qiskit_dt)
	td .+= 0.5*qiskit_dt  # Use the ctrl value at the center of each interval
    else
        td = collect((0:nsteps-1).*qiskit_dt)     # Use the ctrl value at the beginning of each interval
    end
    D  = zeros(length(td),2)

    # output the B-splines without carrier waves
    bc_flag = params.use_bcarrier
    params.use_bcarrier = false

    # Controls (in rad/ns)
    p,q = Juqbox.evalctrl(params, pcof, td, q_ind) # evalctrl uses a 1-based indexing of q_ind
    D[:,1] = copy(p)
    D[:,2] = copy(q)

    println("Making a plot of the ctrl functions")
    plc = plot(td,p,lab="p(t)")
    plot!(td,q,lab="q(t)")

    # reset the bcarrier flag
    params.use_bcarrier = bc_flag

    # Save signal to delimited file
    if length(fileName)>0
    	open(fileName, "w") do io
	    writedlm(io, D)
	end
        println("Saved control signal in Qiskit compatible format to delimited file '", fileName, "'");
    end

    return plc
end

ext: jl


text:
"""
`SimpleTools` is a collection of useful functions that don't have an
obvious home, but are useful for my `SimpleWorld` modules.
"""
module SimpleTools

include("dict_compose.jl")
include("make_continuous.jl")
include("flush_print.jl")
include("sage.jl")
include("dcat.jl")
include("my_eye.jl")
include("cofactor_det.jl")
include("char_poly.jl")
end # end of module


text:
module Integrals

using Reexport,  MonteCarloIntegration, QuadGK, HCubature
@reexport using SciMLBase
using Zygote, ReverseDiff, ForwardDiff , LinearAlgebra

import ChainRulesCore
import ChainRulesCore: NoTangent
import ZygoteRules

struct QuadGKJL <: SciMLBase.AbstractIntegralAlgorithm end
struct HCubatureJL <: SciMLBase.AbstractIntegralAlgorithm end
struct VEGAS <: SciMLBase.AbstractIntegralAlgorithm
    nbins::Int
    ncalls::Int
end

VEGAS(;nbins = 100,ncalls = 1000) = VEGAS(nbins,ncalls)

abstract type QuadSensitivityAlg end
struct ReCallVJP{V}
    vjp::V
end

abstract type IntegralVJP end
struct ZygoteVJP end
struct ReverseDiffVJP
    compile::Bool
end

function scale_x!(_x,ub,lb,x)
    _x .= (ub .- lb) .* x .+ lb
    _x
end

function scale_x(ub,lb,x)
    (ub .- lb) .* x .+ lb
end

function v_inf(t)
	return t ./ (1 .- t.^2)
end

function v_semiinf(t , a , upto_inf)
	if upto_inf == true
		return a .+ (t ./ (1 .- t))
	else
		return a .+ (t ./ (1 .+ t))
	end
end
function transfor_inf_number(t , p , f , lb , ub)
	if lb == -Inf && ub == Inf
		j = (1 .+ t.^2 )/(1 .- t.^2).^2
		return f(v_inf(t) , p)*(j)
	elseif lb != -Inf && ub == Inf
		a = lb
		j = 1 ./ ((1 .- t).^2)
		return f(v_semiinf(t , a , 1) , p)*(j)
	elseif lb == -Inf && ub != Inf
		a = ub
		j = 1 ./ ((1 .+ t ).^2)
		return f(v_semiinf(t , a , 0) , p)*(j)
	end
end

function transform_inf(t , p , f , lb , ub)
	if lb isa Number && ub isa Number
		return transfor_inf_number(t , p , f , lb , ub)
	end

	lbb = lb .== -Inf
	ubb = ub .== Inf
	_none = .!lbb .& .!ubb
	_inf = lbb .& ubb
	semiup = .!lbb .& ubb
	semilw = lbb  .& .!ubb

	function v(t)
		return t.*_none + v_inf(t).*_inf + v_semiinf(t , lb , 1).*semiup + v_semiinf(t , ub , 0).*semilw
	end
    jac = Zygote.@ignore ForwardDiff.jacobian(x -> v(x), t)
	j = det(jac)
	f(v(t) , p)*(j)
end

function transformation_if_inf(prob)
	if (prob.lb isa Number && prob.ub isa Number && (prob.ub == Inf || prob.lb == -Inf))
		if prob.lb == -Inf && prob.ub == Inf
			g = prob.f
			h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = -1.00
			ub = 1.00
	        _prob = IntegralProblem(h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		elseif prob.lb != -Inf && prob.ub == Inf
			g = prob.f
			h_(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = 0.00
			ub = 1.00
	        _prob = IntegralProblem(h_ , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		elseif prob.lb == -Inf && prob.ub != Inf
			g = prob.f
			_h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = -1.00
			ub = 0.00
	        _prob = IntegralProblem(_h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		end
	end
	if -Inf in prob.lb || Inf in prob.ub
		lbb = prob.lb .== -Inf
		ubb = prob.ub .== Inf
		_none = .!lbb .& .!ubb
		_inf = lbb .& ubb
		_semiup = .!lbb .& ubb
		_semilw = lbb  .& .!ubb

		_lb = 0.00.*_semiup + -1.00.*_inf + -1.00.*_semilw +  _none.*prob.lb
		_ub = 1.00.*_semiup + 1.00.*_inf  + 0.00.*_semilw  + _none.*prob.ub
		g = prob.f
		hs(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
		_prob = IntegralProblem(hs, _lb ,_ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
		return _prob
	end
	return prob
end
function SciMLBase.solve(prob::IntegralProblem,::Nothing,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8, kwargs...)
    if lb isa Number
        __solve(prob,QuadGKJL();reltol=reltol,abstol=abstol,kwargs...)
    elseif length(lb) > 8 && reltol < 1e-4 || abstol < 1e-4
        __solve(prob,VEGAS();reltol=reltol,abstol=abstol,kwargs...)
    else
        __solve(prob,HCubatureJL();reltol=reltol,abstol=abstol,kwargs...)
    end
end

function SciMLBase.solve(prob::IntegralProblem,
                            alg::SciMLBase.AbstractIntegralAlgorithm,
                            args...; sensealg = ReCallVJP(ZygoteVJP()), kwargs...)
   prob = transformation_if_inf(prob)
  __solvebp(prob,alg,sensealg,prob.lb,prob.ub,prob.p,args...;kwargs...)
end

# Give a layer to intercept with AD
__solvebp(args...;kwargs...) = __solvebp_call(args...;kwargs...)

function __solvebp_call(prob::IntegralProblem,::QuadGKJL,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    if isinplace(prob) || lb isa AbstractArray || ub isa AbstractArray
        error("QuadGKJL only accepts one-dimensional quadrature problems.")
    end
    @assert prob.batch == 0
    @assert prob.nout == 1
    p = p
    f = x -> prob.f(x,p)
    val,err = quadgk(f, lb, ub,
                     rtol=reltol, atol=abstol,
                     kwargs...)
    SciMLBase.build_solution(prob,QuadGKJL(),val,err,retcode = :Success)
end

function __solvebp_call(prob::IntegralProblem,::HCubatureJL,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    p = p

    if isinplace(prob)
        dx = zeros(prob.nout)
        f = (x) -> (prob.f(dx,x,p); dx)
    else
        f = (x) -> prob.f(x,p)
    end
    @assert prob.batch == 0

    if lb isa Number
        val,err = hquadrature(f, lb, ub;
                              rtol=reltol, atol=abstol,
                              maxevals=maxiters, kwargs...)
    else
        val,err = hcubature(f, lb, ub;
                            rtol=reltol, atol=abstol,
                            maxevals=maxiters, kwargs...)
    end
    SciMLBase.build_solution(prob,HCubatureJL(),val,err,retcode = :Success)
end

function __solvebp_call(prob::IntegralProblem,alg::VEGAS,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    p = p
    @assert prob.nout == 1
    if prob.batch == 0
        if isinplace(prob)
          dx = zeros(prob.nout)
          f = (x) -> (prob.f(dx,x,p); dx)
        else
          f = (x) -> prob.f(x,p)
        end
    else
        if isinplace(prob)
          dx = zeros(prob.batch)
          f = (x) -> (prob.f(dx,x',p); dx)
        else
          f = (x) -> prob.f(x',p)
        end
    end
    val,err,chi = vegas(f, lb, ub, rtol=reltol, atol=abstol,
                        maxiter = maxiters, nbins = alg.nbins,
                        ncalls = alg.ncalls, batch=prob.batch != 0, kwargs...)
    SciMLBase.build_solution(prob,alg,val,err,chi=chi,retcode = :Success)
end

function ChainRulesCore.rrule(::typeof(__solvebp),prob,alg,sensealg,lb,ub,p,args...;kwargs...)
    out = __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
    function quadrature_adjoint(Δ)
        y = typeof(Δ) <: Array{<:Number,0} ? Δ[1] : Δ
        if isinplace(prob)
            dx = zeros(prob.nout)
            _f = (x) -> prob.f(dx,x,p)
            if sensealg.vjp isa ZygoteVJP
                dfdp = function (dx,x,p)
                    _,back = Zygote.pullback(p) do p
                        _dx = Zygote.Buffer(x, prob.nout, size(x,2))
                        prob.f(_dx,x,p)
                        copy(_dx)
                    end

                    z = zeros(size(x,2))
                    for idx in 1:size(x,2)
                        z[1] = 1
                        dx[:,idx] = back(z)[1]
                        z[idx]=0
                    end
                end
            elseif sensealg.vjp isa ReverseDiffVJP
                error("TODO")
            end
        else
            _f = (x) -> prob.f(x,p)
            if sensealg.vjp isa ZygoteVJP
                if prob.batch > 0
                    dfdp = function (x,p)
                        _,back = Zygote.pullback(p->prob.f(x,p),p)

                        out = zeros(length(p),size(x,2))
                        z = zeros(size(x,2))
                        for idx in 1:size(x,2)
                            z[idx] = 1
                            out[:,idx] = back(z)[1]
                            z[idx]=0
                        end
                        out
                    end
                else
                    dfdp = function (x,p)
                        _,back = Zygote.pullback(p->prob.f(x,p),p)
                        back(y)[1]
                    end
                end

            elseif sensealg.vjp isa ReverseDiffVJP
                error("TODO")
            end
        end

        dp_prob = IntegralProblem(dfdp,lb,ub,p;nout=length(p),batch = prob.batch,kwargs...)

        if p isa Number
            dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...)[1]
        else
            dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...).u
        end

        if lb isa Number
            dlb = -_f(lb)
            dub = _f(ub)
            return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),dlb,dub,dp,ntuple(x->NoTangent(),length(args))...)
        else
            return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),dp,ntuple(x->NoTangent(),length(args))...)
        end
    end
    out,quadrature_adjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(
                        sol::SciMLBase.IntegralSolution,::Val{:u})
    sol.u,Δ->(SciMLBase.build_solution(sol.prob,sol.alg,Δ,sol.resid),)
end

### Forward-Mode AD Intercepts

# Direct AD on solvers with QuadGK and HCubature
function __solvebp(prob,alg::QuadGKJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
end

function __solvebp(prob,alg::HCubatureJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
end

# Manually split for the pushforward
function __solvebp(prob,alg,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    primal = __solvebp_call(prob,alg,sensealg,lb,ub,ForwardDiff.value.(p),args...;kwargs...)

    nout = prob.nout*P

    if isinplace(prob)
        dfdp = function (out,x,p)
            dualp = reinterpret(ForwardDiff.Dual{T,V,P}, p)
            if prob.batch > 0
                dx = similar(dualp, prob.nout, size(x,2))
            else
                dx = similar(dualp, prob.nout)
            end
            prob.f(dx,x,dualp)

            ys = reinterpret(ForwardDiff.Dual{T,V,P}, dx)
            idx=0
            for y in ys
                for p in ForwardDiff.partials(y)
                    out[idx+=1] = p
                end
            end
            return out
        end
    else
        dfdp = function (x,p)
            dualp = reinterpret(ForwardDiff.Dual{T,V,P}, p)
            ys = prob.f(x,dualp)
            if prob.batch > 0
                out = similar(p, V, nout, size(x,2))
            else
                out = similar(p, V, nout)
            end

            idx=0
            for y in ys
                for p in ForwardDiff.partials(y)
                    out[idx+=1] = p
                end
            end

            return out
        end
    end
    rawp = copy(reinterpret(V, p))

    dp_prob = IntegralProblem(dfdp,lb,ub,rawp;nout=nout,batch = prob.batch,kwargs...)
    dual = __solvebp_call(dp_prob,alg,sensealg,lb,ub,rawp,args...;kwargs...)
    res = similar(p, prob.nout)
    partials = reinterpret(typeof(first(res).partials), dual.u)
    for idx in eachindex(res)
        res[idx] = ForwardDiff.Dual{T,V,P}(primal.u[idx], partials[idx])
    end
    if primal.u isa Number
        res = first(res)
    end
    SciMLBase.build_solution(prob,alg,res,primal.resid)
end

export QuadGKJL, HCubatureJL, VEGAS
end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.263105417837706 <= q <= 2.4766737340949954)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.7635571507341847 + -0.35250618093605723*1/(1 + exp(-(0.5051094873509578 + -0.17934163301535788*1/(1 + exp(-(-0.331542222824047 + 0.7450661875335634*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + -0.10668129539960924*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))) + -0.9610191438924418*1/(1 + exp(-(0.6781688636512797 + -0.22580923517577522*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + 0.2965219706242044*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))))))))) + 0.1709251594821839*1/(1 + exp(-(0.058205314222445015 + 0.7954955842799771*1/(1 + exp(-(-0.331542222824047 + 0.7450661875335634*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + -0.10668129539960924*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))) + -0.9762293687340509*1/(1 + exp(-(0.6781688636512797 + -0.22580923517577522*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + 0.2965219706242044*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))))))))) + 1/(1 + exp(-(-0.29663194049169395 + 0.013502946151282558*1/(1 + exp(-(0.5051094873509578 + -0.17934163301535788*1/(1 + exp(-(-0.331542222824047 + 0.7450661875335634*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + -0.10668129539960924*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))) + -0.9610191438924418*1/(1 + exp(-(0.6781688636512797 + -0.22580923517577522*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + 0.2965219706242044*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))))))))) + -0.2459643040595232*1/(1 + exp(-(0.058205314222445015 + 0.7954955842799771*1/(1 + exp(-(-0.331542222824047 + 0.7450661875335634*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + -0.10668129539960924*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))) + -0.9762293687340509*1/(1 + exp(-(0.6781688636512797 + -0.22580923517577522*1/(1 + exp(-(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])))) + 0.2965219706242044*1/(1 + exp(-(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# types for dispatch
abstract type RewardCall end
abstract type jit <: RewardCall end
abstract type pre <: RewardCall end
abstract type pre_partial <: RewardCall end

struct SolverOptions{RC<:RewardCall,  ID<:IntDim, M<:Union{Bool,Vector{Bool}},
		T<:Union{Nothing, Array{Float64,2}}}
	disp::Bool
	disp_each_iter::Int
	max_iter::Int
	epsilon::Float64
	rewardcall::Type{RC}
	monotonicity::M
	concavity::M
	numquadnodes::Vector{Int}
	mTransition::T
	intdim::Type{ID}
end

function SolverOptions(p::DDP; disp::Bool = false,
	disp_each_iter::Int = 10,
	max_iter::Int = 500,
	epsilon::Float64 = 1e-3,
	rewardcall::Symbol = :jit,
	monotonicity::Union{Bool,Vector{Bool}} = false,
	concavity::Union{Bool,Vector{Bool}} = false,
	numquadnodes::Vector{Int} = 5*ones(Int64, length(p.shockdist.μ)),
	mTransition::Union{Nothing, Array{Float64,2}} = nothing,
	intdim::Type{T} = p.intdim) where # if want to override type in the model, mostly for testing)
		T<:IntDim

	SolverOptions(disp, disp_each_iter, max_iter, epsilon, eval(rewardcall),
		monotonicity, concavity, numquadnodes, mTransition, intdim)
end

"""
Solve the dynamic programming problem.

##### Parameters
- `p::DDP` : Object that contains the dynamic optimization problem. See `createDDP` for details.
- `;disp::Bool(false)`: turn on/off displaying of intermeditate value function iteration steps
- `;disp_each_iter::Int(10)`: how often to display output
- `;max_iter::Int(500)`: Maximum number of iterations
- `;epsilon::Float64(1e-3)`: Value for epsilon-optimality. A lower value means the solution will be more accurate.
- `;rewardcall::Symbol(:jit)`:  Symbol specifying the handling of the reward in the VFI.
	Acceptable arguments are `:jit` if the reward matrix should not be prebuilt,
	`:pre` if it should be calculated beforehand, and `:pre_partial` if
	only part of the reward matrix that only depends on states should be calculated beforehand.
- `;monotonicity::Union{Bool,Vector{Bool}}(false)`: Option that can speed-up the solver
if optimal actions are monotonically increasing in their respective state variable.
Assumes that the action is next period's state.
- `;concavity::Union{Bool,Vector{Bool}}(false)`:  exploit the concavity of the value
function in the choice of next period's state.
- `:numquadnodes::Vector{Int}`: number of quadrature nodes to use for calculating expectations.

##### Returns

TODO

"""
function solve(p::DDP; kwargs...)

	opts = SolverOptions(p; kwargs...)
	# @unpack_SolverOptions opts
	@unpack mTransition, rewardcall, intdim, disp, numquadnodes = opts

    if rewardcall == pre || intdim == All
        mReward = rewardmatrix(p)
    elseif rewardcall == pre_partial
        mReward = rewardmatrix_partial(p)
    elseif rewardcall == jit
        mReward = nothing
    else
        error("supplied wrong rewardcall option")
    end

    if mTransition == nothing
        mTransition = Array(transitionmatrix(p, numquadnodes = numquadnodes))
    end

    if intdim == All
        meshValFun, tmeshPolFun = _solve(p, intdim, mTransition, mReward, disp)
    else
        # if typeof(p) <: DDP{NS,2} where NS
        #     # if user only supplied monotonivity/concavity for one option, extend to vector
        #     if typeof(monotonicity) == Bool; monotonicity = fill!(Vector{Bool}(undef, 2), monotonicity) end
        #     if typeof(concavity) == Bool; concavity = fill!(Vector{Bool}(undef, 2), concavity) end
        # end

        meshValFun, tmeshPolFun = _solve(p, mTransition, mReward, opts)
    end

    createsolution(p, meshValFun, tmeshPolFun)
end

function display_iter(iter::Int, disp_each_iter::Int, max_diff::Float64)
	if mod(iter, disp_each_iter)==0 || iter == 1
		println(" Iteration = ", iter, " Sup Diff = ", round(max_diff, sigdigits=5))
	end
end


text:
# This file defines a hierarchy of typed expressions

##########################################################################
#
#   Typed expressions
#
##########################################################################

# modes: the role/functionality of an expression

abstract TMode
type ScalarMode <: TMode end
type EWiseMode{D} <: TMode end
type ReducMode <: TMode end
type ColwiseReducMode <: TMode end
type RowwiseReducMode <: TMode end

# types

abstract TExpr

# expressions that can be an argument for element-wise computation
abstract TEWise <: TExpr

abstract TGeneralVar <: TEWise

# expressions that is sured to be a scalar
abstract TScalar <: TEWise

type TEmpty <: TExpr end

# simple expressions

type TNum{T<:Number} <: TScalar
    val::T
end
ju_expr(tx::TNum) = tx.val

==(a::TNum, b::TNum) = (a.val == b.val)
!=(a::TNum, b::TNum) = !(a == b)

type TScalarVar <: TScalar
    name::Symbol
end
ju_expr(tx::TScalarVar) = tx.name

scalar(x::Number) = x

==(a::TScalarVar, b::TScalarVar) = (a.name == b.name)
!=(a::TScalarVar, b::TScalarVar) = !(a == b)

type TVar <: TGeneralVar
    name::Symbol
end
ju_expr(tx::TVar) = tx.name

==(a::TVar, b::TVar) = (a.name == b.name)
!=(a::TVar, b::TVar) = !(a == b)

# qualified variable, e.g. abc.x
type TQVar <: TGeneralVar
    form::Expr
end
ju_expr(tx::TQVar) = tx.form

==(a::TQVar, b::TQVar) = (a.form == b.form)
!=(a::TQVar, b::TQVar) = !(a == b)

type TGeneralScalar <: TScalar
    form::Expr
end
ju_expr(tx::TGeneralScalar) = tx.form

as_scalar(a::TScalar) = a
as_scalar(a::TVar) = TScalarVar(a.name)
as_scalar(a::TQVar) = TGeneralScalar(a.form)
as_scalar(a::TGeneralVar) = TGeneralScalar(ju_expr(a))
as_scalar(a::TExpr) = throw(DeError("Expression of type $(typeof(a)) can not be tagged as scalar."))

# references

abstract TRef <: TEWise

@compat typealias TIndex Union{Int,Symbol}
abstract TRange

type TColon <: TRange end
@compat type TInterval <: TRange
    first::TIndex
    last::Union{TIndex,Void}
end

==(a::TInterval, b::TInterval) = (a.first == b.first) && (a.last == b.last)
!=(a::TInterval, b::TInterval) = !(a == b)

# integer index
type TScalarRef1 <: TGeneralVar # or TGeneralScalar??
    host::TGeneralVar
    i::Any
end
ju_expr(tx::TScalarRef1) = :( $(ju_expr(tx.host))[$(tx.i)] )

==(a::TScalarRef1, b::TScalarRef1) = (a.host == b.host) && (a.i == b.i)
!=(a::TScalarRef1, b::TScalarRef1) = !(a == b)

type TGeneralRef1 <: TRef
    host::TGeneralVar
    i::Any
end
ju_expr(tx::TGeneralRef1) = :( $(ju_expr(tx.host))[$(tx.i)] )
as_scalar(a::TGeneralRef1) = TGeneralScalar(ju_expr(a))

==(a::TGeneralRef1, b::TGeneralRef1) = (a.host == b.host) && (a.i == b.i)
!=(a::TGeneralRef1, b::TGeneralRef1) = !(a == b)

type TGeneralRef2 <: TGeneralVar
    host::TGeneralVar
    i::Any
    j::Any
end
ju_expr(tx::TGeneralRef2) = :( $(ju_expr(tx.host))[$(tx.i), $(tx.j)] )
as_scalar(a::TGeneralRef2) = TGeneralScalar(ju_expr(a))

==(a::TGeneralRef2, b::TGeneralRef2) = (a.host == b.host) && (a.i == b.i) && (a.j == b.j)
!=(a::TGeneralRef2, b::TGeneralRef2) = !(a == b)

type TRef1D <: TRef
    host::TGeneralVar
    rgn::TRange
end

==(a::TRef1D, b::TRef1D) = (a.host == b.host) && (a.rgn == b.rgn)
!=(a::TRef1D, b::TRef1D) = !(a == b)

type TRefCol <: TRef
    host::TGeneralVar
    rrgn::TRange
    icol::TIndex
end

==(a::TRefCol, b::TRefCol) = (a.host == b.host) && (a.rrgn == b.rrgn) && (a.icol == b.icol)
!=(a::TRefCol, b::TRefCol) = !(a == b)

type TRefRow <: TRef
    host::TGeneralVar
    irow::TIndex
    crgn::TRange
end

==(a::TRefRow, b::TRefRow) = (a.host == b.host) && (a.irow == b.irow) && (a.crgn == b.crgn)
!=(a::TRefRow, b::TRefRow) = !(a == b)

type TRef2D <: TRef
    host::TGeneralVar
    rrgn::TRange
    crgn::TRange
end

==(a::TRef2D, b::TRef2D) = (a.host == b.host) && (a.rrgn == b.rrgn) && (a.crgn == b.crgn)
!=(a::TRef2D, b::TRef2D) = !(a == b)

# function calls

@compat type TMap <: TEWise
    fun::Symbol
    args::Tuple{Vararg{TEWise}}
    mode::TMode
    deps::Union{Array{TExpr},Void}
end

function ju_expr(tx::TMap)
    ju_args = [ju_expr(a) for a in tx.args]
    Expr(:call, tx.fun, ju_args...)
end
as_scalar(tx::TMap) = TGeneralScalar(ju_expr(tx))

@compat type TReduc <: TExpr
    fun::Symbol
    args::Tuple{Vararg{TEWise}}
    arg_mode::TMode
    deps::Union{Array{TExpr},Void}
end

@compat type TColwiseReduc <: TExpr
    fun::Symbol
    args::Tuple{Vararg{TEWise}}
    deps::Union{Array{TExpr},Void}
end

@compat type TRowwiseReduc <: TExpr
    fun::Symbol
    args::Tuple{Vararg{TEWise}}
    deps::Union{Array{TExpr},Void}
end

@compat typealias TFunCall Union{TMap,TReduc,TColwiseReduc,TRowwiseReduc}

function ==(a::TFunCall, b::TFunCall)
    na = length(a.args)
    if a.fun == b.fun && typeof(a) == typeof(b) && na == length(b.args)
        for i = 1 : na
            if !(a.args[i] == b.args[i])
                return false
            end
        end
        return true
    end
    return false
end

!=(a::TFunCall, b::TFunCall) = !(a == b)

# others

@compat typealias TRValue Union{TEWise,TFunCall}
@compat typealias TLValue Union{TGeneralVar,TRef}

type TAssign{Lhs<:TLValue, Rhs<:TRValue} <: TExpr
    lhs::Lhs
    rhs::Rhs
    mode::TMode
end

==(a::TAssign, b::TAssign) = (a.lhs == b.lhs) && (a.rhs == b.rhs)
!=(a::TAssign, b::TAssign) = !(a == b)

is_trivial_assignment(ex::TAssign) = isa(ex.lhs, TGeneralVar) &&
    (isa(ex.rhs, TGeneralVar) || isa(ex.rhs, TNum))

type TBlock <: TExpr
    stmts::Array{TAssign}
    TBlock() = new(TAssign[])
end

==(a::TBlock, b::TBlock) = length(a.stmts) == length(b.stmts) && all(a.stmts .== b.stmts)
!=(a::TBlock, b::TBlock) = !(a == b)

##########################################################################
#
#   Determining the expression mode
#
#   i.e. what kind of things that the whole expression wants to do
#   (ewise, reduc, etc)
#
##########################################################################

tmode_num{D}(::EWiseMode{D}) = D

tmode(ex::TScalar) = ScalarMode()
tmode(ex::TScalarRef1) = EWiseMode{0}()
tmode(ex::TGeneralVar) = EWiseMode{0}()
tmode(ex::TGeneralRef1) = EWiseMode{0}()

tmode(ex::TRef1D) = EWiseMode{1}()
tmode(ex::TRef2D) = EWiseMode{2}()
tmode(ex::TRefCol) = EWiseMode{1}()
tmode(ex::TRefRow) = EWiseMode{1}()

tmode(ex::TMap) = ex.mode
tmode(ex::TReduc) = isa(ex.arg_mode, ScalarMode) ? ScalarMode() : ReducMode()
tmode(ex::TColwiseReduc) = ColwiseReducMode()
tmode(ex::TRowwiseReduc) = RowwiseReducMode()

tmode(ex::TAssign) = ex.mode

promote_ewise_tmode(m::TMode) = m
promote_ewise_tmode(m1::ScalarMode, m2::ScalarMode) = ScalarMode()
promote_ewise_tmode(m1::ScalarMode, m2::EWiseMode) = EWiseMode{tmode_num(m2)}()
promote_ewise_tmode(m1::EWiseMode, m2::ScalarMode) = EWiseMode{tmode_num(m1)}()

function promote_ewise_tmode(m1::EWiseMode, m2::EWiseMode)
    d1 = tmode_num(m1)
    d2 = tmode_num(m2)

    d1 == d2 ? EWiseMode{d1}() :
    d1 == 0 ? EWiseMode{d2}() :
    d2 == 0 ? EWiseMode{d1}() :
    throw(DeError("Incompatible ewise mode."))
end

promote_ewise_tmode(m1::TMode, m2::TMode, m3::TMode) = promote_ewise_tmode(promote_ewise_tmode(m1, m2), m3)

##########################################################################
#
#   functions to construct simple components of a typed expression
#
##########################################################################

type DeError <: Exception
    msg::ASCIIString
end

# number literals

tnum(x::Number) = TNum{typeof(x)}(x)

# variable

tvar(s::Symbol) = TVar(s)
tscalarvar(s::Symbol) = TScalarVar(s)

# qualified variable

function is_valid_tqvar(ex::Expr)
    @assert ex.head == :(.) && length(ex.args) == 2

    a1 = ex.args[1]
    a2 = ex.args[2]

    (isa(a1, Symbol) || is_valid_tqvar(a1)) &&
    (   isa(a2, QuoteNode) ||
        (a2.head == (:quote) && isa(a2.args[1], Symbol))
    )
end

function tqvar(ex::Expr)
    if !is_valid_tqvar(ex)
        throw(DeError("$ex is not a valid form for variable."))
    end
    TQVar(ex)
end

# reference expressions

function check_ref_validity(ex::Expr, c::Bool)
    if !c
        throw(DeError("The ref-expression $ex is not supported yet."))
    end
end

tref_arg(i::Any) = i
tref_arg(i::Int) = i
tref_arg(i::Symbol) = i == :(:) ? TColon() : i

function tref_arg(ex::Expr)
    if ex.head == :(:)
        check_ref_validity(ex, length(ex.args) == 2)

        a1 = ex.args[1]
        a2 = ex.args[2]

        if isa(a1, Int) || (isa(a1, Symbol) && a1 != :(:))
            first = a1
        else
            check_ref_validity(ex, false)
        end

        if isa(a2, Int)
            last = a2
        elseif isa(a2, Symbol)
            last = (a2 == :(:) || a2 == symbol("end")) ? nothing : a2
        else
            check_ref_validity(ex, false)
        end

        TInterval(first, last)
    else
        ex
    end
end

tref(x::TGeneralVar, i::Any) = TGeneralRef1(x, i)
tref(x::TGeneralVar, i::Int) = TScalarRef1(x, i)
tref(x::TGeneralVar, r::TRange) = TRef1D(x, r)

tref(x::TGeneralVar, i::Any, j::Any) = TGeneralRef2(x, i, j)
tref(x::TGeneralVar, i::TIndex, c::TRange) = TRefRow(x, i, c)
tref(x::TGeneralVar, r::TRange, j::TIndex) = TRefCol(x, r, j)
tref(x::TGeneralVar, r::TRange, c::TRange) = TRef2D(x, r, c)

function tref(ex::Expr)
    @assert ex.head == :(ref)

    na = length(ex.args)
    check_ref_validity(ex, na == 2 || na == 3)

    h = ex.args[1]
    if isa(h, Symbol)
        h = tvar(h)
    elseif isa(h, Expr) && h.head == :(.)
        h = tqvar(h)
    else
        check_ref_validity(ex, false)
    end

    if na == 2
        a1 = ex.args[2]
        tref(h, tref_arg(a1))

    else
        a1 = ex.args[2]
        a2 = ex.args[3]
        tref(h, tref_arg(a1), tref_arg(a2))
    end
end

##########################################################################
#
#   function call recognition
#
##########################################################################

is_ewise_call(f::Symbol, N::Int) = isa(get_op_kind(TCallSig{f, N}()), EWiseOp)
is_reduc_call(f::Symbol, N::Int) = isa(get_op_kind(TCallSig{f, N}()), ReducOp)

function check_funcall_args(args::TExpr...)
    na = length(args)
    deps = nothing

    pargs = Array(TExpr, na)

    for i = 1 : na
        a = args[i]
        if isa(a, TEWise)
            pargs[i] = a
        elseif isa(a, TReduc) || isa(a, TColwiseReduc) || isa(a, TRowwiseReduc)
            dep_sym = gensym("dep")
            pargs[i] = isa(a, TReduc) ? TScalarVar(dep_sym) : TVar(dep_sym)
            if deps == nothing
                deps = TExpr[]
            end
            push!(deps, tassign(TVar(dep_sym), a))
        else
            throw(DeError("Arguments in unsupported form."))
        end
    end
    (tuple(pargs...), deps)
end

# due to intricate semantics of partial reduction in Julia syntax
# we here hand-coded the specific form to be recognized as partial reduction

function recognize_partial_reduction(f::Symbol, a::TExpr...)
    if f == (:sum) || f == (:mean) || f == (:maximum) || f == (:minimum)
        if length(a) == 2 && isa(a[2], TNum{Int})
            fargs, deps = check_funcall_args(a[1])
            dim = a[2].val
            if dim == 1
                TColwiseReduc(f, fargs, deps)
            elseif dim == 2
                TRowwiseReduc(f, fargs, deps)
            else
                throw(DeError("Devectorize supports either colwise or rowwise reduction."))
            end
        end
    end
end

function tcall(f::Symbol, args)
    n = length(args)

    if f == (:scalar) && n == 1
        as_scalar(args[1])

    elseif is_ewise_call(f, n)
        fargs, deps = check_funcall_args(args...)
        mode = promote_ewise_tmode([tmode(a) for a in fargs]...)
        TMap(f, fargs, mode, deps)

    elseif is_reduc_call(f, n)
        fargs, deps = check_funcall_args(args...)
        arg_mode = promote_ewise_tmode([tmode(a) for a in fargs]...)
        TReduc(f, fargs, arg_mode, deps)

    elseif f == :(*)
        throw(DeError("Devectorize does not support *, please use .* to express element-wise multiplication."))

    elseif f == :(/)
        throw(DeError("Devectorize does not support /, please use ./ to express element-wise division."))

    elseif f == :(^)
        throw(DeError("Devectorize does not support ^, please use .^ to express element-wise power."))

    else
        ex = recognize_partial_reduction(f, args...)
        if ex == nothing
            throw(DeError("Unrecognized function $f with $n arguments (in Devectorize)"))
        end
        return ex
    end
end

function tcall(ex::Expr)
    @assert ex.head == :(call)
    fsym = ex.args[1]
    if !isa(fsym, Symbol)
        throw(DeError("call-expressions with non-symbol function name: $fsym"))
    end
    tcall(fsym, map(texpr, ex.args[2:end]))
end

function tcomparison(ex::Expr)
    @assert ex.head == :(comparison)
    opsym = ex.args[2]
    tcall(opsym, map(texpr, [ex.args[1], ex.args[3]]))
end

##########################################################################
#
#   assignment construction
#
##########################################################################

function decide_assign_tmode(lhs::TExpr, rhs::TExpr)

    # TScalarSym can only be created internally, which would never
    # be placed on the left hand side
    @assert !isa(lhs, TScalarVar)

    if isa(lhs, TGeneralVar)
        mode = tmode(rhs)

    elseif isa(lhs, TRef)
        @assert isa(lhs, TEWise)
        rmode = tmode(rhs)
        if isa(rmode, EWiseMode) || isa(rmode, ScalarMode)
            mode = promote_ewise_tmode(tmode(lhs), rmode)
        elseif isa(rmode, ColwiseReducMode) || isa(rmode, RowwiseReducMode)
            mode = rmode
        else
            throw(DeError("Incompatible lhs and rhs."))
        end

    else
        throw(DeError("Incompatible modes between lhs and rhs."))
    end
end

function tassign(lhs::TExpr, rhs::TExpr)
    if isa(rhs, TAssign)
        throw(DeError("chained assignment is not supported by Devectorize"))
    end
    mode = decide_assign_tmode(lhs, rhs)
    TAssign{typeof(lhs),typeof(rhs)}(lhs, rhs, mode)
end

function topassign(aop::Symbol, lhs::TExpr, rhs::TExpr)
    op = extract_assign_op(TFun{aop}())
    new_rhs = tcall(op, [lhs, rhs])
    tassign(lhs, new_rhs)
end

function tassign(ex::Expr)
    @assert ex.head == :(=)
    @assert length(ex.args) == 2

    lhs = ex.args[1]
    rhs = ex.args[2]
    tassign(texpr(lhs), texpr(rhs))
end

function topassign(ex::Expr)
    topassign(ex.head, texpr(ex.args[1]), texpr(ex.args[2]))
end

##########################################################################
#
#   block construction
#
##########################################################################

function tblock(blk_ex::Expr)

    @assert blk_ex.head == :(block)

    blk = TBlock()
    for e in blk_ex.args
        if isa(e, LineNumberNode) || e.head == (:line)
            continue
        end
        if !(e.head == :(=) || is_opassign(e.head))
            throw(DeError("Each statement in a block must be an assignment or op-assignment."))
        end
        push!(blk.stmts, texpr(e))
    end
    blk
end

##########################################################################
#
#   AST construction:  Expr ==> TExpr
#
##########################################################################

texpr(x::Number) = tnum(x)
texpr(x::Symbol) = tvar(x)

is_empty_tuple(ex::Expr) = ex.head == :(tuple) && isempty(ex.args)

texpr(ex::Expr) =
    ex.head == :(.) ? tqvar(ex) :
    ex.head == :(call) ? tcall(ex) :
    ex.head == :(comparison) ? tcomparison(ex) :
    ex.head == :(ref) ? tref(ex) :
    ex.head == :(=) ? tassign(ex) :
    ex.head == :(block) ? tblock(ex) :
    is_empty_tuple(ex) ? TEmpty() :
    is_opassign(ex.head) ? topassign(ex) :
    throw(DeError("Unrecognized expression: $ex"))


text:
using Test, LibCEED, LinearAlgebra, StaticArrays

@testset "LibCEED Development Tests" begin end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text: Dict{String, Any}("true" => Dict{String, Any}("value" => "1", "type" => "integer"), "false" => Dict{String, Any}("value" => "false", "type" => "bool"), "nan" => Dict{String, Any}("value" => "ceci n'est pas un nombre", "type" => "string"), "inf" => Dict{String, Any}("value" => "100000000", "type" => "integer"))


text:
function gradientdescent(f, α, ρ, c, current...; maxiter=200)
    ok = true
    i = 1
    while i < maxiter
        v, back = pullback(f, current...)
        grad = back(1.0)
        println(i)
        println("Grad: ", grad)
        println("Before: ", current)
        current = backtrack(f, α, ρ, c, current, grad, v)
        println("After: ", current)
        i += 1
    end
    return current
end

function backtrack(f, α, ρ, c, start, grad, v)
    ngrad = sum(x -> norm(x) ^ 2, grad)
    current = start .- α .* grad
    v′ = f(current...)
    while v′ > v - c * α * ngrad
        α *= ρ
        current = start .- α .* grad
        v′ = f(current...)
    end
    return current
end

function make_fg(f, unpack)
    return function(F, G, x)
        v, back = pullback(x) do x
            f(unpack(x)...)
        end
        if G != nothing
            G .= back(1.0)[1]
        end
        if F != nothing
            return v
        end
        return nothing
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module MDCDL # Multi-Dimensional Convolutional Dictionary Learning

using LinearAlgebra
using ImageFiltering.Algorithm: Alg, FIR, FFT
using ComputationalResources: AbstractResource, CPU1

import Base: promote_rule, eltype, ndims, similar, length
import Random: rand, rand!

export rand, rand!

export Shapes, Optimizers, SparseCoders

export  PolyphaseVector, FilterBank, PolyphaseFB, AbstractNsolt,
        Cnsolt, CnsoltTypeI, CnsoltTypeII, Rnsolt, RnsoltTypeI, RnsoltTypeII,
        ParallelFilters, Multiscale, AbstractOperator, TransformSystems

export  istype1, istype2, permdctmtx, cdftmtx, analyze, synthesize,
        adjoint_synthesize, upsample, downsample, serialize, deserialize,
        analysisbank, decimations, orders, nchannels, nstages, kernels,
        analysiskernels, synthesiskernels, kernelsize, getrotations,
        setrotations!, mdarray2polyphase, polyphase2mdarray, iht,
        createTransform, loadfb, savefb, copy_params!

include("sparsecoders/SparseCoders.jl")
include("optimizers/Optimizers.jl")
using MDCDL.Optimizers: iterations

struct PolyphaseVector{T,D}
    data::AbstractMatrix{T}
    nBlocks::NTuple{D, Int}
end

promote_rule(::Type{PolyphaseVector{TA,D}}, ::Type{PolyphaseVector{TB,D}}) where {TA,TB,D} = PolyphaseVector{promote_type(TA,TB), D}

abstract type CodeBook{T,D} end
abstract type FilterBank{T,D} <: CodeBook{T,D} end
abstract type PolyphaseFB{T,D} <: FilterBank{T,D} end

eltype(::Type{CB}) where {T,D,CB<:CodeBook{T,D}} = T
ndims(::Type{CB}) where {T,D,CB<:CodeBook{T,D}} = D
ndims(cb::CodeBook) = ndims(typeof(cb))

decimations(fb::FilterBank) = fb.decimationFactor
nchannels(fb::FilterBank) = sum(fb.nChannels)
orders(fb::FilterBank) = fb.polyphaseOrder

kernels(fb::FilterBank) = fb.kernels
kernelsize(fb::FilterBank) = decimations(fb) .* (1 .+ orders(fb))

struct Multiscale
    filterbanks::Tuple
    Multiscale(fbs...) = new(fbs)
end

length(ms::Multiscale) = length(ms.filterbanks)

include("filterbanks/nsolt.jl")
include("filterbanks/parameterization.jl")
include("filterbanks/parallelFilters.jl")

include("orthonormalMatrixSystem.jl")
include("polyphaseMatrices.jl")

include("shapes.jl")
include("transformSystem.jl")

include("analysisSystem.jl")
include("synthesisSystem.jl")

include("dictionaryLearning/dictionaryLearning.jl")
include("dictionaryLearning/sparseCodingStage.jl")
include("dictionaryLearning/dictionaryUpdateStage.jl")
include("dictionaryLearning/logger.jl")

include("recipes.jl")
include("io.jl")

include("utils.jl")
include("compatibles.jl")

end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using ForwardDiff

import Rotations: jacobian, ∇rotate, ∇composition1, ∇composition2
import Rotations: kinematics, _pure_quaternion, params
import Rotations: vmat, rmult, lmult, hmat, tmat

@testset "Unit Quaternions" begin
    q1 = rand(UnitQuaternion)
    q2 = rand(UnitQuaternion)
    r = @SVector rand(3)
    ω = @SVector rand(3)

    # Constructors
    @test UnitQuaternion(1.0, 0.0, 0.0, 0.0) isa UnitQuaternion{Float64}
    @test UnitQuaternion(1.0, 0, 0, 0) isa UnitQuaternion{Float64}
    @test UnitQuaternion(1.0, 0, 0, 0) isa UnitQuaternion{Float64}
    @test UnitQuaternion(1, 0, 0, 0) isa UnitQuaternion{Int}
    @test UnitQuaternion(1.0f0, 0, 0, 0) isa UnitQuaternion{Float32}

    q = normalize(@SVector rand(4))
    q32 = SVector{4,Float32}(q)
    @test UnitQuaternion(q) isa UnitQuaternion{Float64}
    @test UnitQuaternion(q32) isa UnitQuaternion{Float32}

    r = @SVector rand(3)
    r32 = SVector{3,Float32}(r)
    @test _pure_quaternion(r) isa Quaternion{Float64}
    @test _pure_quaternion(r32) isa Quaternion{Float32}
    @test real(_pure_quaternion(r)) == 0

    @test UnitQuaternion{Float64}(1, 0, 0, 0) isa UnitQuaternion{Float64}
    @test UnitQuaternion{Float32}(1, 0, 0, 0) isa UnitQuaternion{Float32}
    @test UnitQuaternion{Float32}(1.0, 0, 0, 0) isa UnitQuaternion{Float32}
    @test UnitQuaternion{Float64}(1f0, 0, 0, 0) isa UnitQuaternion{Float64}

    # normalization
    @test UnitQuaternion(2.0, 0, 0, 0, true) == one(UnitQuaternion)
    @test UnitQuaternion(2q, true) ≈ UnitQuaternion(q)

    # Copy constructors
    q = rand(UnitQuaternion)
    @test UnitQuaternion(q) === q
    @test UnitQuaternion{Float32}(q) isa UnitQuaternion{Float32}
    UnitQuaternion{Float32}(q)

    # rand
    @test rand(UnitQuaternion) isa UnitQuaternion{Float64}
    @test rand(UnitQuaternion{Float32}) isa UnitQuaternion{Float32}

    # Test math
    @test UnitQuaternion(I) isa UnitQuaternion{Float64}

    ϕ = inv(ExponentialMap())(q1)
    @test expm(ϕ * 2) ≈ q1
    q = Rotations._pure_quaternion(ϕ)
    @test UnitQuaternion(exp(q)) ≈ q1
    @test exp(q) isa Quaternion

    q = UnitQuaternion((@SVector [1, 2, 3, 4.0]), false)
    @test 2 * q == 2 * Matrix(q)
    @test q * 2 == 2 * Matrix(q)

    # Axis-angle
    ϕ = 0.1 * @SVector [1, 0, 0]
    q = expm(ϕ)
    @test q isa UnitQuaternion
    @test logm(expm(ϕ)) ≈ ϕ
    @test expm(logm(q1)) ≈ q1
    @test rotation_angle(q) ≈ 0.1
    @test rotation_axis(q) == [1, 0, 0]

    @test norm(q1 * ExponentialMap()(ϕ)) ≈ √3
    @test q1 ⊖ q2 isa StaticVector{3}
    @test (q1 * CayleyMap()(ϕ)) ⊖ q1 ≈ ϕ

    # Test inverses
    q3 = q2 * q1
    @test q2 \ q3 ≈ q1
    @test q3 / q1 ≈ q2
    @test inv(q1) * r ≈ q1 \ r
    @test r ≈ q3 \ (q2 * q1 * r)
    @test q3 ⊖ q2 ≈ inv(CayleyMap())(q1)

    q = q1
    rhat = Rotations._pure_quaternion(r)
    @test q * r ≈ vmat() * lmult(q) * rmult(q)' * vmat()'r
    @test q * r ≈ vmat() * lmult(q) * rmult(q)' * hmat(r)
    @test q * r ≈ vmat() * lmult(q) * lmult(rhat) * tmat() * params(q)
    @test q * r ≈ vmat() * rmult(q)' * rmult(rhat) * params(q)
    @test q * r ≈ hmat()' * rmult(q)' * rmult(rhat) * params(q)

    @test rmult(params(q)) == rmult(q)
    @test lmult(params(q)) == lmult(q)
    @test hmat(r) == params(_pure_quaternion(r))

    # Test Jacobians
    @test ForwardDiff.jacobian(q -> UnitQuaternion(q, false) * r, params(q)) ≈
          ∇rotate(q, r)

    @test ForwardDiff.jacobian(q -> params(q2 * UnitQuaternion(q, false)), params(q1)) ≈
          ∇composition1(q2, q1)
    @test ForwardDiff.jacobian(q -> params(UnitQuaternion(q, false) * q1), params(q2)) ≈
          ∇composition2(q2, q1)

    b = @SVector rand(4)
    qval = params(q1)
    ForwardDiff.jacobian(q -> ∇composition1(q2, UnitQuaternion(q))'b,
        @SVector [1, 0, 0, 0.0,])
    diffcomp = ϕ -> params(q2 * CayleyMap()(ϕ))
    ∇diffcomp(ϕ) = ForwardDiff.jacobian(diffcomp, ϕ)
    @test ∇diffcomp(@SVector zeros(3)) ≈ Rotations.∇differential(q2)
    @test ForwardDiff.jacobian(ϕ -> ∇diffcomp(ϕ)'b, @SVector zeros(3)) ≈
          Rotations.∇²differential(q2, b)

    @test lmult(q) ≈ ∇composition1(q, q2)

    ϕ = @SVector zeros(3)
    @test Rotations.∇differential(q) ≈ lmult(q) * jacobian(QuatVecMap(), ϕ)
    @test Rotations.∇differential(q) ≈ lmult(q) * jacobian(ExponentialMap(), ϕ)
    @test Rotations.∇differential(q) ≈ lmult(q) * jacobian(CayleyMap(), ϕ)
    @test Rotations.∇differential(q) ≈ lmult(q) * jacobian(MRPMap(), ϕ)

    # Check ops with Float32
    ϕ = SA_F32[1, 2, 3]
    @test expm(SA_F32[1, 2, 3]) isa UnitQuaternion{Float32}

    q32 = rand(UnitQuaternion{Float32})
    @test Rotations._log_as_quat(q32) isa Quaternion{Float32}
    @test log(q32) isa SMatrix
    @test eltype(logm(q32)) == Float32
    @test expm(logm(q32)) ≈ q32

    @test normalize(q32) isa SMatrix{3,3,Float32}

    ω = @SVector rand(3)
    ω32 = Float32.(ω)
    @test Rotations.kinematics(q, ω) isa SVector{4,Float64}
    @test Rotations.kinematics(q32, ω32) isa SVector{4,Float32}
    @test Rotations.kinematics(q32, ω) isa SVector{4,Float32}
    @test Rotations.kinematics(q32, [1, 2, 3]) isa SVector{4,Float32}

    @test eltype(lmult(q32)) == Float32
    @test eltype(lmult(q)) == Float64

    @test eltype(tmat()) == Float64
    @test eltype(tmat(Int)) == Int
    @test eltype(vmat(Float32)) == Float32
    @test eltype(hmat(Float32)) == Float32
end


text:
mutable struct SwitchUpdateState{T}
    weight::Float64
    score::Float64
    noise::Float64
    prev_trace::Trace
    trace::Trace
    index::Int
    discard::ChoiceMap
    updated_retdiff::Diff
    SwitchUpdateState{T}(weight::Float64, score::Float64, noise::Float64, prev_trace::Trace) where T = new{T}(weight, score, noise, prev_trace)
end

function update_recurse_merge(prev_choices::ChoiceMap, choices::ChoiceMap)
    prev_choice_submap_iterator = get_submaps_shallow(prev_choices)
    prev_choice_value_iterator = get_values_shallow(prev_choices)
    choice_submap_iterator = get_submaps_shallow(choices)
    choice_value_iterator = get_values_shallow(choices)
    new_choices = DynamicChoiceMap()

    # Add (address, value) to new_choices from prev_choices if address does not occur in choices.
    for (address, value) in prev_choice_value_iterator
        address in keys(choice_value_iterator) && continue
        set_value!(new_choices, address, value)
    end

    # Add (address, submap) to new_choices from prev_choices if address does not occur in choices.
    # If it does, enter a recursive call to update_recurse_merge.
    for (address, node1) in prev_choice_submap_iterator
        if address in keys(choice_submap_iterator)
            node2 = get_submap(choices, address)
            node = update_recurse_merge(node1, node2)
            set_submap!(new_choices, address, node)
        else
            set_submap!(new_choices, address, node1)
        end
    end

    # Add (address, value) from choices to new_choices. This is okay because we've excluded any conflicting addresses from the prev_choices above.
    for (address, value) in choice_value_iterator
        set_value!(new_choices, address, value)
    end

    sel, _ = zip(prev_choice_submap_iterator...)
    comp = complement(select(sel...))
    for (address, node) in get_submaps_shallow(get_selected(choices, comp))
        set_submap!(new_choices, address, node)
    end
    return new_choices
end

@doc(
"""
update_recurse_merge(prev_choices::ChoiceMap, choices::ChoiceMap)

Returns choices that are in constraints, merged with all choices in the previous trace that do not have the same address as some choice in the constraints."
""", update_recurse_merge)

function update_discard(prev_choices::ChoiceMap, choices::ChoiceMap, new_choices::ChoiceMap)
    discard = choicemap()
    for (k, v) in get_submaps_shallow(prev_choices)
        new_submap = get_submap(new_choices, k)
        choices_submap = get_submap(choices, k)
        sub_discard = update_discard(v, choices_submap, new_submap)
        set_submap!(discard, k, sub_discard)
    end
    for (k, v) in get_values_shallow(prev_choices)
        if (!has_value(new_choices, k) || has_value(choices, k))
            set_value!(discard, k, v)
        end
    end
    discard
end

@doc(
"""
update_discard(prev_choices::ChoiceMap, choices::ChoiceMap, new_choices::ChoiceMap)

Returns choices from previous trace that:
   1. have an address which does not appear in the new trace.
   2. have an address which does appear in the constraints.
""", update_discard)

@inline update_discard(prev_trace::Trace, choices::ChoiceMap, new_trace::Trace) = update_discard(get_choices(prev_trace), choices, get_choices(new_trace))

function process!(gen_fn::Switch{C, N, K, T},
        index::Int,
        index_argdiff::UnknownChange,
        args::Tuple,
        kernel_argdiffs::Tuple,
        choices::ChoiceMap,
        state::SwitchUpdateState{T}) where {C, N, K, T, DV}

    # Generate new trace.
    merged = update_recurse_merge(get_choices(state.prev_trace), choices)
    branch_fn = getfield(gen_fn.branches, index)
    new_trace, weight = generate(branch_fn, args, merged)
    weight -= get_score(state.prev_trace)
    state.discard = update_discard(state.prev_trace, choices, new_trace)

    # Set state.
    state.index = index
    state.weight = weight
    state.noise = project(new_trace, EmptySelection()) - project(state.prev_trace, EmptySelection())
    state.score = get_score(new_trace)
    state.trace = new_trace
    state.updated_retdiff = UnknownChange()
end

function process!(gen_fn::Switch{C, N, K, T},
        index::Int,
        index_argdiff::NoChange, # TODO: Diffed wrapper?
        args::Tuple,
        kernel_argdiffs::Tuple,
        choices::ChoiceMap,
        state::SwitchUpdateState{T}) where {C, N, K, T}

    # Update trace.
    new_trace, weight, retdiff, discard = update(getfield(state.prev_trace, :branch), args, kernel_argdiffs, choices)

    # Set state.
    state.index = index
    state.weight = weight
    state.noise = project(new_trace, EmptySelection()) - project(state.prev_trace, EmptySelection())
    state.score = get_score(new_trace)
    state.trace = new_trace
    state.updated_retdiff = retdiff
    state.discard = discard
end

@inline process!(gen_fn::Switch{C, N, K, T}, index::C, index_argdiff::Diff, args::Tuple, kernel_argdiffs::Tuple, choices::ChoiceMap, state::SwitchUpdateState{T}) where {C, N, K, T} = process!(gen_fn, getindex(gen_fn.cases, index), index_argdiff, args, kernel_argdiffs, choices, state)

function update(trace::SwitchTrace{T},
        args::Tuple,
        argdiffs::Tuple,
        choices::ChoiceMap) where T
    gen_fn = trace.gen_fn
    index, index_argdiff = args[1], argdiffs[1]
    state = SwitchUpdateState{T}(0.0, 0.0, 0.0, trace)
    process!(gen_fn, index, index_argdiff, args[2 : end], argdiffs[2 : end], choices, state)
    return SwitchTrace(gen_fn, state.index, state.trace, get_retval(state.trace), args, state.score, state.noise), state.weight, state.updated_retdiff, state.discard
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
"""
    solve(problem::MatrixFormalism, gradient; ninterval = 500)

Solve for magnetization using Matrix Formalism.
"""
function solve(problem::MatrixFormalism, gradient; ninterval = 500)
    (; model, matrices, lap_eig) = problem
    (; γ) = model
    (; M) = matrices

    # Laplace eigenmodes
    λ = lap_eig.values
    ϕ = lap_eig.funcs
    Ax = lap_eig.moments
    T = lap_eig.massrelax
    neig = length(λ)

    ρ = initial_conditions(model)

    # Project initial spin density onto Laplace eigenfunction basis
    ν = ϕ' * (M * ρ)

    # Laplace operator in Laplace eigenfunction basis
    L = diagm(λ)

    # Bloch-Torrey operator in Laplace eigenfunction basis for given gradient
    function K!(K, g⃗)
        @. K = L + T + im * γ * (g⃗[1] * Ax[1] + g⃗[2] * Ax[2] + g⃗[3] * Ax[3])
    end

    K = zeros(eltype(ρ), neig, neig)

    if isa(gradient, ScalarGradient) && isconstant(gradient.profile)
        # The gradient is constant on each interval
        t = intervals(gradient.profile)
        for i = 1:length(t)-1
            δᵢ = t[i+1] - t[i]
            tᵢ = (t[i+1] + t[i]) / 2
            K!(K, gradient(tᵢ))
            expmv!(-δᵢ, K, ν)
        end
    else
        t = LinRange(0, echotime(gradient), ninterval + 1)
        for i = 1:ninterval
            δᵢ = t[i+1] - t[i]
            g⃗ᵢ = (gradient(t[i+1]) + gradient(t[i])) / 2
            K!(K, g⃗ᵢ)
            expmv!(-δᵢ, K, ν)
        end
    end

    # Final magnetization coefficients in finite element nodal basis
    ϕ * ν
end


text:
const timeAttributes = ["age", "tenure", "time in node"]
const histAttributes = ["had transition", "started as", "was"]
const relationEntries = Dict{Function, String}( Base.:(==) => "IS",
    Base.:(!=) => "IS NOT", Base.:∈ => "IN", Base.:∉ => "NOT IN",
    Base.:> => ">", Base.:(>=) => ">=", Base.:< => "<", Base.:(<=) => "<=" )

function conditionToSQLite( condition::MPcondition, mpSim::MPsim )

    sqlValue = ""

    if condition.value isa Real
        sqlValue = string( condition.value )
    elseif condition.value isa String
        sqlValue = string( "'", condition.value, "'" )
    else
        sqlValue = string( "( '", join( condition.value, "', '" ), "' )" )
    end  # if condition.value isa Real

    sqlCond = ""

    if condition.attribute == "had transition"
        sqlCond = string( " AND transition ",
            relationEntries[condition.operator], " ", sqlValue )
    elseif condition.attribute == "started as"
        sqlCond = string( " AND ", mpSim.sNode, " IS NULL AND ", mpSim.tNode,
            " ", relationEntries[condition.operator], " ", sqlValue )
    elseif condition.attribute == "was"
        sqlCond = string( " AND ", mpSim.sNode, " ",
            relationEntries[condition.operator], " ", sqlValue )
    else  # Regular conditions.
        sqlCond = string( "`", condition.attribute, "` ",
            relationEntries[condition.operator], " ", sqlValue )
    end  # if condition.attribute == "had transition"

    return sqlCond

end  # conditionToSQLite( condition, mpSim )


text:
using Plots;
gr();
using LaTeXStrings

function velocity_verlet(x₀, v₀, a, τ, N, α)
    t = range(0, stop = N * τ, length = N + 1)
    x = Float64[x₀]
    v = Float64[v₀]

    for n = 1:N
        push!(x, x[n] + τ * v[n] + τ^2 * a(x[n], α) / 2)
        push!(v, v[n] + τ * (a(x[n], α) + a(x[n+1], α)) / 2)
    end

    return t, x, v
end

function acceleration(θ, α)
    return 9.81 * sin(θ) - (α / 2) * tan(θ / 2) * sec(θ / 2)
end

θs = range(0, stop = 2π, length = 1000)
xs = cos.(θs);
ys = sin.(θs);
function plot_system(x, y, α)
    global xs
    global ys
    scatter(
        xs,
        ys,
        title = (latexstring("\\alpha = $(α)")),
        aspect_ratio = true,
        marker = (:black, 1),
        legend = false,
        xaxis = (L"x"),
        yaxis = (L"y"),
        dpi = 144,
    )
    scatter!([0], [-1], marker = (:red, 8))
    scatter!([x], [y], marker = (:red, 8))
end

function animate_system(t, θ, α)
    interval = 500
    θ = θ[1:interval:end]
    t = collect(t)[1:interval:end]
    x = sin.(θ)
    y = cos.(θ)

    i_last = nothing
    anim = @animate for i ∈ 1:length(x)
        if i > length(x) / 2 && round(θ[i], digits = 2) ≈ round(θ₀, digits = 2)
            i_last = i
            break
        end
        plot_system(x[i], y[i], α)
    end

    fps = i_last / t[i_last]
    gif(anim, savepath * "animation_$(α).gif", fps = fps)
end

function plot_θ(t, θ, α)
    plot(
        t,
        θ,
        dpi = 144,
        title = (latexstring("\\alpha = $(α)")),
        legend = false,
        xaxis = (L"t"),
        yaxis = (
            L"\theta",
            (minimum([0.9 * minimum(θ), 1.1 * minimum(θ)]), 1.1 * maximum(θ)),
        ),
    )
    savefig(savepath * "oscillation_$(α)")
end

function plot_system_(x, y)
    global xs
    global ys
    scatter(
        xs,
        ys,
        aspect_ratio = true,
        marker = (:black, 1),
        legend = false,
        xaxis = (L"x"),
        yaxis = (L"y"),
        dpi = 144,
    )
    scatter!(
        xs[floor(Int, length(xs) / 8):floor(Int, length(xs) / 4)] .* 0.1,
        ys[floor(Int, length(xs) / 8):floor(Int, length(xs) / 4)] .* 0.1,
        marker = (:black, 1),
    )
    x1 = range(0, stop = x, length = 2)
    y1 = range(0, stop = y, length = 2)
    plot!(x1, y1, line = (:black, 1))
    x2 = range(0, stop = 1e-3, length = 2)
    y2 = range(0, stop = 1, length = 2)
    plot!(x2, y2, line = (:black, 1))
    scatter!([0], [-1], marker = (:red, 8))
    scatter!([x], [y], marker = (:red, 8))
    annotate!((0.08, 0.2, L"\theta"))
    annotate!((x + 0.1, y + 0.1, L"q"))
    annotate!((0, -0.85, L"q'"))
    annotate!((-0.1, 0.5, L"r"))
    savefig(savepath * "system.png")
end

# %%

savepath = "../images/2020-04-22-charged-particles-ring/"
plot_system_(sin(π / 4), cos(π / 4))

θ₀ = π / 4
ω₀ = 0
τ = 1e-4
N = 10^5
for α in [0.1, 1, 10, 100, 1000]
    t, θ, ω = velocity_verlet(θ₀, ω₀, acceleration, τ, N, α)
    animate_system(t, θ, α)
    plot_θ(t, θ, α)
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractSimulatorImpl

export AbstractSimulator

abstract AbstractSimulator

end


text:
# Pooling is so similar, we abstract over meanpooling and maxpooling, simply replacing
# the inner loop operation and a few initialization parameters.
for name in (:max, :mean)
    @eval function $((Symbol("$(name)pool_direct!")))(
                    y::AbstractArray{T, 5}, x::AbstractArray{T, 5},
                    pdims::PoolDims; alpha::T=T(1), beta::T=T(0)) where T
        $((Symbol("$(name)pool_direct!")))(
            y, x, pdims,
            Val(kernel_size(pdims)), Val(channels_out(pdims)),
            Val(padding(pdims)), Val(dilation(pdims)), Val(stride(pdims));
            alpha, beta)
        return y
    end

    @eval function $((Symbol("$(name)pool_direct!")))(
        y::AbstractArray{T,5}, x::AbstractArray{T,5},
        pdims::PoolDims,
        # kernel size, channels out, padding, dilation, stride
        ::Val{K}, ::Val{C}, ::Val{P}, ::Val{D}, ::Val{S};
        alpha::T=T(1), beta::T=T(0),
    ) where {T, K, C, P, D, S}
        @assert beta == T(0) "beta not supported yet"
        check_dims(size(x), size(y), pdims)

        width, height, depth = input_size(pdims)
        kernel_w, kernel_h, kernel_d = K
        pad_w_lo, _, pad_h_lo, _, pad_d_lo, _ = P
        dil_w, dil_h, dil_d = D
        stride_w, stride_h, stride_d = S

        # We use calc_padding_regions to split outselves up into separate regions that may or
        # may not need to worry about padding:
        padded_regions, central_region = calc_padding_regions(pdims)

        # A helper function to project from output (w, h) to input (input_w, input_h)
        @inline project(idx, stride, pad) = (idx - 1) * stride - pad + 1

        # If we're doing mean pooling, we represent division by kernel size by rolling it
        # into the `alpha` multiplier.
        if $(name == :mean)
            alpha = alpha / prod(K)
        end

        # Each loop, we initialize `m` to something, set that here.
        m_init = if $(name == :max)
            T <: AbstractFloat ? nextfloat(typemin(T)) : typemin(T)
        elseif $(name == :mean)
            T(0)
        else
            error("Unimplemented codegen path")
        end

        # Start with the central region
        w_region, h_region, d_region = central_region

        @inbounds for batch_idx in 1:size(x, 5), c in 1:C
            for d in d_region
            pd = project(d, stride_d, pad_d_lo)
            for h in h_region
            ph = project(h, stride_h, pad_h_lo)
            for w in w_region
            pw = project(w, stride_w, pad_w_lo)
            m = m_init

            for kd in 1:kernel_d,
                kh in 1:kernel_h,
                kw in 1:kernel_w

                input_kd = pd + (kd - 1) * dil_d
                input_kh = ph + (kh - 1) * dil_h
                input_kw = pw + (kw - 1) * dil_w

                # This conditional will be optimized away at compile time
                if $(name == :max)
                    xv = x[input_kw, input_kh, input_kd, c, batch_idx]
                    if xv > m
                        m = xv
                    end
                elseif $(name == :mean)
                    m += x[input_kw, input_kh, input_kd, c, batch_idx]
                else
                    error("Unimplemented codegen path")
                end
            end

            y[w, h, d, c, batch_idx] = alpha * m # + beta * y[w, h, d, c, batch_idx]
            end
            end
            end
        end

        # Next, the padded regions
        @inbounds for (w_region, h_region, d_region) in padded_regions
            for batch_idx in 1:size(x, 5), c in 1:C
                for d in d_region
                pd = project(d, stride_d, pad_d_lo)
                for h in h_region
                ph = project(h, stride_h, pad_h_lo)
                for w in w_region
                pw = project(w, stride_w, pad_w_lo)
                m = m_init

                for kd in 1:kernel_d
                    input_kd = pd + (kd - 1) * dil_d
                    if input_kd <= 0 || input_kd > depth
                        # add here condition for handling options for paded value handling
                        continue
                    end

                    for kh in 1:kernel_h
                        input_kh = ph + (kh - 1) * dil_h
                        if input_kh <= 0 || input_kh > height
                            # add here condition for handling options for paded value handling
                            continue
                        end

                        for kw in 1:kernel_w
                            input_kw = pw + (kw - 1) * dil_w
                            if input_kw <= 0 || input_kw > width
                                # add here condition for handling options for paded value handling
                                continue
                            end

                            if $(name == :max)
                                xv = x[input_kw, input_kh, input_kd, c, batch_idx]
                                if xv > m
                                    m = xv
                                end
                            elseif $(name == :mean)
                                m += x[input_kw, input_kh, input_kd, c, batch_idx]
                            else
                                error("Unimplemented codegen path")
                            end
                        end
                    end
                end
                y[w, h, d, c, batch_idx] = alpha * m # + beta * y[w, h, d, c, batch_idx]
                end
                end
                end
            end
        end

        return y
    end

    @eval function $((Symbol("∇$(name)pool_direct!")))(
                    dx::AbstractArray{T,5}, dy::AbstractArray{T,5},
                    y::AbstractArray{T,5}, x::AbstractArray{T,5},
                    pdims::PoolDims; kwargs...) where T
        $((Symbol("∇$(name)pool_direct!")))(
            dx, dy, y, x, pdims, Val(kernel_size(pdims)); kwargs...)
        return dx
    end

    # Same story for gradients, and although this is very similar to the forward pass,
    # it's unfortunately different enough that I think we need a separate function.  :(
    @eval function $((Symbol("∇$(name)pool_direct!")))(
                    dx::AbstractArray{T,5}, dy::AbstractArray{T,5},
                    y::AbstractArray{T,5}, x::AbstractArray{T,5},
                    pdims::PoolDims, ::Val{K}; # == kernel_size(pdims)
                    alpha::T=T(1), beta::T=T(0)) where {T, K}
        check_dims(size(x), size(dy), pdims)

        width, height, depth = input_size(pdims)
        kernel_w, kernel_h, kernel_d = K
        out_c = channels_out(pdims)
        pad_w_lo, _, pad_h_lo, _, pad_d_lo, _ = padding(pdims)
        dil_w, dil_h, dil_d = dilation(pdims)
        stride_w, stride_h, stride_d = stride(pdims)

        # We use calc_padding_regions to split outselves up into separate regions that
        # may or may not need to worry about padding:
        padded_regions, central_region = calc_padding_regions(pdims)

        # A helper function to project from output (w, h) to input (input_w, input_h)
        @inline project(idx, stride, pad) = (idx - 1) * stride - pad + 1

        # If we're doing mean pooling, we represent division by kernel size by rolling
        # it into the `alpha` multiplier.
        if $(name == :mean)
            alpha = alpha / prod(K)
        end

        # Start with the central region
        w_region, h_region, d_region = central_region
        @inbounds for batch_idx in 1:size(x, 5), c in 1:out_c
            for d in d_region
            pd = project(d, stride_d, pad_d_lo)
            for h in h_region
            ph = project(h, stride_h, pad_h_lo)
            for w in w_region
            pw = project(w, stride_w, pad_w_lo)

            # Grab the output at this index for future use
            y_idx = y[w, h, d, c, batch_idx]
            dy_idx = dy[w, h, d, c, batch_idx]
            maxpool_already_chose = false

            for kd in 1:kernel_d,
                kh in 1:kernel_h,
                kw in 1:kernel_w

                input_kd = pd + (kd - 1) * dil_d
                input_kh = ph + (kh - 1) * dil_h
                input_kw = pw + (kw - 1) * dil_w

                # This conditional will be optimized away at compile time,
                # or my name isn't shengdan jingyu
                # x_idxs = (input_kw, input_kh, input_kd, c, batch_idx)
                if $(name == :max)
                    if maxpool_already_chose
                        break
                    end
                    # If it's equal; this is the one we chose. We only choose one per
                    # kernel window, all other elements of dx must be zero.
                    # Uncomment line below if using with non-precise output (e.g. by NNPACK)
                    # if abs(y_idx - x[x_idxs...]) < 1e-5 && !maxpool_already_chose
                    if y_idx ≈ x[input_kw, input_kh, input_kd, c, batch_idx]
                        dx[input_kw, input_kh, input_kd, c, batch_idx] += dy_idx * alpha #+ beta * dx[x_idxs...]
                        maxpool_already_chose = true
                    # Maxpooling does not support `beta` right now.  :(
                    # else
                    #    dx[x_idxs...] = T(0) + beta*dx[x_idxs...]
                    end
                elseif $(name == :mean)
                    # Either does meanpool :(
                    dx[input_kw, input_kh, input_kd, c, batch_idx] += dy_idx * alpha
                else
                    error("Unimplemented codegen path")
                end
            end
            end
            end
            end
        end

        # Next, the padded regions
        @inbounds for (w_region, h_region, d_region) in padded_regions
            for batch_idx in 1:size(x, 5), c in 1:out_c
                for d in d_region
                pd = project(d, stride_d, pad_d_lo)
                for h in h_region
                ph = project(h, stride_h, pad_h_lo)
                for w in w_region
                pw = project(w, stride_w, pad_w_lo)

                # Grab the incoming gradient at this index for future use
                y_idx = y[w, h, d, c, batch_idx]
                dy_idx = dy[w, h, d, c, batch_idx]
                maxpool_already_chose = false

                # In these loops, we have to check that we're not reaching off the edge,
                # we do so by putting in a bunch of conditionals.  :/
                for kd in 1:kernel_d
                    input_kd = pd + (kd - 1) * dil_d
                    if input_kd <= 0 || input_kd > depth
                        continue
                    end

                    for kh in 1:kernel_h
                        input_kh = ph + (kh - 1) * dil_h
                        if input_kh <= 0 || input_kh > height
                            continue
                        end

                        for kw in 1:kernel_w
                            input_kw = pw + (kw - 1) * dil_w
                            if input_kw <= 0 || input_kw > width
                                continue
                            end

                            # Same as above
                            # x_idxs = (input_kw, input_kh, input_kd, c, batch_idx)
                            if $(name == :max)
                                if maxpool_already_chose
                                    break
                                end
                                # Uncomment line below if using with non-precise output
                                # if abs(y_idx - x[x_idxs...]) < 1e-5 && !maxpool_already_chose
                                if y_idx ≈ x[input_kw, input_kh, input_kd, c, batch_idx]
                                    dx[input_kw, input_kh, input_kd, c, batch_idx] += dy_idx * alpha #+ beta * dx[x_idxs...]
                                    maxpool_already_chose = true
                                # else
                                #    dx[x_idxs...] = T(0) + beta*dx[x_idxs...]
                                end
                            elseif $(name == :mean)
                                dx[input_kw, input_kh, input_kd, c, batch_idx] += dy_idx * alpha #+ beta * dx[x_idxs...]
                            else
                                error("Unimplemented codegen path")
                            end
                        end
                    end
                end
            end
            end
            end
            end
        end

        return dx
    end
end


text:
"""
Simulations to estimate coverage of confidence bands with Blasques and Buccheri's methods
To Do :
- profile R interface to see if simulations for PMLE can be faster
-
"""

#region import and models

using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using JLD2
using Distributed
using SharedArrays
using ScoreDrivenERGM

begin
nWorkers = 12
addprocs(nWorkers - nprocs())
@sync @everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
    using ScoreDrivenERGM
    import ScoreDrivenERGM:StaticNets, DynNets
    import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, sample_ergm_sequence, stats_from_mat, array_2_vec_all_par, unrestrict_all_par, number_ergm_par, estimate_filter_and_conf_bands, conf_bands_coverage_parallel, estimate, plot_filtered_and_conf_bands
    using ScoreDrivenERGM.Utilities

    model_mle = DynNets.SdErgmDirBin0Rec0_mle()
    model_pmle = DynNets.SdErgmDirBin0Rec0_pmle()
    indTvPar = trues(2)

    end

end
#endregion

# #region coverage simulations
dgpSetAR, ~, dgpSetSD = ScoreDrivenERGM.DynNets.list_example_dgp_settings(model_mle)

c= Dict{String, Any}()
c["T"] = [100, 300, 600]
c["N"] = [100, 200, 300]
c["dgpSettings"] = [dgpSetSD, dgpSetAR]
c["model"] =[ model_pmle]
c["nSampleCoverage"] = 60

list = dict_list(c)

for d in dict_list(c)

    allCoverFiltPar, allCoverPar, allvEstSdResPar, allfVecT_filt, allParDgpT, allConfBandsFiltPar,allConfBandsPar, allErrFlags =  ScoreDrivenERGM.DynNets.conf_bands_coverage_parallel(d["model"], d["dgpSettings"], d["T"], d["N"],  d["nSampleCoverage"]; quantilesVals = [[0.975, 0.025]])

    res =  (;allCoverFiltPar, allCoverPar, allvEstSdResPar, allfVecT_filt, allParDgpT, allConfBandsFiltPar,allConfBandsPar, allErrFlags ) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

    #res = load(datadir("sims", "cofBandsCover", saveName))

    resDict = merge(res, d)

    saveName = replace.( savename(d, "jld";allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "")

   @tagsave( datadir("sims", "cofBandsCover", saveName), resDict)

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
"""
Overview: collections of functions required to solve the ε¹-invariance equation
"""

"""
> solve_homological_veps!(η,Cp,Cp⁺,ndofs,p,sys_mat,sys_rhs,sys_res,style,ic,M,C,K)
It solves al homological equations for a given order of the asymtptotic developmet of the ε¹-invariance.
- η : eigenvalue associated to the non-autonomous forcing
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- sys_mat : homological equation matrix
- sys_rhs : homological equation right hand side
- sys_res : homological equation residual
- style : parametrisation style
- ic : index combination
- M : mass matrix
- C : damping matrix
- K : stiffness matrix
"""
function solve_homological_veps!(η,Cp,Cp⁺,ndofs,p,sys_mat,sys_rhs,
                                 sys_res,style,ic,M,C,K)
  #
  resonant_modes = zeros(Bool,ndofs)
  nm = Int(ndofs/2)
  #
  for i = 1:Cp⁺[p+1].nc
    #
    σ⁺ = η
    #
    for j = 1:p
      ic[j] = Cp[p+1].comb[j,i]
      σ⁺ += Cp[2].f[ic[j],ic[j]]
    end
    #
    check_resonances!(Cp,σ⁺,ndofs,style,resonant_modes)
    assembly_sys_mat!(Cp,sys_mat.data,M.data,C.data,K.data,σ⁺,resonant_modes,ndofs)
    assembly_sys_rhs_veps!(Cp,Cp⁺,sys_rhs,sys_res,M.data,C.data,i,p,ndofs,σ⁺,resonant_modes)
    #
    sys_res = sparse(sys_mat)\sys_rhs
    #
    @inbounds for j = 1:M.data.m
      Cp⁺[p+1].W[j,i] = Cp⁺[p+1].W[j+M.data.m,i] + σ⁺*sys_res[j]
      Cp⁺[p+1].W[j+M.data.m,i] = sys_res[j]
    end
    for j = 1:ndofs
      if (resonant_modes[j]>0)
        Cp⁺[p+1].f[j,i] = sys_res[M.data.m+j]
      end
    end
    #
    for j = 1:ndofs
      @inbounds for jj = 1:M.data.m
        Cp⁺[p+1].W[jj,i] += Cp⁺[p+1].f[j,i]*Cp[2].W[jj+M.data.m,j]
      end
    end
    #
  end
  #
  return nothing
  #
end

"""
> assembly_sys_rhs!(Cp,Cp⁺,sys_rhs,sys_res,M,C,comb,p,ndofs,σ,resonant_modes)
It assembles the right hand side of a given homological equation of the ε⁰-developmet
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- sys_rhs : right hand side of the homological equations
- sys_res : residual of the homological equations
- M : mass matrix
- C : damping matrix
- comb : reference monomial index
- p : order of the asymptotic development
- ndofs : number of degrees of freedom
- σ : eingevalues λ summation
- resonant_modes : book-keeping of the set R of resonant modes
"""
function assembly_sys_rhs_veps!(Cp,Cp⁺,sys_rhs,sys_res,
                                M,C,comb,p,ndofs,
                                σ,resonant_modes)
  #
  fill!(sys_rhs,0.0)
  fill!(sys_res,0.0)
  nm = Int(ndofs/2)
  #
  @inbounds for irow = 1:M.m
    pos1 = M.colptr[irow]
    pos2 = M.colptr[irow+1]-1

    for j = pos1:pos2
      jcol = M.rowval[j]
      #
      sys_rhs[irow] -= M.nzval[j]*Cp⁺[p+1].W[jcol,comb] +
                       (σ*M.nzval[j]+C.nzval[j])*Cp⁺[p+1].W[jcol+M.m,comb]
      sys_res[irow] -= M.nzval[j]*Cp⁺[p+1].W[jcol+M.m,comb]
      #
      if (irow!=jcol)
        #
        sys_rhs[jcol] -= M.nzval[j]*Cp⁺[p+1].W[irow,comb] +
                         (σ*M.nzval[j]+C.nzval[j])*Cp⁺[p+1].W[irow+M.m,comb]
        sys_res[jcol] -= M.nzval[j]*Cp⁺[p+1].W[irow+M.m,comb]
        #
      end
    end
  end
  # resonant modes
  for i = 1:ndofs
    if (resonant_modes[i]>0)
      if (i<=nm)
        @inbounds for j = 1:M.m
          sys_rhs[M.m+i] += Cp[2].W[M.m+j,i]*sys_res[j]
        end
      else
        @inbounds for j = 1:M.m
          sys_rhs[M.m+i] += Cp[2].W[M.m+j,i-nm]*sys_res[j]
        end
      end
    end
  end
  # elastic nl
  @inbounds for i = 1:M.m
    sys_rhs[i] -= Cp⁺[p+1].nlr[i,comb]
  end
  #
  return nothing
  #
end

"""
> recursive_assembly_veps!(Cp,Cp⁺,ndofs,p,mesh,U,Cic,cc,nls=0)
It assembles right hand sides of all homological equations for a given order p of the expansion.
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- mesh : grid data structure
- U : displacement field
- Cic : monomials index combination
- cc : recursion depth
- nls : nonlinear static solution. 0 = no, 1 = yes.
"""
function recursive_assembly_veps!(Cp,Cp⁺,ndofs,p,mesh,U,Cic,cc,nls=0)
  #
  for i = 1:ndofs
    Cic[cc] = i
    if (cc<p)
      cc += 1
      recursive_assembly_veps!(Cp,Cp⁺,ndofs,p,mesh,U,Cic,cc,nls)
      cc -= 1
    else
      #
      pos = find_index_global(Cic,ndofs,p)
      entry = Cp⁺[p+1].cmap[pos]
      #
      assembly_quadratic_nl_veps!(Cp,Cp⁺,entry,Cic,ndofs,p,mesh,U,nls)
      assembly_cubic_nl_veps!(Cp,Cp⁺,entry,Cic,ndofs,p,mesh,U)
      assembly_μ_ν_resid_veps!(Cp,Cp⁺,entry,Cic,ndofs,p,U.neq)
      #
    end
  end
  #
  return nothing
  #
end

"""
> recursive_assembly_veps!(Cp,Cp⁺,ndofs,p,mesh,U,Cic,cc,nls=0)
It assembles right hand sides of all homological equations for a given order p of the expansion.
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- entry : array entry of a given monomial
- ic : index combination associated to a monomial
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- neq : number of equations
"""
function assembly_μ_ν_resid_veps!(Cp,Cp⁺,entry,ic,ndofs,p,neq)
  #
  ent = abs(entry)
  #
  for k = 2:p
    for l = 0:p-k
      posf = find_index_rdyn(k,l,ndofs,ic)
      cmf = Cp[k+1].cmap[posf]
      #
      if (cmf>0)
        f = @view Cp[k+1].f[:,cmf]
        #
        for s = 1:ndofs
          posw = find_index_map(k,l,s,p,ndofs,ic)
          cmw = Cp⁺[p-k+2].cmap[posw]
          if (cmw>0)
            W⁺ = Cp⁺[p-k+2].W[:,cmw]
            @inbounds for i = 1:neq*2
              Cp⁺[p+1].W[i,ent] += f[s]*W⁺[i]
            end
          end
        end
      end
    end
  end
  #
  for k = 0:p-1
    for l = 0:p-k
      #
      posf = find_index_rdyn(k,l,ndofs,ic)
      #
      cmf = Cp⁺[k+1].cmap[posf]
      #
      if (cmf>0)
        f⁺ = @view Cp⁺[k+1].f[:,cmf]
        for s = 1:ndofs
          #
          posw = find_index_map(k,l,s,p,ndofs,ic)
          cmw = Cp[p-k+2].cmap[posw]
          #
          if (cmw>0)
            W = @view Cp[p-k+2].W[:,cmw]
            @inbounds for i = 1:neq*2
              Cp⁺[p+1].W[i,ent] += f⁺[s]*W[i]
            end
          end
        end
      end
    end
  end
  #
  return nothing
  #
end

"""
> assembly_cubic_nl_veps!(Cp,Cp⁺,entry,Cic,ndofs,p,mesh,U)
It computes cubic nonlinearity vectors for all monomials of order p of the non-autonomous asymptotic development
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- entry : array entry of a given monomial
- ic : index combination associated to a monomial
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- mesh : grid data structure
- U : displacement field
"""
function assembly_cubic_nl_veps!(Cp,Cp⁺,entry,ic,ndofs,p,mesh,U)
  #
  for k = 1:p-1
    for l = 1:p-k
      #
      pos1,pos2,pos3 = find_position_h(k,l,p,ndofs,ic)
      #
      cm1 = Cp[k+1].cmap[pos1]
      cm2 = Cp[l+1].cmap[pos2]
      cm3 = Cp⁺[p-k-l+1].cmap[pos3]
      #
      if (cm1>0 && cm2>0 && cm3>0)
        #
        Ψ₁ = @view Cp[k+1].W[U.neq+1:2*U.neq,cm1]
        Ψ₂ = @view Cp[l+1].W[U.neq+1:2*U.neq,cm2]
        Ψ₃ = @view Cp⁺[p-k-l+1].W[U.neq+1:2*U.neq,cm3]
        #
        assembly_H_nl!(Cp⁺[p+1], abs(entry), Ψ₁, Ψ₂, Ψ₃, mesh, U, 3.0)
        #
      end
      #
    end
  end
  #
  return nothing
  #
end

"""
> assembly_cubic_nl_veps!(Cp,Cp⁺,entry,Cic,ndofs,p,mesh,U)
It computes quadratic nonlinearity vectors for all monomials of order p of the non-autonomous asymptotic development
- Cp : autonomous parametrisation data structure
- Cp⁺ :  non-auonomous parametrisation data structure
- entry : array entry of a given monomial
- ic : index combination associated to a monomial
- ndofs : number of degrees of freedom
- p : order of the asymptotic development
- mesh : grid data structure
- U : displacement field
- nls : nonlinear static solution. 0 = no, 1 = yes
"""
function assembly_quadratic_nl_veps!(Cp,Cp⁺,entry,ic,ndofs,p,mesh,U,nls)
  #
  for k = 1:p
    pos1, pos2 = find_position_g(k,p,ndofs,ic)
    #
    cm1 = Cp[k+1].cmap[pos1]
    cm2 = Cp⁺[p-k+1].cmap[pos2]
    #
    #
    if (cm1>0 && cm2>0)
      #
      Ψ₁ = @view Cp[k+1].W[U.neq+1:2*U.neq,cm1]
      Ψ₂ = @view Cp⁺[p-k+1].W[U.neq+1:2*U.neq,cm2]
      #
      assembly_G_nl!(Cp⁺[p+1], abs(entry), Ψ₁, Ψ₂, mesh, U, 2.0)
      #
      if (nls>0)
        assembly_H_nl!(Cp⁺[p+1], abs(entry), Ψ₁, Ψ₂, Cp[1].W, mesh, U, 6.0)
      end
      #
    end
    #
  end
  #
  return nothing
  #
end


text:
module WeakValueDicts

export WeakValueDict

using Base: Callable, ValueIterator, secret_table_token

"""
    WeakValueDict([itr])

`WeakValueDict()` constructs a hash table where the values are weak references
to objects which may be garbage collected even when referenced in a hash table.

See [`Dict`](@ref) for further help.
"""
mutable struct WeakValueDict{K,V} <: AbstractDict{K,V}
    ht::Dict{K,WeakRef}
    lock::ReentrantLock
    finalizer

    # Constructors mirror Dict's
    function WeakValueDict{K,V}() where {V} where {K}
        wvd = new(Dict{Any,V}(), ReentrantLock())
        wvd.finalizer = function (k, v)
            # When a weak value is finalized, remove from dictionary if it is
            # still there.
            if islocked(wvd)
                # If locked, we add another finalizer and defer the deletion
                # until the new finalizer is invoked.
                finalizer((v) -> wvd.finalizer(k, v), v)
                return nothing
            end
            delete!(wvd, k)
        end
        return wvd
    end
end

WeakValueDict() = WeakValueDict{Any,Any}()
WeakValueDict(kv::Tuple{}) = WeakValueDict()
WeakValueDict(ps::Pair{K,V}...) where {K,V} = WeakValueDict{K,V}(ps)
WeakValueDict(ps::Pair{K}...) where {K} = WeakValueDict{K,Any}(ps)
WeakValueDict(ps::(Pair{K,V} where {K})...) where {V} = WeakValueDict{Any,V}(ps)
WeakValueDict(ps::Pair...) = WeakValueDict{Any,Any}(ps)

function WeakValueDict{K,V}(kv) where {V} where {K}
    h = WeakValueDict{K,V}()
    for (k, v) in kv
        h[k] = v
    end
    return h
end

function WeakValueDict{K,V}(p::Pair) where {V} where {K}
    setindex!(WeakValueDict{K,V}(), p.second, p.first)
end

function WeakValueDict{K,V}(ps::Pair...) where {V} where {K}
    h = WeakValueDict{K,V}()
    sizehint!(h, length(ps))
    for p in ps
        h[p.first] = p.second
    end
    return h
end

function WeakValueDict(kv)
    try
        Base.dict_with_eltype((K, V) -> WeakValueDict{K,V}, kv, eltype(kv))
    catch
        if !Base.isiterable(typeof(kv)) || !all(x -> isa(x, Union{Tuple,Pair}), kv)
            throw(ArgumentError(
                "WeakValueDict(kv): " *
                    "kv needs to be an iterator of tuples or pairs",
            ))
        else
            rethrow()
        end
    end
end

Base.copy(d::WeakValueDict) = WeakValueDict(d)

function Base.empty(d::WeakValueDict, ::Type{K}, ::Type{V}) where {K,V}
    return WeakValueDict{K,V}()
end

Base.islocked(wvd::WeakValueDict) = islocked(wvd.lock)
Base.lock(f, wvd::WeakValueDict) = lock(f, wvd.lock)

function Base.getindex(wvd::WeakValueDict{K, V}, key::K)::V where {K, V}
    return lock(wvd) do
        return getindex(wvd.ht, key).value::V
    end
end

function Base.setindex!(wvd::WeakValueDict{K, V}, value::V, key::K) where {K,V}
    finalizer((value) -> wvd.finalizer(key, value), value)
    lock(wvd) do
        wvd.ht[key] = WeakRef(value)
    end
    return wvd
end

function Base.getkey(wvd::WeakValueDict{K}, kk, default) where {K}
    return lock(wvd) do
        k = getkey(wvd.ht, kk, secret_table_token)
        k === secret_table_token && return default
        return k
    end
end

function Base.map!(f, iter::ValueIterator{<:WeakValueDict})
    dict = iter.dict.ht

    # This code is adapted from the implementation of `map!` for WeakKeyDict.
    # Here be `Dict` internals and violations of encapsulation.
    # Beware all ye who enter here.
    vals = dict.vals
    # @inbounds is here so the it gets propigated to isslotfiled
    @inbounds for i = dict.idxfloor:lastindex(vals)
        if Base.isslotfilled(dict, i)
            # Hold an explicit reference here so that we avoid race conditions
            # where the GC can trigger between checking if the value is
            # `nothing` and actually mapping it.
            value = vals[i].value
            if value === nothing
                continue
            end
            vals[i].value = f(value)
        end
    end
    return iter
end

function Base.get(wvd::WeakValueDict{K}, key, default) where {K}
    return get(() -> default, wvd, key)
end

function Base.get(default::Callable, wvd::WeakValueDict{K, V}, key) where {K, V}
    return lock(wvd) do
        if haskey(wvd.ht, key)
            return wvd.ht[key].value::V
        end
        return default()
    end
end

function Base.get!(wvd::WeakValueDict{K, V}, key::K, default::V) where {K,V}
    return get!(() -> default, wvd, key)
end

function Base.get!(default::Callable, wvd::WeakValueDict{K, V}, key::K) where {K,V}
    return lock(wvd) do
        if haskey(wvd.ht, key)
            value = wvd.ht[key].value
            if value !== nothing
                return value
            end
        end
        return (wvd[key] = default()::V)
    end
end

function Base.pop!(wvd::WeakValueDict{K, V}, key)::V where {K,V}
    return lock(wvd) do
        ref = pop!(wvd.ht, key)
        # TODO: Can this ever be nothing? Can a WeakRef's value be nothing
        # before the finalizers run?
        return ref.value::V
    end
end

function Base.pop!(wvd::WeakValueDict{K, V}, key, default) where {K, V}
    return lock(wvd) do
        if haskey(wvd.ht)
            ref = pop!(wvd.ht, key)
            return ref.value::V
        end
        return default
    end
end

function Base.delete!(wvd::WeakValueDict, key)
    lock(wvd) do
        delete!(wvd.ht, key)
        return wvd
    end
end
Base.empty!(wvd::WeakValueDict) = (lock(() -> empty!(wvd.ht), wvd); wvd)
Base.haskey(wvd::WeakValueDict{K}, key) where {K} =
    lock(() -> haskey(wvd.ht, key), wvd)

Base.isempty(wvd::WeakValueDict) = isempty(wvd.ht)
Base.length(t::WeakValueDict) = length(t.ht)
Base.sizehint!(wvd::WeakValueDict, newsz) = sizehint!(wvd.ht, newsz)

function Base.iterate(t::WeakValueDict{K,V}, state...) where {K,V}
    y = lock(() -> iterate(t.ht, state...), t)
    y === nothing && return nothing
    wkv, newstate = y
    key, value_ref = wkv
    kv = Pair{K,V}(key, value_ref.value::V)
    return (kv, newstate)
end

Base.filter!(f, d::WeakValueDict) = Base.filter_in_one_pass!(f, d)

end # module