text:
struct PhysicalConfiguration
    Dom::Domain
    S::AbstractMaterial
end

struct PhysicalBoundaryConditions{T}
    position
    values::AbstractVector
    variable
end

struct SemiconductorProblem
    physicalConfiguration::Vector{PhysicalConfiguration}
    #bcs::AbstractVector{PhysicalBoundaryConditions{T}}
end

function DoppingFunction(dom::Domain,C::Number)
    return f(x) = (x ∈ dom)*C
end

# function RelaxStep(S::Semiconductor1D,x::AbstractVector)
#     applyPsiProfile(S,x)
#     return relaxationCostFunction(S)
# end

# function Relax(S::Semiconductor1D)
#     error(x) = RelaxStep(S,x)
#     return optimize(error,0.0*collect(1:length(S)),GradientDescent())
# end

function GummelPoisson(PP::Semiconductor1D,ψ₀,ψₗ,α=10.0;abserror=1e-12,itermax=100)
    n_i = ni(PP[1])
    ϕₜ = kb*Temperature(PP)
    Cdoping = C(PP)
    Cscal = maximum(vcat([n_i],abs.(Cdoping)))
    lscal = Rₗ(PP)
    Nnodes = length(PP)
    ui = -Ei(PP)/ϕₜ
    return GummelPoisson(ψ₀/ϕₜ + ui[1],ψₗ/ϕₜ + ui[end],ui,Cdoping,Nnodes,n_i,lscal,ϕₜ,r(PP),α;abserror=abserror,itermax=itermax)*ϕₜ
    return r0*lscal,u*ϕₜ,e
end

function GummelPoisson(u₀,uₗ,ui,C,Nnodes,n_i,lscal,Vscal,r0,α=1.0;abserror=1e-12,itermax=100,relerror=1e-12)
    Cscal = maximum(vcat([n_i],abs.(C)))
    λ = ϵ₀*11.7*Vscal/(lscal^2*q*Cscal)
    r0 ./=lscal
    C = C./Cscal
    δ = n_i/Cscal
    dr = diff(r0)
    #Operators
    ∇² = CenteredDifference(2,5,dr,Nnodes-2)
    Bc =  DirichletBC(u₀,uₗ)
    φBC = DirichletBC(0.0,0.0)
    ## Initialization
    u0 = ui
    u0[1] = u₀
    u0[end] = uₗ
    ## Step
    u1 = copy(u0)
    limexp(x) = x > 1e2 ? exp(1e2) : exp(x)
    n(u) = δ*limexp(u)
    p(u) = δ*limexp(-u)
    R(u) = -λ*(∇²*Bc*u)+δ*limexp.(u)-δ*limexp.(-u) - C[2:end-1]
    e = Inf64
    iter = 1
    Rmin = Inf64
    uf = similar(u0)
    nR = Inf64
    while e > abserror && iter < itermax && nR > relerror
        res = R(u0[2:end-1])
        SGL(u) = -λ*(∇²*φBC*u) + (n.(u)+p.(u)).*u + res
        ϕ = zeros(length(u0))
        ϕ[2:end-1] .= nlsolve(SGL,ϕ[2:end-1],autodiff = :forward).zero
        b = abs.(ϕ) .> 1.0
        # This damping procedure is inspired from Solid-State Electronics, vol. 19,
        # pp. 991-992 (1976).
        ϕ[b] = log.(1.0.+abs.(ϕ[b]).*1.72).*sign.(ϕ[b])
        u1[2:end-1] .+= ϕ[2:end-1]
        nR = sum(abs.(ϕ[2:end-1]))/length(res)
        e = norm(res)
        #println("i:: $(iter) norm:$(nR) res:$(e)")
        if(nR < Rmin)
            uf = u0
            Rmin = nR
        end
        u0 = u1
        iter = iter+1

    end
    return uf
end

function scallingNv(S::SemiconductorData,L)
    n_i = ni(S)
    ϕₜ = kb*Temperature(S)
    Vscal = ϕₜ
    Cscal = max(n_i,abs.(C(S)),Nv(S))
    lscal = L
    λ = ϵ₀*11.7*Vscal/(lscal^2*q*Cscal)
    δ = n_i/Cscal
    μscal = 1000
    Ld0 = sqrt(ϵ₀*11.7* ϕₜ/q/abs(C(S)))
    Rscal = ϕₜ*μscal*Cscal/lscal/lscal
    tscal = lscal^2/ϕₜ/μscal
    Jscal = q*Vscal*μscal*Cscal/lscal
    Dscal = μscal*Vscal
    return Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i

end

function scallingLD(S::SemiconductorData,L)
    n_i = ni(S)
    ϕₜ = kb*Temperature(S)
    Vscal = ϕₜ
    Cscal = max(n_i,abs.(C(S)))
    lscal = sqrt(ϵ₀*11.7* ϕₜ/q/abs(C(S)))
    λ = ϵ₀*11.7*Vscal/(lscal^2*q*Cscal)
    δ = n_i/Cscal
    μscal = 1000
    Ld0 = sqrt(ϵ₀*11.7* ϕₜ/q/abs(C(S)))
    Rscal = ϕₜ*μscal*Cscal/lscal/lscal
    tscal = lscal^2/ϕₜ/μscal
    Jscal = q*Vscal*μscal*Cscal/lscal
    Dscal = μscal*Vscal
    return Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i

end

function scallingNvLD(S::SemiconductorData,L)
    n_i = ni(S)
    ϕₜ = kb*Temperature(S)
    Vscal = ϕₜ
    Cscal = max(n_i,abs.(C(S)),Nv(S))
    lscal = sqrt(ϵ₀*11.7* ϕₜ/q/abs(C(S)))
    λ = ϵ₀*11.7*Vscal/(lscal^2*q*Cscal)
    δ = n_i/Cscal
    μscal = 1000
    Ld0 = sqrt(ϵ₀*11.7* ϕₜ/q/abs(C(S)))
    Rscal = ϕₜ*μscal*Cscal/lscal/lscal
    tscal = lscal^2/ϕₜ/μscal
    Jscal = q*Vscal*μscal*Cscal/lscal
    Dscal = μscal*Vscal
    return Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i

end

function limexp(x,lm=-50,lM=50)
    if x < lm
        return exp(lm)*((x-lm)+1.0)
    elseif x > lM
        return exp(lM)*((x-lM)+1.0)
    else
        return exp(x)
    end

end
function DDVoronoiMOSFET(X,Ec,Ev,Nc,Nv,Vbi,Vg,Vd,lscal,Rscal,Cscal,Vscal,λ,Jn,Jp,R,C,fetIn::MOSFETInputDeck;u0=nothing,n0=nothing,p0=nothing,ξ₀=1e-5,verbose=false)
    ## Selective Space Discretization

    #Lint = maximum(X)
    ## Boltzman Dist
    F(x) = limexp(x)
    ηₙ(u,un) = u-Ec/Vscal
    ncf(u,un) = Nc*F(ηₙ(u,un))/Cscal
    ηₚ(u,up) = -u+Ev/Vscal
    pcf(u,up) = Nv*F(ηₚ(u,up))/Cscal
    nc(u,un) = un
    pc(u,up) = up

    #Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])/meas(edge)

                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))/meas(edge)
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))/meas(edge)
                # @info "u $(map(x->x.value,u.val)) - m: $(meas(edge)) - f:$(map(x->x.value,f.val))"
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal,node[2]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(X,physics)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    Vs = 0.0
    us = 0.0
    ug = (Vbi+Vg)/Vscal
    ud = (Vd)/Vscal
    uns = ncf(us,0.0)*exp(-Vs/Vscal)
    ung = ncf(ug,0.0)*exp(-Vg/Vscal)
    und = ncf(ud,0.0)exp(-Vd/Vscal)
    ups = pcf(us,0.0)*exp(Vs/Vscal)
    upg = pcf(ug,0.0)exp(Vg/Vscal)
    upd = pcf(ud,0.0)exp(Vd/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,1,us)
    boundary_dirichlet!(sys,1,6,us)
    boundary_dirichlet!(sys,1,7,ug)
    boundary_dirichlet!(sys,1,8,ud)

    boundary_dirichlet!(sys,2,1,uns)
    boundary_dirichlet!(sys,2,6,uns)
    boundary_dirichlet!(sys,2,7,ung)
    boundary_dirichlet!(sys,2,8,und)

    boundary_dirichlet!(sys,3,1,ups)
    boundary_dirichlet!(sys,3,6,ups)
    boundary_dirichlet!(sys,3,7,upg)
    boundary_dirichlet!(sys,3,8,upd)

    uv = [us,ug,ud]
    nv = [uns,ung,und]
    pv = [ups,upg,upd]
    Cords = X[Coordinates]
    # Create a solution array
    inival=unknowns(sys,inival=0.1)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        for i in 1:size(inival)[2]
            inival[1,i] = dot(belongVector(Cords[1,i]*lscal,Cords[2,i]*lscal,fetIn),uv)
        end
    else
        inival[1,:] .= u0
    end
    if isnothing(n0)
        for i in 1:size(inival)[2]
            inival[2,i] = dot(belongVector(Cords[1,i]*lscal,Cords[2,i]*lscal,fetIn),nv)
        end
    else
        inival[2,:] = n0
    end
    if isnothing(p0)
        for i in 1:size(inival)[2]
            inival[3,i] = dot(belongVector(Cords[1,i]*lscal,Cords[2,i]*lscal,fetIn),pv)
        end
    else
        inival[3,:] .=p0
    end

    solution=unknowns(sys)

    # Create solver control info
    control=VoronoiFVM.NewtonControl()
    #control.max_lureuse = 0
    control.verbose=verbose
    control.max_iterations = 500
    #control.Δu_opt=100
    control.damp_initial= ξ₀
    # Stationary solution of the problem
    VoronoiFVM.solve!(solution,inival,sys, control=control)

    return solution[1,:]*Vscal,solution[2,:]*Cscal,solution[3,:]*Cscal
end

function DDVoronoiMOSFET(X,Ec,Ev,Nc,Nv,Vbi,Vg,Vd,lscal,Rscal,Cscal,Vscal,Ld0,L,λ,nld,Jn,Jp,nn,R,C;u0=nothing,n0=nothing,p0=nothing,ξ₀=1e-5,verbose=false)
    ## Selective Space Discretization

    #Lint = maximum(X)
    ## Boltzman Dist
    F(x) = abs(x) > 100 ? exp(sign(x)*100) : exp(x)
    ηₙ(u,un) = u-Ec/Vscal
    ncf(u,un) = Nc*F(ηₙ(u,un))/Cscal
    ηₚ(u,up) = -u+Ev/Vscal
    pcf(u,up) = Nv*F(ηₚ(u,up))/Cscal
    nc(u,un) = un
    pc(u,up) = up

    #Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])
                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal,node[2]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(X,physics,unknown_storage=:sparse)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    Vs = 0.0
    us = 0.0
    ug = (Vbi+Vg)/Vscal
    ud = (Vd)/Vscal
    uns = ncf(us,0.0)*exp(-Vs/Vscal)
    ung = ncf(ug,0.0)*exp(-Vg/Vscal)
    und = ncf(ud,0.0)exp(-Vd/Vscal)
    ups = pcf(us,0.0)*exp(Vs/Vscal)
    upg = pcf(ug,0.0)exp(Vg/Vscal)
    upd = pcf(ud,0.0)exp(Vd/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,1,us)
    boundary_dirichlet!(sys,1,6,us)
    boundary_dirichlet!(sys,1,7,ug)
    boundary_dirichlet!(sys,1,8,ud)

    boundary_dirichlet!(sys,2,1,uns)
    boundary_dirichlet!(sys,2,6,uns)
    boundary_dirichlet!(sys,2,7,ung)
    boundary_dirichlet!(sys,2,8,und)

    boundary_dirichlet!(sys,3,1,ups)
    boundary_dirichlet!(sys,3,6,ups)
    boundary_dirichlet!(sys,3,7,upg)
    boundary_dirichlet!(sys,3,8,upd)

    # Create a solution array
    inival=unknowns(sys,inival=0.1)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        inival[1,1:L] .= us
        inival[1,L+1:2L] .= ud
    else
        inival[1,:] .= u0
    end
    if isnothing(n0)
        inival[2,1:L] .= uns
        inival[2,L+1:2L] .= und

    else
        inival[2,:] .= n0
    end
    if isnothing(p0)
        inival[3,1:L] .= ups
        inival[3,L+1:2L] .= upd
    else
        inival[3,:] .= p0
    end

    solution=unknowns(sys)

    # Create solver control info
    control=VoronoiFVM.NewtonControl()
    #control.max_lureuse = 0
    control.verbose=verbose
    control.max_iterations = 500
    #control.Δu_opt=100
    control.damp_initial= ξ₀
    # Stationary solution of the problem
    VoronoiFVM.solve!(solution,inival,sys, control=control)

    return solution[1,:]*Vscal,solution[2,:]*Cscal,solution[3,:]*Cscal
end

function DDVoronoiHC2D(X,Ec,Ev,Nc,Nv,Vbi,Va,Vb,lscal,Rscal,Cscal,Vscal,Ld0,L,λ,nld,Jn,Jp,nn,R,C,u0=nothing,n0=nothing,p0=nothing)
    ## Selective Space Discretization

    #Lint = maximum(X)
    ## Boltzman Dist
    F(x) = abs(x) > 100 ? exp(sign(x)*100) : exp(x)
    ηₙ(u,un) = u-Ec/Vscal
    ncf(u,un) = Nc*F(ηₙ(u,un))/Cscal
    ηₚ(u,up) = -u+Ev/Vscal
    pcf(u,up) = Nv*F(ηₚ(u,up))/Cscal
    nc(u,un) = un
    pc(u,up) = up

    #Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])/meas(edge)
                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))/meas(edge)
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))/meas(edge)
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(X,physics,unknown_storage=:sparse)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    ul = 0.0+Va/Vscal
    ur = (Vbi+Vb)/Vscal
    unr = ncf(ul,0.0)*exp(-Va/Vscal)
    unl = ncf(ur,0.0)exp(-Vb/Vscal)
    upr = pcf(ul,0.0)*exp(Va/Vscal)
    upl = pcf(ur,0.0)exp(Vb/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,4,ul)
    boundary_dirichlet!(sys,1,2,ur)
    boundary_dirichlet!(sys,2,4,unr)
    boundary_dirichlet!(sys,2,2,unl)
    boundary_dirichlet!(sys,3,4,upr)
    boundary_dirichlet!(sys,3,2,upl)

    # Create a solution array
    inival=unknowns(sys,inival=0.1)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        inival[1,1:L] .= ul
        inival[1,L+1:2L] .= ur
    else
        inival[1,:] .= u0
    end
    if isnothing(n0)
        inival[2,:] .= n0
    else
        inival[2,1:L] .= unr
        inival[2,L+1:2L] .= unl
    end
    if isnothing(p0)
        inival[3,1:L] .= upr
        inival[3,L+1:2L] .= upl
    else

    end

    solution=unknowns(sys)

    # Create solver control info
    control=VoronoiFVM.NewtonControl()
    #control.max_lureuse = 0
    control.verbose=true
    control.max_iterations = 500
    #control.Δu_opt=100
    control.damp_initial=1e-10
    # Stationary solution of the problem
    VoronoiFVM.solve!(solution,inival,sys, control=control)

    return solution[1,:]*Vscal,solution[2,:],solution[3,:]
end

function DDVoronoiHC(Ec,Ev,Nc,Nv,Vbi,Va,Vb,lscal,Rscal,Cscal,Vscal,Ld0,L,λ,nld,Jn,Jp,nn,R,C,u0=nothing)
    ## Selective Space Discretization
    dx0 = 3*Ld0/L
    dxld = dx0/nld
    dxS0 = dx0/nn
    p0 = 0
    p1 = dx0
    p2 = 0.5-dx0
    p3 = 0.5+dx0
    p4 = 1.0-dx0
    p5 = 1.0
    X1 = p0:dxld:p1
    X2 = p1:dxS0:p2
    X3 = p2:dxld:p3
    X4 = p3:dxS0:p4
    X5 = p4:dxld:p5
    X = vcat(X1,X2[2:end],X3[2:end],X4[2:end],X5[2:end])*L/lscal
    Lint = maximum(X)
    ## Boltzman Dist
    F(x) = abs(x) > 100 ? exp(sign(x)*100) : exp(x)
    ηₙ(u,un) = u-Ec/Vscal
    ncf(u,un) = Nc*F(ηₙ(u,un))/Cscal
    ηₚ(u,up) = -u+Ev/Vscal
    pcf(u,up) = Nv*F(ηₚ(u,up))/Cscal
    nc(u,un) = un
    pc(u,up) = up

    Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])/meas(edge)
                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))/meas(edge)
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))/meas(edge)
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(Gx,physics,unknown_storage=:sparse)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    ul = 0.0+Va/Vscal
    ur = (Vbi+Vb)/Vscal
    unr = ncf(ul,0.0)*exp(-Va/Vscal)
    unl = ncf(ur,0.0)exp(-Vb/Vscal)
    upr = pcf(ul,0.0)*exp(Va/Vscal)
    upl = pcf(ur,0.0)exp(Vb/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,1,ul)
    boundary_dirichlet!(sys,1,2,ur)
    boundary_dirichlet!(sys,2,1,unr)
    boundary_dirichlet!(sys,2,2,unl)
    boundary_dirichlet!(sys,3,1,upr)
    boundary_dirichlet!(sys,3,2,upl)

    # Create a solution array
    inival=unknowns(sys,inival=0.1)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        inival[1,1:L] .= ul
        inival[1,L+1:2L] .= ur
    else
        inival[1,:] .= u0
    end
    inival[2,1:L] .= unr
    inival[2,L+1:2L] .= unl
    inival[3,1:L] .= upr
    inival[3,L+1:2L] .= upl

    solution=unknowns(sys)

    # Create solver control info
    control=VoronoiFVM.NewtonControl()
    #control.max_lureuse = 0
    control.verbose=true
    control.max_iterations = 500
    #control.Δu_opt=100
    control.damp_initial=1e-10
    # Stationary solution of the problem
    VoronoiFVM.solve!(solution,inival,sys, control=control)

    return X*lscal,solution[1,:]*Vscal,solution[2,:],solution[3,:]
end

function DDVoronoiH(Ec,Ev,Nc,Nv,Vbi,Va,Vb,lscal,Ld0,L,λ,nld,Jn,Jp,nn,R,C,u0=nothing)
    ## Selective Space Discretization
    dx0 = 3*Ld0/L
    dxld = dx0/nld
    dxS0 = dx0/nn
    p0 = 0
    p1 = dx0
    p2 = 0.5-dx0
    p3 = 0.5+dx0
    p4 = 1.0-dx0
    p5 = 1.0
    X1 = p0:dxld:p1
    X2 = p1:dxS0:p2
    X3 = p2:dxld:p3
    X4 = p3:dxS0:p4
    X5 = p4:dxld:p5
    X = vcat(X1,X2[2:end],X3[2:end],X4[2:end],X5[2:end])*L/lscal
    Lint = maximum(X)
    ## Boltzman Dist
    F(x) = abs(x) > 100 ? exp(sign(x)*100) : exp(x)
    ηₙ(u,un) = u-Ec/Vscal
    nc(u,un) = Nc*F(ηₙ(u,un))*un/Cscal
    ηₚ(u,up) = -u+Ev/Vscal
    pc(u,up) = Nv*F(ηₚ(u,up))*up/Cscal

    Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])/meas(edge)
                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))/meas(edge)
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))/meas(edge)
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(Gx,physics,unknown_storage=:sparse)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    ul = 0.0+Va/Vscal
    ur = (Vbi+Vb)/Vscal
    unr = exp(-Va/Vscal)
    unl = exp(-Vb/Vscal)
    upr = exp(Va/Vscal)
    upl = exp(Vb/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,1,ul)
    boundary_dirichlet!(sys,1,2,ur)
    boundary_dirichlet!(sys,2,1,unr)
    boundary_dirichlet!(sys,2,2,unl)
    boundary_dirichlet!(sys,3,1,upr)
    boundary_dirichlet!(sys,3,2,upl)

    # Create a solution array
    inival=unknowns(sys,inival=0.1)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        inival[1,1:L] .= ul
        inival[1,L+1:2L] .= ur
    else
        inival[1,:] .= u0
    end
    inival[2,1:L] .= unr
    inival[2,L+1:2L] .= unl
    inival[3,1:L] .= upr
    inival[3,L+1:2L] .= upl

    solution=unknowns(sys)

    # Create solver control info
    control=VoronoiFVM.NewtonControl()
    #control.max_lureuse = 0
    control.verbose=true
    control.max_iterations = 500
    #control.Δu_opt=100
    control.damp_initial=1e-10
    # Stationary solution of the problem
    solve!(solution,inival,sys, control=control)

    return X*lscal,solution[1,:]*Vscal,-log.(solution[2,:])*Vscal,log.(solution[3,:])*Vscal
end

function DDVoronoi(Ec,Ev,Nc,Nv,Vbi,Va,Vb,lscal,Ld0,L,λ,nld,Jn,Jp,nn,R,C,control,u0=nothing)
    ## Selective Space Discretization
    dx0 = 3*Ld0/L
    dxld = dx0/nld
    dxS0 = dx0/nn
    p0 = 0
    p1 = dx0
    p2 = 0.5-dx0
    p3 = 0.5+dx0
    p4 = 1.0-dx0
    p5 = 1.0
    X1 = p0:dxld:p1
    X2 = p1:dxS0:p2
    X3 = p2:dxld:p3
    X4 = p3:dxS0:p4
    X5 = p4:dxld:p5
    X = vcat(X1,X2[2:end],X3[2:end],X4[2:end],X5[2:end])*L/lscal
    Lint = maximum(X)
    ## Boltzman Dist
    F(x) = abs(x) > 100 ? exp(sign(x)*100) : exp(x)
    ηₙ(u,un) = u-un-Ec/Vscal
    nc(u,un) = Nc*F(ηₙ(u,un))/Cscal
    ηₚ(u,up) = up-u+Ev/Vscal
    pc(u,up) = Nv*F(ηₚ(u,up))/Cscal

    Gx = VoronoiFVM.Grid(X)

    #Create PhysicalModel

    physics=VoronoiFVM.Physics(
        flux= function flux!(f,u,edge)
                f[1]=λ*(u[1,1]-u[1,2])/meas(edge)
                l = 2
                k = 1
                bm = fbernoulli(-(u[1,l]-u[1,k]))
                bp = fbernoulli((u[1,l]-u[1,k]))
                f[2]= -Jn*(bm*nc(u[1,k],u[2,k])-bp*nc(u[1,l],u[2,l]))/meas(edge)
                f[3]= Jp*(bp*pc(u[1,k],u[3,k])-bm*pc(u[1,l],u[3,l]))/meas(edge)
                #bpp,bmp=fbernoulli_pm(u[1,1]-u[1,-1])
                #f[3]= 0.0#Jp*(bp*pc(u[1,1],u[3,1])-bm*pc(u[1,2],u[3,2]))
        end,
        source=function (f,node)
            f[1] = C(node[1]*lscal)/Cscal
            f[2] = 0.0
            f[3] = 0.0
        end,
        reaction=function reaction!(f,u,node)
            f[1] = nc(u[1],u[2])-pc(u[1],u[3])
            f[2] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
            f[3] = R(nc(u[1],u[2])*Cscal,pc(u[1],u[3])*Cscal)/Rscal
        end)

    sys=VoronoiFVM.System(Gx,physics,unknown_storage=:sparse)

    # Add species 1 to region 1
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_species!(sys,3,[1])

    ul = 0.0+Va/Vscal
    ur = (Vbi+Vb)/Vscal
    unr = exp(-Va/Vscal)
    unl = exp(-Vb/Vscal)
    upr = exp(Va/Vscal)
    upl = exp(Vb/Vscal)
    # Set boundary conditions
    boundary_dirichlet!(sys,1,1,ul)
    boundary_dirichlet!(sys,1,2,ur)
    boundary_dirichlet!(sys,2,1,Va/Vscal)
    boundary_dirichlet!(sys,2,2,Vb/Vscal)
    boundary_dirichlet!(sys,3,1,Va/Vscal)
    boundary_dirichlet!(sys,3,2,Vb/Vscal)

    # Create a solution array
    inival=unknowns(sys,inival=0.001)
    print(size(inival))
    L = length(inival[1,:])÷2
    if isnothing(u0)
        inival[1,1:L] .= ul
        inival[1,L+1:2L] .= ur
    else
        inival[1,:] .= u0
    end

    solution=unknowns(sys)

    # Create solver control info

    # Stationary solution of the problem
    VoronoiFVM.solve!(solution,inival,sys, control=control)

    return X*lscal,solution[1,:]*Vscal,-log.(solution[2,:])*Vscal,log.(solution[3,:])*Vscal
end

function PoissonVoronoi(Ec,Ev,Nc,Nv,Vbi,Va,Vb,lscal,Ld0,L,λ,nld,nn,C,Vscal,Cscal)
    ## Selective Space Discretization
   ## Selective Space Discretization
   dx0 = 3*Ld0/L
   dxld = dx0/nld
   dxS0 = dx0/nn
   p0 = 0
   p1 = dx0
   p2 = 0.5-dx0
   p3 = 0.5+dx0
   p4 = 1.0-dx0
   p5 = 1.0
   X1 = p0:dxld:p1
   X2 = p1:dxS0:p2
   X3 = p2:dxld:p3
   X4 = p3:dxS0:p4
   X5 = p4:dxld:p5
   X = vcat(X1,X2[2:end],X3[2:end],X4[2:end],X5[2:end])*L/lscal
   Lint = maximum(X)
   ## Boltzman Dist
   F(x) = abs(x) > 500 ? exp(sign(x)*500) : exp(x)
   ηₙ(u) = u-Ec/Vscal
   nc(u) = Nc*F(ηₙ(u))/Cscal
   ηₚ(u) = -u+Ev/Vscal
   pc(u) = Nv*F(ηₚ(u))/Cscal

   Gx = VoronoiFVM.Grid(X)

   #Create PhysicalModel

   physics=VoronoiFVM.Physics(
       flux= function flux!(f,u,edge)
               f[1]=λ*(u[1,1]-u[1,2])
       end,
       source=function (f,node)
           f[1] = C(node[1]*lscal)/Cscal
       end,
       reaction=function reaction!(f,u,node)
           f[1] = nc(u[1])-pc(u[1])
       end)

   sys=VoronoiFVM.System(Gx,physics,unknown_storage=:sparse)

   # Add species 1 to region 1
   enable_species!(sys,1,[1])

   ul = 0.0+Va/Vscal
   ur = (Vbi+Vb)/Vscal

   # Set boundary conditions
   boundary_dirichlet!(sys,1,1,ul)
   boundary_dirichlet!(sys,1,2,ur)

   # Create a solution array
   inival=unknowns(sys,inival=1.0)
   print(size(inival))
   L = length(inival[1,:])÷2
   inival[1,1:L] .= ul
   inival[1,L+1:2L] .= ur

   solution=unknowns(sys)

   # Create solver control info
   control=VoronoiFVM.NewtonControl()
   #control.max_lureuse = 0
   control.verbose=true
   control.max_iterations = 500
   #control.Δu_opt=100
   control.damp_initial=1e-2
   control.tol_absolute = 1e-4
   control.tol_relative = 1e-4
   control.damp_growth = 2.0
   # Stationary solution of the problem
   VoronoiFVM.solve!(solution,inival,sys, control=control)

    return X*lscal,solution[1,:]*Vscal
end


text:
"""
   LaguerreWeight(α)

is a quasi-vector representing `x^α * exp(-x)` on `0..Inf`.
"""
struct LaguerreWeight{T} <: Weight{T}
    α::T
end

LaguerreWeight{T}() where T = LaguerreWeight{T}(zero(T))
LaguerreWeight() = LaguerreWeight{Float64}()
axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),)
function getindex(w::LaguerreWeight, x::Number)
    x ∈ axes(w,1) || throw(BoundsError())
    x^w.α * exp(-x)
end

sum(L::LaguerreWeight{T}) where T = gamma(L.α + 1)

struct Laguerre{T} <: OrthogonalPolynomial{T}
    α::T
    Laguerre{T}(α) where T = new{T}(convert(T, α))
end
Laguerre{T}() where T = Laguerre{T}(zero(T))
Laguerre() = Laguerre{Float64}()
Laguerre(α::T) where T = Laguerre{float(T)}(α)
orthogonalityweight(L::Laguerre)= LaguerreWeight(L.α)

==(L1::Laguerre, L2::Laguerre) = L1.α == L2.α
axes(::Laguerre{T}) where T = (Inclusion(HalfLine{T}()), oneto(∞))

"""
     laguerrel(n, α, z)

computes the `n`-th generalized Laguerre polynomial, orthogonal with
respec to `x^α * exp(-x)`, at `z`.
"""
laguerrel(n::Integer, α, z::Number) = Base.unsafe_getindex(Laguerre{promote_type(typeof(α), typeof(z))}(α), z, n+1)

"""
     laguerrel(n, z)

computes the `n`-th Laguerre polynomial, orthogonal with
respec to `exp(-x)`, at `z`.
"""
laguerrel(n::Integer, z::Number) = laguerrel(n, 0, z)

# L_{n+1} = (-1/(n+1) x + (2n+α+1)/(n+1)) L_n - (n+α)/(n+1) L_{n-1}
# - (n+α) L_{n-1} + (2n+α+1)* L_n -(n+1) L_{n+1} = x  L_n
# x*[L_0 L_1 L_2 …] = [L_0 L_1 L_2 …] * [(α+1)    -(α+1); -1  (α+3)     -(α+2);0  -2   (α+5) -(α+3); …]
function jacobimatrix(L::Laguerre{T}) where T
    α = L.α
    Tridiagonal(-(1:∞), (α+1):2:∞, -(α+1:∞))
end

recurrencecoefficients(L::Laguerre{T}) where T = ((-one(T)) ./ (1:∞), ((L.α+1):2:∞) ./ (1:∞), (L.α:∞) ./ (1:∞))

##########
# Derivatives
##########

@simplify function *(D::Derivative, L::Laguerre)
    T = promote_type(eltype(D),eltype(L))
    D = _BandedMatrix(Fill(-one(T),1,∞), ∞, -1,1)
    Laguerre(L.α+1)*D
end

ext: jl


text:
# Quante.jl
## A fork of PqQuante in Julia; experimenting with quantum chemistry in the Julia language.
module Quante

# Pull in necessary STDLIB stuff for Julia 0.7+
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
	using Test # so we can have tests distributed in source
	using Distributed
	using SharedArrays
	using LinearAlgebra
	using SpecialFunctions
end

include("physicalconstants.jl")

include("utility.jl")
include("gamma.jl")
include("basis.jl")
include("atomsmolecules.jl")

include("basissets.jl")

include("overlap.jl")
include("nuclear.jl")
include("kinetic.jl")

include("2eints.jl")
include("hrr.jl")

include("rhf.jl")

end # module


text:
export cwt, icwt, wct, wxs

"""

    wxs(ref, cur, t, window, freqbands, dj, unwrapflag)

Compute wavelet cross spectrum for two signals and a given array of frequency bands

# Arguments
- `ref::AbstractArray`: Reference time series
- `cur::AbstractArray`: Current time series
- `t::AbstractArray`: Time vector common to ref and cur
- `window::AbstractArray`: Times over which to compute dv/v
- `fs::Float64`: Sampling frequency
- `freqbands::AbstractArray`: Frequency bands over which to compute dv/v
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12
- `unwrapflag::Bool`: whether or not to unwrap phase

# Returns
- `freqbands::AbstractArray`: fmin and fmax for each iteration of the dv/v algorithm
- `dvv::Float64`: dv/v for current correlation
- `dvv_err::Float64`: Error for calculation of dv/v
- `int::Float64`: Intercept for regression calculation
- `int_err::Float64`: Error for intercept
- `dvv0::Float64`: dv/v for current correlation forced through the origin
- `dvv0_err::Float64`: Error for calculation of dvv0
"""
function wxs(ref::AbstractArray, cur::AbstractArray, t::AbstractArray, window::AbstractArray, fs::Float64, freqbands::AbstractArray; dj::AbstractFloat=1/12, unwrapflag::Bool=false)
    # define sample interval
    dt = 1/fs

    # perform wavelet coherence transform
    WXS, WXA, WCT, aWCT, coi, freqs = wct(cur, ref, dt, dj, freqbands)

    # do inverse cwt for different frequency bands
    if unwrapflag==true
        phase = unwrap(aWCT, dims=ndims(aWCT))
    else
        phase = aWCT
    end

    # if a frequency window is given (instead of a set of frequency bands), we assume
    # dv/v should be calculated for each frequency. We construct a 2D array of the
    # form [f1 f1; f2 f2; ...], which can be treated the same as a 2D array of frequency bands
    if ndims(freqbands)==1
        freqbands = hcat(freqs, freqs)
    end
    # number of frequency bands
    (nbands,_) = size(freqbands)

    # dt vector will be filled by regression of phase/frequency, and then
    # will be used to find dv/v by regression of dt/t
    delta_t = zeros(nbands, length(window))
    delta_t_err = zeros(nbands, length(window))

    dvv = zeros(nbands)
    dvv_err = similar(dvv)
    int = similar(dvv)
    int_err = similar(dvv)
    dvv0 = similar(dvv)
    dvv0_err = similar(dvv)

    # iterate over frequency bands
    for iband=1:nbands
        (fmin, fmax) = freqbands[iband, :]
        # frequency checks
        if fmax < fmin
            println("Error: please make sure columns 1 and 2 are the correct frequency limits in freqbands!")
        end
        freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)
        iphase = phase[freq_ind, window]

        # get dt by regression of phase delays/frequency band
        for itime=1:length(window)
            if fmin==fmax
                # simple quotient instead of regression, since we have only 1 point
                delta_t[iband, itime] = iphase[itime]/(2π*fmax)
            else
                # get weights
                w = 1 ./ WCT[freq_ind, itime]
                infNaN = findall(x->(isnan.(x) || isinf.(x)), w)
                if length(infNaN)!=0
                    w[infNaN] .= 1.0
                end
                # WLS inversion
                # This does NOT force the best fit line through the origin
                model = glm(@formula(Y ~ X),DataFrame(X=freqs[freq_ind]*2π,Y=iphase[:,itime]),Normal(),IdentityLink(),wts=w)
                delta_t[iband, itime] = coef(model)[2]
                delta_t_err[iband, itime] = stderror(model)[2]
            end
        end

        # regression in time to get time shift
        w2 = 1 ./ mean(WCT[freq_ind, window], dims=1)
        infNaN =findall(x->(isnan.(x) || isinf.(x)), w2[:])
        if length(infNaN)!=0
            w2[:, infNaN] .= 1.0
        end

        # find slope of dt/t to find -dv/v
        dvv[iband], dvv_err[iband], int[iband], int_err[iband], dvv0[iband], dvv0_err[iband] = dvv_lstsq(t[window], delta_t[iband, :], w=w2[:])
    end

    return freqbands, dvv, dvv_err, int, int_err, dvv0, dvv0_err
end

"""

  cwt(signal,dt,freqmin,freqmax)

Continuous wavelet transform of the signal at specified scales.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `signal::AbstractArray`: Input signal array.
- `dt::AbstractFloat`: Sampling interval [s].
- `freqmin::AbstractFloat`: Minimum frequency for cwt [Hz].
- `freqmax::AbstractFloat`: Maximum frequency for cwt [Hz].
- `f0::Real`: Nondimensional frequency from Torrence & Campo, 1998 eq. 1 [Hz].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `W::AbstractArray`: Wavelet transform from Morlet wavelet.
- `freqs::AbstractArray`: Fourier frequencies for wavelet scales [Hz].
- `coi::AbstractArray`: Cone of influence - maximum period (in s) of useful
    information at that particular time. Periods greater than
    those are subject to edge effects.

This is a Julia translation of the cwt in pycwt https://github.com/regeirk/pycwt
"""
function cwt(signal::AbstractArray{T,1},dt::AbstractFloat,freqmin::AbstractFloat,
             freqmax::AbstractFloat;f0=6.,dj=1/12) where T <: AbstractFloat
    n0 = length(signal)
    flambda = T.((4 .* π) ./ (f0 .+ sqrt(2 .+ f0^2)))
    s0 = 2 * dt / flambda
    J = convert(Int,round(log2(n0 .* dt ./ s0) ./ dj))

    # The scales as of Mallat 1999
    sj = T.(s0 .* 2 .^ ((0:J) .* dj))
    freqs = 1 ./ (flambda .* sj)

    # subset by freqmin & freqmax
    ind = findall((freqs .>= freqmin) .& (freqs .<= freqmax))
    sj = sj[ind]
    freqs = freqs[ind]

    # signal fft
    signal_ft = fft(signal,1)
    N = length(signal_ft)
    # Fourier angular frequencies
    ftfreqs = T.(2 .* π * FFTW.fftfreq(n0,1/dt))

    # Creates wavelet transform matrix as outer product of scaled transformed
    # wavelets and transformed signal according to the convolution theorem.
    # (i)   Transform scales to column vector for outer product;
    # (ii)  Calculate 2D matrix [s, f] for each scale s and Fourier angular
    #       frequency f;
    # (iii) Calculate wavelet transform;
    psi_ft_bar = ((sj' .* ftfreqs[2] .* T(N)) .^ T(0.5)) .* conj.(psi_ft(ftfreqs * sj',f0))
    W = ifft(signal_ft .* psi_ft_bar,1)

    # Checks for NaN in transform results and removes them from the scales if
    # needed, frequencies and wavelet transform. Trims wavelet transform at
    # length `n0`.
    sel = findall(.!all(isnan.(W),dims=1)[:])
    if length(sel) > 0
        sj = sj[sel]
        freqs = freqs[sel]
        W = W[:,sel]
    end

    # Determines the cone-of-influence. Note that it is returned as a function
    # of time in Fourier periods. Uses triangualr Bartlett window with
    # non-zero end-points.
    coi = n0 ./ 2 .- abs.((0:n0-1) .- (n0 - 1.) ./2.)
    coi = T.(flambda ./ sqrt(2.) .* dt .* coi)
    return W, sj, freqs, coi
end

"""

  icwt(W,sj,dt)

Inverse continuous wavelet transform at specified scales.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `W::AbstractArray`: Wavelet transform, the result of the `cwt` function.
- `sj::AbstractArray`: Scale indices as returned by the `cwt` function.
- `dt::AbstractFloat`: Sampling interval [s].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `iW::AbstractArray`: Inverse wavelet transform from Morlet wavelet.

This is a Julia translation of the icwt in pycwt https://github.com/regeirk/pycwt
Note that the pycwt version has incorrect scaling.
"""
function icwt(W::AbstractArray, sj::AbstractArray, dt::AbstractFloat;dj=1/12)
    T = real(eltype(W))
    # As of Torrence and Compo (1998), eq. (11)
    iW = T(dj .* sqrt(dt) ./ 0.776 .* (pi ^ 0.25)) .* sum(real.(W) ./ (sj .^ T(0.5))',dims=2)
    return iW
end

"""

  icwt(W,sj,dt)

Inverse continuous wavelet transform at specified scale.
Note: uses the Morlet wavelet ONLY.

# Arguments
- `W::AbstractArray`: Wavelet transform, the result of the `cwt` function.
- `sj::AbstractFloat`: Scale index as returned by the `cwt` function.
- `dt::AbstractFloat`: Sampling interval [s].
- `dj::AbstractFloat`: Spacing between discrete scales. Default value is 1/12.

# Returns
- `iW::AbstractArray`: Inverse wavelet transform from Morlet wavelet.

This is a Julia translation of the icwt in pycwt https://github.com/regeirk/pycwt
Note that the pycwt version has incorrect scaling.
"""
function icwt(W::AbstractArray, sj::AbstractFloat, dt::AbstractFloat;dj=1/12)
    T = real(eltype(W))
    # As of Torrence and Compo (1998), eq. (11)
    iW = T(dj .* sqrt(dt) ./ 0.776 .* (pi ^ 0.25)) .* real.(W) ./ sj ^ T(0.5)
    return iW
end

function psi_ft(A::AbstractArray{T},f0::Real) where T <: AbstractFloat
    return exp.(T(-0.5) .* (A .- T(f0)) .^2) .* T(π ^ -0.25)
end

"""

    smooth(W, dt, dj, scales)

Smooth wavelet spectrum.

# Arguments
- `W::AbstractArray`: wavelet spectrum
- `dt::Float64`: sampling interval in time
- `dj::Float64`: spacing between discrete scales
- `scales::AbstractArray`: wavelet scales

# Returns
- `T::AbstractArray`: Smoothed wavelet spectrum
"""
function smooth(W::AbstractArray, dt::Float64, dj::Float64, scales::AbstractArray)
    # The smoothing is performed by using a filter given by the absolute value
    # of the wavelet function at each scale, normalized to have a total weight
    # of unity, as per suggestions by Torrence &W ebster (1999) and by Grinsted et al. (2004).
    (m, n) = size(W)

    # Filter in time
    k = 2π*FFTW.fftfreq(length(W[1,:]))
    k2 = k.^2
    snorm = scales ./ dt

    # Smoothing by Gaussian window (absolute value of wavelet function)
    # using the convolution theorem: multiplication by Gaussian curve in
    # Fourier domain for each scale, outer product of scale and frequency
    F = exp.(-0.5 .* (snorm.^2) .* k2') # outer product
    smooth = (ifft(F .* fft(W,2), 2))

    T = smooth[:, 1:n] # Remove possible padded region due to FFTW

    # Filter in scale. For the Morlet wavelet, this is simply a boxcar with 0.6 width
    # construct boxcar
    wsize = convert(Int64, round(0.60 / dj * 2))
    if wsize % 2 == 0
        wsize+=1
    end
    halfWin = div(wsize,2)

    # iterate over 'horizontal' and smooth in the 'vertical'
    # this could also be done by adding an axis to the transpose and performing a 2d convolution
    for i=1:size(T,2)
        # pad signal for 'same' padding
        paddedT = vcat(T[1,i]*ones(halfWin), T[:,i], T[end,i]*ones(halfWin))
        # smooth signal
        win = Array{eltype(paddedT), 1}(undef, wsize)
        win.=ones(wsize)/wsize
        smoothT = conv(paddedT, win)
        # trim signal
        T[:,i] = smoothT[2*halfWin+1:end-2*halfWin]
    end

    return T
end

"""

    wct(y1, y2, dt, dj, freqbands)

Wavelet coherence transform, which finds regions in the wavelet domain where the two time
series co-vary but don't necessarily have high power.

# Arguments
- `y1::AbstractArray`: input signal
- `y2::AbstractArray`: input signal
- `dt::AbstractFloat`: sampling interval in time
- `dj::AbstractFloat`: spacing between discrete scales
- `freqbands::AbstractArray`: Frequency bands over which to compute dv/v
- `f0::Real`: Nondimensional frequency from Torrence & Campo, 1998 eq. 1 [Hz].
- `norm::Bool`: Whether or not to normalize signals before cwt

# Returns
- `WXS::AbstractArray`: Wavelet cross-spectrum
- `WXA::AbstractArray`: Amplitude of wavelet cross-spectrum
- `rWCT::AbstractArray`: Real part of the wavelet coherence transform
- `aWCT::AbstractArray`: Angle of wavelet coherence transform
- `coi::AbstractArray`: Cone of influence - maximum period (in s) of useful
    information at that particular time. Periods greater than
    those are subject to edge effects.
- `freqs::AbstractArray`: Fourier frequencies for wavelet scales [Hz].

Original from pycwt and modified from Python implementation by Congcong Yuan
"""
function wct(y1::AbstractArray, y2::AbstractArray, dt::AbstractFloat, dj::AbstractFloat, freqbands::AbstractArray; f0=6.,norm::Bool=true)
    # normalize signals
    if norm
        y1 = (y1 .- mean(y1)) ./ std(y1)
        y2 = (y2 .- mean(y2)) ./ std(y2)
    end

    # calculate the CWT of the time series, using identical parameters for both calculations
    W1, sj, freqs, coi = cwt(y1, dt, minimum(freqbands), maximum(freqbands))
    W2, sj, freqs, coi = cwt(y2, dt, minimum(freqbands), maximum(freqbands))
    W1=W1'
    W2=W2'

    scales = Array{Float64,2}(undef, size(W1))
    for i=1:size(scales, 2)
        scales[:,i] .= sj
    end

    # smooth wavelet spectra before truncating
    S1 = smooth(abs.(W1).^2 ./ scales, dt, dj, sj)
    S2 = smooth(abs.(W2).^2 ./ scales, dt, dj, sj)

    # compute cross wavelet transform
    W12 = W1 .* conj.(W2)

    S12 = smooth(W12 ./ scales, dt, dj, sj)
    WCT = abs.(S12).^2 ./ (S1 .* S2)
    aWCT = angle.(W12)

    # calculate cross spectrum and its amplitude
    WXS = W12
    WXA = abs.(S12)

    return WXS, WXA, real.(WCT), aWCT, coi, freqs
end


text:
using YaoBenchmarks
using Test

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


text:
using LinearAlgebra
using QuantumOptics
using WignerSymbols

function Kfactor(F,I=1.5,J=0.5)
    return F*(F+1)-I*(I+1)-J*(J+1)
end

function LandeFactor(J,L=0,S=0.5)
    return 1+(J*(J+1)+S*(S+1)-L*(L+1))/(2*J*(J+1))
end

struct CGmatrix{T<:Rational}
    cgMatrix::Array{Float64}
    function CGmatrix{T}(j1::T,j2::T) where T<:Rational
        if !(( j1.den == 1 || j1.den == 2 ) && ( j2.den == 1  || j2.den ==2 ))
            error("spins are not half-integer")
        end

        row=1
        col=1
        dims = Int64((j1*2+1) * (j2*2+1))
        cgMatrix=zeros(Float64,dims,dims)
        j3max=maximum((abs(j1+j2), abs(j1-j2)))
        j3min=minimum((abs(j1+j2), abs(j1-j2)))
        for j3=j3max:-1:j3min
            for m3=j3:-1:-j3
                for m2=j2:-1:-j2
                    for m1=j1:-1:-j1
                        cgMatrix[row,col]=Float64(clebschgordan(j1,m1,j2,m2,j3,m3))
                        col+=1
                    end

                end
                row+=1
                col=1
            end
        end
        new(cgMatrix)
    end
end
CGmatrix(a::T,b::T) where T<: Rational = CGmatrix{T}(a,b)

function getCGmatrix(self::CGmatrix)
    return self.cgMatrix
end

struct HyperfineManifold{T<:Rational}
    fmBaseKets::Array{Ket}
    fvalues::Array{Float64}
    function HyperfineManifold{T}(spinI::T,spinJ::T) where T<:Rational
        if !(( spinI.den == 1 || spinI.den == 2 ) && ( spinJ.den == 1  || spinJ.den ==2 ))
            error("spins are not half-integer")
        end

        mIbasis=SpinBasis(spinI)
        mJbasis=SpinBasis(spinJ)
        mImJbasis=CompositeBasis(mIbasis,mJbasis)

        mIbasisSz= 0.5*sigmaz(mIbasis)               # Sz operator for S1 basis: 0.5 \hbar [1 0; 0 -1]
        mJbasisSz = 0.5*sigmaz(mJbasis)               # Sz operator for S2 basis: 0.5 \hbar [1 0; 0 -1]
        mIbasisS₊ = sigmap(mIbasis)                # Splus operator for S1 basis: \hbar [0 1; 0  0]
        mIbasisS₋ = sigmam(mIbasis)               # Sminus operator for S1 basis: \hbar [0 0; 1  0]
        mJbasisS₊ = sigmap(mJbasis)                # Splus operator for S2 basis: \hbar [0 1; 0  0]
        mJbasisS₋ = sigmam(mJbasis)               # Sminus operator for S1 basis: \hbar [0 0; 1  0]
        mIbasisIdentity = identityoperator(mIbasis)   # Identity operator for S1 basis: [1 0; 0 1]
        mJbasisIdentity = identityoperator(mJbasis)   # Identity operator for S2 basis: [1 0; 0 1]
        mIbasisS² = spinI*(spinI+1)*mIbasisIdentity              # S^2 operator for S1 basis: I(I+1) \hbar^2 [1 0; 0 1]
        mJbasisS² = spinJ*(spinJ+1)*mJbasisIdentity              # S^2 operator for S2 basis: J(J+1) \hbar^2 [1 0; 0 1]

        Jz=tensor(mIbasisIdentity,mJbasisSz)     # S1z in S1S2 basis
        Iz=tensor(mIbasisSz,mJbasisIdentity)     # S2z in S1S2 basis
        J₊=tensor(mIbasisIdentity,mJbasisS₊)     # S1plus in S1S2 basis
        J₋=tensor(mIbasisIdentity,mJbasisS₋)     # S1minus in S1S2 basis
        I₊=tensor(mIbasisS₊,mJbasisIdentity)     # S2plus in S1S2 basis
        I₋=tensor(mIbasisS₋,mJbasisIdentity)     # S2minus in S1S2 basis
        mImJbasisI²=tensor(mIbasisS²,mJbasisIdentity)
        mImJbasisJ²=tensor(mIbasisIdentity,mJbasisS²)
        IJ=Iz*Jz+0.5*(I₊*J₋+I₋*J₊)        # S1̇ S2 in S1S2 basis
        mImJbasisF²= mImJbasisI² + mImJbasisJ² + 2*Iz*Jz +I₊*J₋ +I₋*J₊

        C=getCGmatrix(CGmatrix(spinI,spinJ))

        TransformMatrix=DenseOperator(mImJbasis,C)
        TransformMatrixDagger=dagger(TransformMatrix)

        Izfmbasis=TransformMatrix*(Iz)*TransformMatrixDagger    # Eq. E-10 in Chapter XII.E.2.a. of vol. 2 of Cohen-Tannoudji, Quantum Mechanics.
        Jzfmbasis=TransformMatrix*(Jz)*TransformMatrixDagger

        K1=Kfactor(1)
        K2=Kfactor(2)
        FMhfsDiagonal=Diagonal([K2, K2, K2, K2, K2, K1, K1, K1])
        IJfmbasis2=0.5*DenseOperator(mImJbasis,FMhfsDiagonal)

        fmBaseKets=Array{Ket}(undef,length(mImJbasis))
        fvalues=Array{Rational}(undef, length(mImJbasis))
        mvalues=Array{Rational}(undef, length(mImJbasis))

        for i = 1:length(mImJbasis)
            fmBaseKets[i] = Ket(mImJbasis,C[i,:])
            mᵢ = real(dagger(fmBaseKets[i])*Iz*fmBaseKets[i])
            mⱼ = real(dagger(fmBaseKets[i])*Jz*fmBaseKets[i])
            mf = mᵢ+mⱼ
            F² = real(dagger(fmBaseKets[i])*mImJbasisF²*fmBaseKets[i])
            F = (-1 + sqrt(1+4*F²))/2
            fvalues[i] = Int64(round(2*F))//2
            mvalues[i] = Int64(round(2*mf))//2
        end

        new(fmBaseKets,fvalues)
    end
end
HyperfineManifold(a::T,b::T) where T<: Rational = HyperfineManifold{T}(a,b)

max_stars_repo_name: LU-LaserCenter/LC-Julia-Notebooks

max_issues_repo_name: LU-LaserCenter/LC-Julia-Notebooks

max_forks_repo_name: LU-LaserCenter/LC-Julia-Notebooks


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ram
using AWS.Compat
using AWS.UUIDs

"""
    AcceptResourceShareInvitation()

Accepts an invitation to a resource share from another AWS account.

# Required Parameters
- `resourceShareInvitationArn`: The Amazon Resource Name (ARN) of the invitation.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
"""
accept_resource_share_invitation(resourceShareInvitationArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/acceptresourceshareinvitation", Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn); aws_config=aws_config)
accept_resource_share_invitation(resourceShareInvitationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/acceptresourceshareinvitation", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn), args)); aws_config=aws_config)

"""
    AssociateResourceShare()

Associates the specified resource share with the specified principals and resources.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `principals`: The principals.
- `resourceArns`: The Amazon Resource Names (ARN) of the resources.
"""
associate_resource_share(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/associateresourceshare", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
associate_resource_share(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/associateresourceshare", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    AssociateResourceSharePermission()

Associates a permission with a resource share.

# Required Parameters
- `permissionArn`: The ARN of the AWS RAM permission to associate with the resource share.
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `replace`: Indicates whether the permission should replace the permissions that are currently associated with the resource share. Use true to replace the current permissions. Use false to add the permission to the current permission.
"""
associate_resource_share_permission(permissionArn, resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/associateresourcesharepermission", Dict{String, Any}("permissionArn"=>permissionArn, "resourceShareArn"=>resourceShareArn); aws_config=aws_config)
associate_resource_share_permission(permissionArn, resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/associateresourcesharepermission", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("permissionArn"=>permissionArn, "resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    CreateResourceShare()

Creates a resource share.

# Required Parameters
- `name`: The name of the resource share.

# Optional Parameters
- `allowExternalPrincipals`: Indicates whether principals outside your AWS organization can be associated with a resource share.
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `permissionArns`: The ARNs of the permissions to associate with the resource share. If you do not specify an ARN for the permission, AWS RAM automatically attaches the default version of the permission for each resource type.
- `principals`: The principals to associate with the resource share. The possible values are IDs of AWS accounts, the ARN of an OU or organization from AWS Organizations.
- `resourceArns`: The Amazon Resource Names (ARN) of the resources to associate with the resource share.
- `tags`: One or more tags.
"""
create_resource_share(name; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/createresourceshare", Dict{String, Any}("name"=>name); aws_config=aws_config)
create_resource_share(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/createresourceshare", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("name"=>name), args)); aws_config=aws_config)

"""
    DeleteResourceShare()

Deletes the specified resource share.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
"""
delete_resource_share(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("DELETE", "/deleteresourceshare", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
delete_resource_share(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("DELETE", "/deleteresourceshare", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    DisassociateResourceShare()

Disassociates the specified principals or resources from the specified resource share.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `principals`: The principals.
- `resourceArns`: The Amazon Resource Names (ARNs) of the resources.
"""
disassociate_resource_share(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/disassociateresourceshare", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
disassociate_resource_share(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/disassociateresourceshare", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    DisassociateResourceSharePermission()

Disassociates an AWS RAM permission from a resource share.

# Required Parameters
- `permissionArn`: The ARN of the permission to disassociate from the resource share.
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
"""
disassociate_resource_share_permission(permissionArn, resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/disassociateresourcesharepermission", Dict{String, Any}("permissionArn"=>permissionArn, "resourceShareArn"=>resourceShareArn); aws_config=aws_config)
disassociate_resource_share_permission(permissionArn, resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/disassociateresourcesharepermission", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("permissionArn"=>permissionArn, "resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    EnableSharingWithAwsOrganization()

Enables resource sharing within your AWS Organization. The caller must be the master account for the AWS Organization.

"""
enable_sharing_with_aws_organization(; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/enablesharingwithawsorganization"; aws_config=aws_config)
enable_sharing_with_aws_organization(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/enablesharingwithawsorganization", args; aws_config=aws_config)

"""
    GetPermission()

Gets the contents of an AWS RAM permission in JSON format.

# Required Parameters
- `permissionArn`: The ARN of the permission.

# Optional Parameters
- `permissionVersion`: The identifier for the version of the permission.
"""
get_permission(permissionArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getpermission", Dict{String, Any}("permissionArn"=>permissionArn); aws_config=aws_config)
get_permission(permissionArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getpermission", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("permissionArn"=>permissionArn), args)); aws_config=aws_config)

"""
    GetResourcePolicies()

Gets the policies for the specified resources that you own and have shared.

# Required Parameters
- `resourceArns`: The Amazon Resource Names (ARN) of the resources.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `principal`: The principal.
"""
get_resource_policies(resourceArns; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourcepolicies", Dict{String, Any}("resourceArns"=>resourceArns); aws_config=aws_config)
get_resource_policies(resourceArns, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourcepolicies", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArns"=>resourceArns), args)); aws_config=aws_config)

"""
    GetResourceShareAssociations()

Gets the resources or principals for the resource shares that you own.

# Required Parameters
- `associationType`: The association type. Specify PRINCIPAL to list the principals that are associated with the specified resource share. Specify RESOURCE to list the resources that are associated with the specified resource share.

# Optional Parameters
- `associationStatus`: The association status.
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `principal`: The principal. You cannot specify this parameter if the association type is RESOURCE.
- `resourceArn`: The Amazon Resource Name (ARN) of the resource. You cannot specify this parameter if the association type is PRINCIPAL.
- `resourceShareArns`: The Amazon Resource Names (ARN) of the resource shares.
"""
get_resource_share_associations(associationType; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshareassociations", Dict{String, Any}("associationType"=>associationType); aws_config=aws_config)
get_resource_share_associations(associationType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshareassociations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("associationType"=>associationType), args)); aws_config=aws_config)

"""
    GetResourceShareInvitations()

Gets the invitations for resource sharing that you've received.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `resourceShareArns`: The Amazon Resource Names (ARN) of the resource shares.
- `resourceShareInvitationArns`: The Amazon Resource Names (ARN) of the invitations.
"""
get_resource_share_invitations(; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshareinvitations"; aws_config=aws_config)
get_resource_share_invitations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshareinvitations", args; aws_config=aws_config)

"""
    GetResourceShares()

Gets the resource shares that you own or the resource shares that are shared with you.

# Required Parameters
- `resourceOwner`: The type of owner.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `name`: The name of the resource share.
- `nextToken`: The token for the next page of results.
- `resourceShareArns`: The Amazon Resource Names (ARN) of the resource shares.
- `resourceShareStatus`: The status of the resource share.
- `tagFilters`: One or more tag filters.
"""
get_resource_shares(resourceOwner; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshares", Dict{String, Any}("resourceOwner"=>resourceOwner); aws_config=aws_config)
get_resource_shares(resourceOwner, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/getresourceshares", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceOwner"=>resourceOwner), args)); aws_config=aws_config)

"""
    ListPendingInvitationResources()

Lists the resources in a resource share that is shared with you but that the invitation is still pending for.

# Required Parameters
- `resourceShareInvitationArn`: The Amazon Resource Name (ARN) of the invitation.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
"""
list_pending_invitation_resources(resourceShareInvitationArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listpendinginvitationresources", Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn); aws_config=aws_config)
list_pending_invitation_resources(resourceShareInvitationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listpendinginvitationresources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn), args)); aws_config=aws_config)

"""
    ListPermissions()

Lists the AWS RAM permissions.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `resourceType`: Specifies the resource type for which to list permissions. For example, to list only permissions that apply to EC2 subnets, specify ec2:Subnet.
"""
list_permissions(; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listpermissions"; aws_config=aws_config)
list_permissions(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listpermissions", args; aws_config=aws_config)

"""
    ListPrincipals()

Lists the principals that you have shared resources with or that have shared resources with you.

# Required Parameters
- `resourceOwner`: The type of owner.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `principals`: The principals.
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.
- `resourceShareArns`: The Amazon Resource Names (ARN) of the resource shares.
- `resourceType`: The resource type. Valid values: codebuild:Project | codebuild:ReportGroup | ec2:CapacityReservation | ec2:DedicatedHost | ec2:Subnet | ec2:TrafficMirrorTarget | ec2:TransitGateway | imagebuilder:Component | imagebuilder:Image | imagebuilder:ImageRecipe | license-manager:LicenseConfiguration I resource-groups:Group | rds:Cluster | route53resolver:ResolverRule
"""
list_principals(resourceOwner; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listprincipals", Dict{String, Any}("resourceOwner"=>resourceOwner); aws_config=aws_config)
list_principals(resourceOwner, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listprincipals", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceOwner"=>resourceOwner), args)); aws_config=aws_config)

"""
    ListResourceSharePermissions()

Lists the AWS RAM permissions that are associated with a resource share.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
"""
list_resource_share_permissions(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresourcesharepermissions", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
list_resource_share_permissions(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresourcesharepermissions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    ListResourceTypes()

Lists the shareable resource types supported by AWS RAM.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
"""
list_resource_types(; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresourcetypes"; aws_config=aws_config)
list_resource_types(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresourcetypes", args; aws_config=aws_config)

"""
    ListResources()

Lists the resources that you added to a resource shares or the resources that are shared with you.

# Required Parameters
- `resourceOwner`: The type of owner.

# Optional Parameters
- `maxResults`: The maximum number of results to return with a single call. To retrieve the remaining results, make another call with the returned nextToken value.
- `nextToken`: The token for the next page of results.
- `principal`: The principal.
- `resourceArns`: The Amazon Resource Names (ARN) of the resources.
- `resourceShareArns`: The Amazon Resource Names (ARN) of the resource shares.
- `resourceType`: The resource type. Valid values: codebuild:Project | codebuild:ReportGroup | ec2:CapacityReservation | ec2:DedicatedHost | ec2:Subnet | ec2:TrafficMirrorTarget | ec2:TransitGateway | imagebuilder:Component | imagebuilder:Image | imagebuilder:ImageRecipe | license-manager:LicenseConfiguration I resource-groups:Group | rds:Cluster | route53resolver:ResolverRule
"""
list_resources(resourceOwner; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresources", Dict{String, Any}("resourceOwner"=>resourceOwner); aws_config=aws_config)
list_resources(resourceOwner, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/listresources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceOwner"=>resourceOwner), args)); aws_config=aws_config)

"""
    PromoteResourceShareCreatedFromPolicy()

Resource shares that were created by attaching a policy to a resource are visible only to the resource share owner, and the resource share cannot be modified in AWS RAM. Use this API action to promote the resource share. When you promote the resource share, it becomes:   Visible to all principals that it is shared with.   Modifiable in AWS RAM.

# Required Parameters
- `resourceShareArn`: The ARN of the resource share to promote.

"""
promote_resource_share_created_from_policy(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/promoteresourcesharecreatedfrompolicy", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
promote_resource_share_created_from_policy(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/promoteresourcesharecreatedfrompolicy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)

"""
    RejectResourceShareInvitation()

Rejects an invitation to a resource share from another AWS account.

# Required Parameters
- `resourceShareInvitationArn`: The Amazon Resource Name (ARN) of the invitation.

# Optional Parameters
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
"""
reject_resource_share_invitation(resourceShareInvitationArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/rejectresourceshareinvitation", Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn); aws_config=aws_config)
reject_resource_share_invitation(resourceShareInvitationArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/rejectresourceshareinvitation", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareInvitationArn"=>resourceShareInvitationArn), args)); aws_config=aws_config)

"""
    TagResource()

Adds the specified tags to the specified resource share that you own.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.
- `tags`: One or more tags.

"""
tag_resource(resourceShareArn, tags; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/tagresource", Dict{String, Any}("resourceShareArn"=>resourceShareArn, "tags"=>tags); aws_config=aws_config)
tag_resource(resourceShareArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/tagresource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn, "tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes the specified tags from the specified resource share that you own.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.
- `tagKeys`: The tag keys of the tags to remove.

"""
untag_resource(resourceShareArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/untagresource", Dict{String, Any}("resourceShareArn"=>resourceShareArn, "tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceShareArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/untagresource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn, "tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateResourceShare()

Updates the specified resource share that you own.

# Required Parameters
- `resourceShareArn`: The Amazon Resource Name (ARN) of the resource share.

# Optional Parameters
- `allowExternalPrincipals`: Indicates whether principals outside your AWS organization can be associated with a resource share.
- `clientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `name`: The name of the resource share.
"""
update_resource_share(resourceShareArn; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/updateresourceshare", Dict{String, Any}("resourceShareArn"=>resourceShareArn); aws_config=aws_config)
update_resource_share(resourceShareArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ram("POST", "/updateresourceshare", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceShareArn"=>resourceShareArn), args)); aws_config=aws_config)


text:
# MIT License

# Copyright (c) 2021 Fran Urbano
# Copyright (c) 2021 Val Lyashov
# Copyright (c) 2020 Bryan Palmer

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
Package: Forecast

    hma(x, n)

Applies the Henderson moving average filter to dataset `x` with `n`-term.

"Henderson moving averages are filters which were derived by Robert Henderson in 1916
for use in actuarial applications. They are trend filters, commonly used in time series
analysis to smooth seasonally adjusted estimates in order to generate a trend estimate.

They are used in preference to simpler moving averages because they can reproduce
polynomials of up to degree 3, thereby capturing trend turning points.

The ABS uses Henderson moving averages to produce trend estimates from a seasonally
adjusted series. The trend estimates published by the ABS are typically derived using
a 13 term Henderson filter for monthly series, and a 7 term Henderson filter for quarterly series.

Henderson filters can be either symmetric or asymmetric. Symmetric moving averages can be applied
at points which are sufficiently far away from the ends of a time series. In this case, the
smoothed value for a given point in the time series is calculated from an equal number of values
on either side of the data point." - Australian Bureau of Statistics (www.abs.gov.au)

# Arguments
- `x`: Observations' support.
- `n`: Observation values, tipically 13 or 7 for quarterly data.

# Returns
An array of Henderson filter smoothed values provided in `s`.

# Examples
```julia-repl
julia> hma(rand(1000), 303)
1000-element Vector{Float64}}:
[...]
```
"""
@inline function hma(x::AbstractVector{<:T}, n::Integer) where T<:Real

    if !isodd(n)
        @warn "Henderson moving averages only accepts odds values for 'n';\nreturning instead the average of hma(x,"*string(n)*"-1) and hma(x,"*string(n)*"+1)"

        hma1 = hma(x,n-1)
        hma2 = hma(x,n+1)
        return (hma1 + hma2) / 2
    end

    @assert n >= 5 "n must be greater or equal 5"
    @assert (lx = length(x)) >= n "length(s) must be greater or equal to n"

    # Promotion to Float type
    P = Base.promote_op(/,T,T)
    x,_ = Base.promote(collect(x),[P(1.0)])

    w = hmaSymmetricWeights(n,P)

    m = (n-1) ÷ 2

    ic = n < 13 ? P(1.0) : (13 <= n < 15 ? P(3.5) : P(4.5))
    b2s2 = P(4.0)/P(pi)/ic^2

    function hmai(i)
        if i - 1 < m
            u = hmaAsymmetricWeights(m + i, w, b2s2)[end:-1:1]
            sum(x[1:i + m] .* u)
        elseif i - 1 + m >= lx
            u = hmaAsymmetricWeights(m + lx - i + 1, w, b2s2)
            sum(x[i-m:lx] .* u)
        else
            sum(x[i-m:i+m] .* w)
        end
    end

    hmav = similar(x)
    for i in 1:lx
        @inbounds hmav[i] = hmai(i)
    end
    hmav

end

"""
package: Smoothers

    hmaSymmetricWeights(n,P)

Caluclate the hma symmetric weights for 'n'

# Arguments
- `n`: Number of symmetric weights

# Returns

 Vector of symmetrical hma weights with type related to the bit size of `n`

# Refenrences
- "A Guide to Interpreting Time Series" ABS (2003), page 41.
"""
@inline function hmaSymmetricWeights(n::Integer,P::DataType)

    m = (P(n)-1)÷2
    m1 = (m+1)^2
    m2 = (m+2)^2
    m3 = (m+3)^2
    logd = log(m+2)+log(m2-1)+log(4*m2-1)+log(4*m2-9)+log(4*m2-25)
    d = 315/(8*exp(logd))

    w = Vector{P}(undef,Int(m)+2)
    for (i,v) in enumerate(0:m+1)
         @inbounds w[i] = real(exp(log(m1-v^2)+log(m2-v^2)+log(m3-v^2)+
                               log(complex(P(3.0)*m2-P(11.0)*v^2-P(16.0)))+log(d)))
    end

    u = vcat(w[end-1:-1:1],w[2:end-1])

    return mod(n, 2) != 0 ? u : vcat(u, missing)

end

"""
package: Smoothers

    hmaAsymmetricWeights(m,w)

Calculate the hma asymmetric end-weights for 'm' given 'w'

# Arguments
- `m`: Number of asymmetric weights (m < length(w))
- `w`: Vector of symmetrical hma weights

# Returns

Vector of asymmetrical hma weights

# References

- Mike Doherty (2001), "The Surrogate Henderson Filters in X-11",Aust, NZ J of Stat. 43(4), 2001, 385–392
"""
@inline function hmaAsymmetricWeights(m::Integer, w::AbstractVector{<:T}, b2s2::T) where T<:Real

    n = length(w)
    # @assert m <= n "The m argument must be less than w"
    # @assert m >= (n-1)÷2 "The m argument must be greater than (n-1)/2"

    sumResidual = sum(w[range(m + 1, n, step = 1)]) / T(m)

    sumEnd = T(0.0)
    m12 = (m+T(1.0))/T(2.0)
    for i in m+1:n
        @inbounds sumEnd += (i-m12)*w[i]
    end

    denominator = T(1.0) + ((m*(m-T(1.0))*(m+T(1.0)) / T(12.0) ) * b2s2)

    aw = Vector{T}(undef,m)
    for r in 1:m
        numerator = (T(r) - (m+T(1.0)) / T(2.0)) * b2s2
        @inbounds aw[r] = w[r] + sumResidual +  numerator / denominator * sumEnd
    end
    aw
end


text:
import HypothesisTests.OneSampleADTest

"""
    OneSampleADTest(uv::UncertainValue, d::UnivariateDistribution,
        n::Int = 1000) -> OneSampleADTest

Perform a one-sample Anderson–Darling test of the null hypothesis that
a draw of `n` realisations of the uncertain value `uv` comes from the
distribution `d` against the alternative hypothesis that the sample is not
drawn from `d`.
"""
function OneSampleADTest(uv::AbstractUncertainValue, d::UnivariateDistribution, n::Int = 1000)
    x = resample(uv, n)
    OneSampleADTest(x, d)
end

"""
    OneSampleADTestPooled(ud::UncertainDataset, d::UnivariateDistribution,
        n::Int = 1000)) -> OneSampleADTest

First, draw `n` realisations of each uncertain value in `ud` and pool them
together. Then perform a one-sample Anderson–Darling test of the null
hypothesis that the pooled values comes from the distribution `d` against the
alternative hypothesis that the sample is not drawn from `d`.
"""
function OneSampleADTestPooled(ud::UncertainDataset, d::UnivariateDistribution, n::Int = 1000)
    x = vcat(resample(ud, n)...,)
    OneSampleADTest(x, d)
end

"""
    OneSampleADTestElementWise(ud::UncertainDataset, d::UnivariateDistribution,
        n::Int = 1000)) -> Vector{OneSampleADTest}
First, draw `n` realisations of each uncertain value in `ud`, keeping one
pool of values for each uncertain value.
Then, perform an element-wise (pool-wise) one-sample Anderson–Darling test of
the null hypothesis that each value pool comes from the distribution `d`
against the alternative hypothesis that the sample is not drawn from `d`.
"""
function OneSampleADTestElementWise(ud::UncertainDataset, d::UnivariateDistribution, n::Int = 1000)
    [OneSampleADTest(resample(ud[i], n), d) for i = 1:length(ud)]
end

export
OneSampleADTest,
OneSampleADTestPooled,
OneSampleADTestElementWise


text:
export BrickletRealTimeClockV2DateTime
struct BrickletRealTimeClockV2DateTime
    year::Integer
    month::Integer
    day::Integer
    hour::Integer
    minute::Integer
    second::Integer
    centisecond::Integer
    weekday::Integer
    timestamp::Integer
end

export BrickletRealTimeClockV2Alarm
struct BrickletRealTimeClockV2Alarm
    month::Integer
    day::Integer
    hour::Integer
    minute::Integer
    second::Integer
    weekday::Integer
    interval::Integer
end

export BrickletRealTimeClockV2SPITFPErrorCount
struct BrickletRealTimeClockV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletRealTimeClockV2Identity
struct BrickletRealTimeClockV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletRealTimeClockV2
"""
Battery-backed real-time clock
"""
mutable struct BrickletRealTimeClockV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletRealTimeClockV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_real_time_clock_v2")
        deviceInternal = package.BrickletRealTimeClockV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_date_time
"""
    $(SIGNATURES)

Sets the current date (including weekday) and the current time.

If the backup battery is installed then the real-time clock keeps date and
time even if the Bricklet is not powered by a Brick.

The real-time clock handles leap year and inserts the 29th of February
accordingly. But leap seconds, time zones and daylight saving time are not
handled.
"""
function set_date_time(device::BrickletRealTimeClockV2, year, month, day, hour, minute, second, centisecond, weekday)
    device.deviceInternal.set_date_time(year, month, day, hour, minute, second, centisecond, weekday)
end

export get_date_time
"""
    $(SIGNATURES)

Returns the current date (including weekday) and the current time of the
real-time.

The timestamp represents the current date and the the current time of the
real-time clock converted to milliseconds and is an offset to 2000-01-01 00:00:00.0000.
"""
function get_date_time(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_date_time()
end

export get_timestamp
"""
    $(SIGNATURES)

Returns the current date and the time of the real-time clock converted to
milliseconds. The timestamp has an effective resolution of hundredths of a
second and is an offset to 2000-01-01 00:00:00.0000.
"""
function get_timestamp(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_timestamp()
end

export set_offset
"""
    $(SIGNATURES)

Sets the offset the real-time clock should compensate for in 2.17 ppm steps
between -277.76 ppm (-128) and +275.59 ppm (127).

The real-time clock time can deviate from the actual time due to the frequency
deviation of its 32.768 kHz crystal. Even without compensation (factory
default) the resulting time deviation should be at most ±20 ppm (±52.6
seconds per month).

This deviation can be calculated by comparing the same duration measured by the
real-time clock (``rtc_duration``) an accurate reference clock
(``ref_duration``).

For best results the configured offset should be set to 0 ppm first and then a
duration of at least 6 hours should be measured.

The new offset (``new_offset``) can be calculated from the currently configured
offset (``current_offset``) and the measured durations as follow::

  new_offset = current_offset - round(1000000 * (rtc_duration - ref_duration) / rtc_duration / 2.17)

If you want to calculate the offset, then we recommend using the calibration
dialog in Brick Viewer, instead of doing it manually.

The offset is saved in the EEPROM of the Bricklet and only needs to be
configured once.
"""
function set_offset(device::BrickletRealTimeClockV2, offset)
    device.deviceInternal.set_offset(offset)
end

export get_offset
"""
    $(SIGNATURES)

Returns the offset as set by :func:`Set Offset`.
"""
function get_offset(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_offset()
end

export set_date_time_callback_configuration
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Date Time` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_date_time_callback_configuration(device::BrickletRealTimeClockV2, period)
    device.deviceInternal.set_date_time_callback_configuration(period)
end

export get_date_time_callback_configuration
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Date Time Callback Configuration`.
"""
function get_date_time_callback_configuration(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_date_time_callback_configuration()
end

export set_alarm
"""
    $(SIGNATURES)

Configures a repeatable alarm. The :cb:`Alarm` callback is triggered if the
current date and time matches the configured alarm.

Setting a parameter to -1 means that it should be disabled and doesn't take part
in the match. Setting all parameters to -1 disables the alarm completely.

For example, to make the alarm trigger every day at 7:30 AM it can be
configured as (-1, -1, 7, 30, -1, -1, -1). The hour is set to match 7 and the
minute is set to match 30. The alarm is triggered if all enabled parameters
match.

The interval has a special role. It allows to make the alarm reconfigure itself.
This is useful if you need a repeated alarm that cannot be expressed by matching
the current date and time. For example, to make the alarm trigger every 23
seconds it can be configured as (-1, -1, -1, -1, -1, -1, 23). Internally the
Bricklet will take the current date and time, add 23 seconds to it and set the
result as its alarm. The first alarm will be triggered 23 seconds after the
call. Because the interval is not -1, the Bricklet will do the same again
internally, take the current date and time, add 23 seconds to it and set that
as its alarm. This results in a repeated alarm that triggers every 23 seconds.

The interval can also be used in combination with the other parameters. For
example, configuring the alarm as (-1, -1, 7, 30, -1, -1, 300) results in an
alarm that triggers every day at 7:30 AM and is then repeated every 5 minutes.
"""
function set_alarm(device::BrickletRealTimeClockV2, month, day, hour, minute, second, weekday, interval)
    device.deviceInternal.set_alarm(month, day, hour, minute, second, weekday, interval)
end

export get_alarm
"""
    $(SIGNATURES)

Returns the alarm configuration as set by :func:`Set Alarm`.
"""
function get_alarm(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_alarm()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletRealTimeClockV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletRealTimeClockV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletRealTimeClockV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletRealTimeClockV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletRealTimeClockV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletRealTimeClockV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletRealTimeClockV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletRealTimeClockV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletRealTimeClockV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
module SatSolver

include("formula.jl")

export Clause, Formula, is_formula_satisfiable

include("nnf_parser/nnf_parser.jl")

export parse_nnf_formula

include("io.jl")

export print_dimacs_cnf

include("unit_propagation.jl")

include("brute_force.jl")
export brute_force

include("dpll.jl")
export dpll

end # module


text:
module gmm

using InferGMM, Distributions
using PyPlot
import ..utils: gaussian_2D_level_curve_pts
export plot_gmm

function plot_gmm(dGMM::GMM; figsize=(10,10), alpha=0.5, bins=50, fill=false, axs=nothing) where T <: AbstractFloat
    d, k = size(dGMM), ncomponents(dGMM)
    (d > 20) && throw("will not plot for d > 20")

    # d == 2 fits on a single axis: deal with this case first
    if d == 2
        if axs == nothing
            fig, axs = PyPlot.subplots(1, 1, figsize=figsize)
        end
        _plot = fill ? axs.fill : axs.plot
        for j in 1:k
            levcurv = gaussian_2D_level_curve_pts(dGMM.mus[j,:], dGMM.sigmas[:,:,j])
            _plot(levcurv[:,1], levcurv[:,2], alpha=alpha*dGMM.pis[j]/maximum(dGMM.pis))
        end
        return axs
    end

    # d > 2:
    if axs == nothing
        fig, axs = PyPlot.subplots(d, d, figsize=figsize)
    end
    rsmp = rand(dGMM, 5000)

    for ix = 1:d, iy = 1:d
        if ix != iy
            for j in 1:k
                _plot = fill ? axs[iy, ix].fill : axs[iy, ix].plot
                levcurv = gaussian_2D_level_curve_pts(dGMM.mus[j,:][[ix,iy]], dGMM.sigmas[[ix,iy],[ix,iy],j])
                _plot(levcurv[:,1], levcurv[:,2], alpha=alpha*dGMM.pis[j]/maximum(dGMM.pis))
            end
        else
            axs[iy, ix].hist(rsmp[ix, :], bins=bins)
        end
    end
    return axs
end

end

max_stars_repo_name: ornithos/AxPlot

max_issues_repo_name: ornithos/AxPlot

max_forks_repo_name: ornithos/AxPlot


text:
"""
    AbstractStochasticStructure{N}

Abstract supertype for the underlying memory structure of a stochastic program. `N` is the number of stages.
"""
abstract type AbstractStochasticStructure{N} end
"""
    StochasticInstantiation

Abstract supertype for the underlying memory structure of a stochastic program. `N` is the number of stages.
"""
abstract type StochasticInstantiation end
"""
    UnspecifiedInstantiation

Default instantiation value, which defers the choice to `default_structure`.

See also: [`default_structure`](@ref)
"""
struct UnspecifiedInstantiation <: StochasticInstantiation end
"""
    Deterministic

Instantiates with the `DeterministicEquivalent` structure.

See also: [`DeterministicEquivalent`](@ref)
"""
struct Deterministic <: StochasticInstantiation end
"""
    Vertical

Instantiates with the `VerticalStructure` on a single core.

See also: [`VerticalStructure`](@ref)
"""
struct Vertical <: StochasticInstantiation end
"""
    Horizontal

Instantiates with the `HorizontalStructure` on a single core.

See also: [`HorizontalStructure`](@ref)
"""
struct Horizontal <: StochasticInstantiation end
"""
    DistributedVertical

Instantiates with the `VerticalStructure` on multiple cores.

See also: [`VerticalStructure`](@ref)
"""
struct DistributedVertical <: StochasticInstantiation end
"""
    DistributedHorizontal

Instantiates with the `HorizontalStructure` on multiple cores.

See also: [`HorizontalStructure`](@ref)
"""
struct DistributedHorizontal <: StochasticInstantiation end
"""
    StochasticStructure(scenario_types::ScenarioTypes{M}, instantiation::StochasticInstantiation) where M

Constructs a stochastic structure over the `M` provided scenario types according to the specified `instantiation`. Should be overrided for every defined stochastic structure.

    StochasticStructure(scenarios::NTuple{M, Vector{<:AbstractScenario}}, instantiation::StochasticInstantiation) where M

Constructs a stochastic structure over the `M` provided scenario sets according to the specified `instantiation`. Should be overrided for every defined stochastic structure.
"""
function StochasticStructure end
"""
    default_structure(instantiation::StochasticInstantiation, optimizer)

Returns a `StochasticInstantiation` based on the provided `instantiation` and `optimizer`. If an explicit `instantiation` is provided it is always prioritized. Otherwise, if `instantiation` is `UnspecifiedInstantiation`, returns whatever structure requested by the optimizer. Defaults to `Deterministic` if no optimizer is provided.
"""
function default_structure(instantiation::StochasticInstantiation, optimizer)
    return instantiation
end
function default_structure(::UnspecifiedInstantiation, optimizer)
    if optimizer isa MOI.AbstractOptimizer
        if optimizer isa AbstractStructuredOptimizer
            # default to vertical structure
            if nworkers() > 1
                # Distribute in memory if Julia processes are available
                return DistributedVertical()
            else
                return Vertical()
            end
        else
            # Default to DEP structure if standard MOI optimizer is given
            return Deterministic()
        end
    end
    # Default to deterministic if no recognized optimizer has been given.
    return Deterministic()
end

# Optimization #
# ========================== #
"""
    UnsupportedStructure{Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure}

Error indicating that an optimizer of type `Opt` does not support the stochastic structure `S`.
"""
struct UnsupportedStructure{Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure} <: Exception end

function Base.showerror(io::IO, err::UnsupportedStructure{Opt, S}) where {Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure}
    print(io, "The stochastic structure $S is not supported by the optimizer $Opt")
end

"""
    UnloadedStructure{Opt <: StochasticProgramOptimizerType}

Error thrown when an optimizer of type `Opt` has not yet loaded a stochastic structure and an operation which requires a structure to be loaded is called.
"""
struct UnloadedStructure{Opt <: StochasticProgramOptimizerType} <: Exception end

function Base.showerror(io::IO, err::UnloadedStructure{Opt}) where Opt <: StochasticProgramOptimizerType
    print(io, "The optimizer $Opt has no loaded structure. Consider `load_structure!`")
end

"""
    UnloadasbleStructure{Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure}

Error thrown when an optimizer of type `Opt` cannot load a structure of type `S`.
"""
struct UnloadableStructure{Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure} <: Exception
    message # Human-friendly explanation why the structure cannot be loaded by the optimizer
end

function Base.showerror(io::IO, err::UnloadableStructure{Opt, S}) where {Opt <: StochasticProgramOptimizerType, S <: AbstractStochasticStructure}
    print(io, "The optimizer $Opt cannot load structure $S")
    m = message(err)
    if Base.isempty(m)
        print(io, ".")
    else
        print(io, ": ", m)
    end
end
message(err::UnloadableStructure) = err.message

# Getters #
# ========================== #
function structure_name(structure::AbstractStochasticStructure)
    return "Unknown"
end
function num_decisions(structure::AbstractStochasticStructure{N}, stage::Integer = 1) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return num_decisions(structure.decisions[stage])
end
function scenario_type(structure::AbstractStochasticStructure, s::Integer = 2)
    return _scenario_type(scenarios(structure, s))
end
function _scenario_type(::Vector{S}) where S <: AbstractScenario
    return S
end
function num_scenarios(structure::AbstractStochasticStructure, stage::Integer = 2)
    return length(scenarios(structure, stage))
end
function probability(structure::AbstractStochasticStructure, stage::Integer, scenario_index::Integer)
    return probability(scenario(structure, stage, scenario_index))
end
function stage_probability(structure::AbstractStochasticStructure, stage::Integer = 2)
    return probability(scenarios(structure, stage))
end
function expected(structure::AbstractStochasticStructure, stage::Integer = 2)
    return expected(scenarios(structure, stage))
end
function distributed(structure::AbstractStochasticStructure, stage::Integer)
    return false
end
function DecisionRef(proxy::JuMP.Model, structure::AbstractStochasticStructure, index::VI, at_stage::Integer, scenario_index::Integer)
    at_stage > 1 || error("There are no scenarios in the first at_stage.")
    n = num_scenarios(structure, at_stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return DecisionRef(proxy, index)
end
# ========================== #

# Printing #
# ========================== #
function _print(io::IO, ::AbstractStochasticStructure)
    # Just give summary as default
    show(io, stochasticprogram)
end
# ========================== #