text:
using NeuralPDE
using Quadrature, Cubature, Cuba
using Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
using DelimitedFiles
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

function level_set(strategy, minimizer, maxIters)

    ##  DECLARATIONS
    @parameters  t x y
    @variables   u(..)

    Dt = Differential(t)
    Dx = Differential(x)
    Dy = Differential(y)

    # Discretization
    xwidth      = 1.0      #ft
    ywidth      = 1.0
    tmax        = 1.0      #min
    xScale      = 1.0
    yScale      = 1.0
    xMeshNum    = 10
    yMeshNum    = 10
    tMeshNum    = 10
    dx  = xwidth/xMeshNum
    dy  = ywidth/yMeshNum
    dt  = tmax/tMeshNum

    domains = [t ∈ Interval(0.0,tmax),
               x ∈ Interval(0.0,xwidth),
               y ∈ Interval(0.0,ywidth)]

    xs = 0.0 : dx : xwidth
    ys = 0.0 : dy : ywidth
    ts = 0.0 : dt : tmax

    # Definitions
    x0    = 0.5
    y0    = 0.5
    Uwind = [0.0, 2.0]  #wind vector

    # Operators
    gn   = (Dx(u(t,x,y))^2 + Dy(u(t,x,y))^2)^0.5  #gradient's norm
    ∇u   = [Dx(u(t,x,y)), Dy(u(t,x,y))]
    n    = ∇u/gn              #normal versor
    #U    = ((Uwind[1]*n[1] + Uwind[2]*n[2])^2)^0.5 #inner product between wind and normal vector

    R0 = 0.112471
    ϕw = 0#0.156927*max((0.44*U)^0.04086,1.447799)
    ϕs = 0
    S  = R0*(1 + ϕw + ϕs)

    # Equation
    eq = Dt(u(t,x,y)) + S*gn ~ 0  #LEVEL SET EQUATION

    initialCondition = (xScale*(x - x0)^2 + (yScale*(y - y0)^2))^0.5 - 0.2   #Distance from ignition

    bcs = [u(0,x,y) ~ initialCondition]  #from literature

    ## NEURAL NETWORK
    n = 10   #neuron number

    chain = FastChain(FastDense(3,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))   #Neural network from Flux library

    indvars = [t,x,y]   #phisically independent variables
    depvars = [u(t,x,y)]       #dependent (target) variable

    dim = length(domains)

    losses = []
    error = []
    times = []

    dx_err = 0.1

    error_strategy = GridTraining(dx_err)

    initθ = Float64.(DiffEqFlux.initial_params(chain))
    eltypeθ = eltype(initθ)
    parameterless_type_θ = DiffEqBase.parameterless_type(initθ)

    phi = NeuralPDE.get_phi(chain,parameterless_type_θ)
    derivative = NeuralPDE.get_numeric_derivative()

    _pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,
                                             phi,derivative,nothing,chain,initθ,error_strategy)

    bc_indvars = NeuralPDE.get_variables(bcs,indvars,depvars)
    _bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,
                                              phi,derivative,nothing,chain,initθ,error_strategy,
                                              bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    train_sets = NeuralPDE.generate_training_sets(domains,dx_err,[eq],bcs,eltypeθ,indvars,depvars)
    train_domain_set, train_bound_set = train_sets

    pde_loss_functions = [NeuralPDE.get_loss_function(_pde_loss_function,
                                          train_domain_set[1],eltypeθ,
                                          parameterless_type_θ,error_strategy)]

    bc_loss_functions = [NeuralPDE.get_loss_function(bc,
                                          train_set,eltypeθ,
                                          parameterless_type_θ,error_strategy) for (train_set,bc) in zip(train_bound_set, _bc_loss_functions)]
    loss_functions =  [pde_loss_functions;bc_loss_functions]
    function loss_function_(θ,p)
        return sum(map(l->l(θ) ,loss_functions))
    end

    cb_ = function (p,l)
        deltaT_s = time_ns() #Start a clock when the callback begins, this will evaluate questo misurerà anche il calcolo degli uniform error

        ctime = time_ns() - startTime - timeCounter #This variable is the time to use for the time benchmark plot
        append!(times, ctime/10^9) #Conversion nanosec to seconds
        append!(losses, l)
        loss_ = loss_function_(p,nothing)
        append!(error, loss_)
        timeCounter = timeCounter + time_ns() - deltaT_s #timeCounter sums all delays due to the callback functions of the previous iterations

        return false
    end

    @named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)

    discretization = NeuralPDE.PhysicsInformedNN(chain, strategy; init_params =initθ)
    prob = NeuralPDE.discretize(pde_system,discretization)

    timeCounter = 0.0
    startTime = time_ns() #Fix initial time (t=0) before starting the training
    res = GalacticOptim.solve(prob, minimizer, cb = cb_, maxiters=maxIters)

    phi = discretization.phi

    params = res.minimizer

    # Model prediction
    domain = [ts, xs, ys]

    u_predict = [reshape([first(phi([t,x,y],res.minimizer)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts]  #matrix of model's prediction

    return [error, params, domain, times, losses] #add numeric solution
end

#level_set(NeuralPDE.QuadratureTraining(algorithm = CubaCuhre(), reltol = 1e-8, abstol = 1e-8, maxiters = 100), ADAM(0.01), 500)

maxIters = [(1,1,1,1000,1000,1000,1000),(1,1,1,500,500,500,500)] #iters for ADAM/LBFGS
# maxIters = [(1,1,1,2,2,2,2),(1,1,1,2,2,2,2)] #iters for ADAM/LBFGS

strategies = [NeuralPDE.QuadratureTraining(quadrature_alg = CubaCuhre(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.QuadratureTraining(quadrature_alg = HCubatureJL(), reltol = 1, abstol = 1e-4, maxiters = 100, batch = 0),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLp(), reltol = 1, abstol = 1e-4, maxiters = 100),
              NeuralPDE.GridTraining(0.1),
              NeuralPDE.StochasticTraining(400 ; bcs_points= 50),
              NeuralPDE.QuasiRandomTraining(400 ; bcs_points= 50)]

strategies_short_name = ["CubaCuhre",
                        "HCubatureJL",
                        "CubatureJLh",
                        "CubatureJLp",
                        "GridTraining",
                        "StochasticTraining",
                        "QuasiRandomTraining"]

minimizers = [ADAM(0.005),
              #BFGS()]
              LBFGS()]

minimizers_short_name = ["ADAM",
                         "LBFGS"]
                        # "BFGS"]

# Run models
prediction_res =  Dict()
error_res =  Dict()
domains = Dict()
params_res = Dict()  #to use same params for the next run
times = Dict()
losses_res = Dict()

## Convergence
for min =1:length(minimizers) # minimizer
      for strat=1:length(strategies) # strategy
            # println(string(strategies_short_name[strat], "  ", minimizers_short_name[min]))
            res = level_set(strategies[strat], minimizers[min], maxIters[min][strat])
            push!(error_res, string(strat,min)     => res[1])
            push!(params_res, string(strat,min) => res[2])
            push!(domains, string(strat,min)        => res[3])
            push!(times, string(strat,min)        => res[4])
            push!(losses_res, string(strat,min)        => res[5])
      end
end

#Plotting the first strategy with the first minimizer out from the loop to initialize the canvas
current_label = string(strategies_short_name[1], " + " , minimizers_short_name[1])
error = Plots.plot(times["11"], error_res["11"], yaxis=:log10, label = current_label)# xlims = (0,10))#legend = true)#, size=(1200,700))
plot!(error, times["21"], error_res["21"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[1]))
plot!(error, times["31"], error_res["31"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[1]))
plot!(error, times["41"], error_res["41"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[1]))
plot!(error, times["51"], error_res["51"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[1]))
plot!(error, times["61"], error_res["61"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[1]))
plot!(error, times["71"], error_res["71"], yaxis=:log10, label = string(strategies_short_name[7], " + " , minimizers_short_name[1]))

plot!(error, times["12"], error_res["12"], yaxis=:log10, label = string(strategies_short_name[1], " + " , minimizers_short_name[2]))
plot!(error, times["22"], error_res["22"], yaxis=:log10, label = string(strategies_short_name[2], " + " , minimizers_short_name[2]))
plot!(error, times["32"], error_res["32"], yaxis=:log10, label = string(strategies_short_name[3], " + " , minimizers_short_name[2]))
plot!(error, times["42"], error_res["42"], yaxis=:log10, label = string(strategies_short_name[4], " + " , minimizers_short_name[2]))
plot!(error, times["52"], error_res["52"], yaxis=:log10, label = string(strategies_short_name[5], " + " , minimizers_short_name[2]))
plot!(error, times["62"], error_res["62"], yaxis=:log10, label = string(strategies_short_name[6], " + " , minimizers_short_name[2]))
plot!(error, times["72"], error_res["72"], yaxis=:log10, title = string("Level Set convergence ADAM/LBFGS"), ylabel = "log(error)", xlabel = "t", label = string(strategies_short_name[7], " + " , minimizers_short_name[2]))

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
"""
    Semiconductor1D

    Represents a series on SemiconductorData objects with positions
"""
mutable struct Semiconductor1D{K} <: AbstractArray{SemiconductorPoint{K},1}
    points::AbstractArray{SemiconductorPoint{K},1}
end

Base.size(S::Semiconductor1D) = size(S.points)
Base.push!(S::Semiconductor1D,e::SemiconductorPoint{Float64}) = (S.points = vcat(S.points,[e]))
Base.length(S::Semiconductor1D) = length(S.points)
Base.getindex(S::Semiconductor1D,i) = S.points[i]
Base.lastindex(S::Semiconductor1D) = length(S)
Base.setindex!(S::Semiconductor1D,v,i) = (S.points[i] = v)

#Semiconductor1D(V::AbstractArray )= Semiconductor1D(V)
"""
    Semiconductor1D(S::SemiconductorData,r,psi)

    Create a Semiconductor1D with S as the base semiconductor at positions r and electrostatic potential psi (referenced at the Fermi Level) at positions r (in cm)

    ```jldoctest
    julia>
        T = 300 # K
        N_a = 5e17 # cm^-3
        r = LinRange(0,1e-4,101) # in cm
        psi = LinRange(0.0,1.0,101) # in V
        semi1D = Semiconductor1D(SemicondorData(T,BoltzmanDist(),PSilicon(N_a)),r,psi)
    ```
"""
Semiconductor1D(S::SemiconductorData,r,psi) = Semiconductor1D(SemiconductorPoint1D.(Ref(S),r,psi))

"""
    Semiconductor1D(S::SemiconductorData,r,psi,ϕn,ϕp)

    Create a Semiconductor1D with S as the base semiconductor at positions r and potentials psi, electron quasi-Fermi level ϕn, and holes quasi-Fermi level ϕn (referenced at the Fermi Level) at positions r
"""
Semiconductor1D(S::SemiconductorData,r,psi,ϕn,ϕp) = Semiconductor1D(SemiconductorPoint1D.(Ref(S),r,psi,ϕn,ϕp))

Semiconductor1D(S::SemiconductorData,r::AbstractVector{T}) where T <: Number = Semiconductor1D(SemiconductorPoint1D.(Ref(S),r,zeros(length(r))))

#Semiconductor1D() = Semiconductor1D(Vector{SemiconductorPoint{Float64}})
"""
    Eg(S::Semiconductor1D) Bandgap at different positions
"""
Eg(S::Semiconductor1D) = Eg.(S.points)
"""
    Ec(S::Semiconductor1D) Conduction Band at different positions referenced to the Fermi Level
"""
Ec(S::Semiconductor1D) = Ec.(S.points)

"""
    EV(S::Semiconductor1D) Valence Band at different positions referenced to the Fermi Level
"""
Ev(S::Semiconductor1D) = Ev.(S.points)
Ef(S::Semiconductor1D) = Ef.(S.points)
"""
    Ei(S::Semiconductor1D) Intrinsic Fermi Level at different positions referenced to the Fermi Level
"""
Ei(S::Semiconductor1D) = Ei.(S.points)
"""
    n(S::Semiconductor1D) Electrons concentration at different positions referenced to the Fermi Level
"""
n(S::Semiconductor1D) = n.(S.points)
"""
    p(S::Semiconductor1D) Holes concentration at different positions in cm^-3
"""
p(S::Semiconductor1D) = p.(S.points)
n(S::Semiconductor1D,ψ) = n.(S.points,ψ)
p(S::Semiconductor1D,ψ) = p.(S.points,ψ)
"""
    Na(S::Semiconductor1D) Acceptors concentration at different positions in cm^-3
"""
Na(S::Semiconductor1D) = Na.(S.points)
"""
    Nd(S::Semiconductor1D) Donors concentration at different positions in cm^-3
"""
Nd(S::Semiconductor1D) = Nd.(S.points)
N(S::Semiconductor1D) = N.(S.points)
N(S::Semiconductor1D,ψ) = N.(S.points,ψ)
P(S::Semiconductor1D,ψ) = P.(S.points,ψ)
P(S::Semiconductor1D) = P.(S.points)
"""
    ni(S::Semiconductor1D) Intrinsic carriers concentration at different positions in cm^-3
"""
ni(S::Semiconductor1D) = ni.(S.points)
"""
    φf(S::Semiconductor1D) Ef-Ei at different positions in V
"""
φf(S::Semiconductor1D) = φf.(S.points)
"""
    r(S::Semiconductor1D) Positions in which the 1D Semiconductor is defined in cm
"""
r(S::Semiconductor1D) = r.(S.points)
"""
    Rₗ(S::Semiconductor1D) End positions in which the 1D Semiconductor is defined in cm
"""
Rₗ(S::Semiconductor1D) = maximum(r(S))
"""
    R₀(S::Semiconductor1D) Begin positions in which the 1D Semiconductor is defined in cm
"""
R₀(S::Semiconductor1D) = minimum(r(S))
"""
    L(S::Semiconductor1D) Semiconductor1D length in cm
"""
L(S::Semiconductor1D) = Rₗ(S)-R₀(S)
"""
    ρ(S::Semiconductor1D) Semiconductor1D electric charge density
"""
ρ(S::Semiconductor1D) = ρ.(S.points)
#ρ(S::Semiconductor1D,ψ) = -q/ϵ₀/ϵᵣ(S[1]) *  (P(S,ψ).-N(S,ψ)) *S.W*S.W
"""
    ξ(S::Semiconductor1D) Semiconductor1D electric field calculated analytically
"""
ξ(S::Semiconductor1D) = ξ.(S.points)
ξ(S::Semiconductor1D,ψ) = ξ.(S.points,ψ)
ξ(S::Semiconductor1D,ψ,ϕ) = ξ.(S.points,ψ,ϕ)
"""
    ψ(S::Semiconductor1D) Electrostatic potential at the different positions
"""
ψ(S::Semiconductor1D) = ψ.(S.points)
Δχ(S::Semiconductor1D,δχ) = Δχ.(S.points,δχ)
"""
    Temperature(S::Semiconductor1D) Latice Temperature
"""
Temperature(S::Semiconductor1D) = Temperature(S.points[1])
"""
    ϵᵣ(S::Semiconductor1D) Relative permitivity
"""
ϵᵣ(S::Semiconductor1D)=ϵᵣ(S.points[1])
"""
    χ(S::Semiconductor1D) Work function
"""
χ(S::Semiconductor1D) = χ.(S.points)
"""
    Nv(S::Semiconductor1D) Effective density of states on the Valence Band
"""
Nv(S::Semiconductor1D) = Nv(S.points[1])
"""
    Nc(S::Semiconductor1D) Effective density of states on the Conduction Band
"""
Nc(S::Semiconductor1D) = Nc(S.points[1])
Pt(S::Semiconductor1D) = Pt.(S.points)
Nt(S::Semiconductor1D) = Nt.(S.points)
"""
    DistFun(S::Semiconductor1D) Distribution Function of the Semiconductor1D
"""
DistFun(S::Semiconductor1D,E) = DistFun(S.points,E)
"""
    C(S::Semiconductor1D) Signed doping concentration of the semiconductor
"""
C(S::Semiconductor1D) = C.(S.points)
"""
    φfₑ(S::Semiconductor1D) Electron Quasi-Fermi level
"""
φfₑ(S::Semiconductor1D) = φfₑ.(S.points)
"""
    φfₚ(S::Semiconductor1D) holes Quasi-Fermi level
"""
φfₚ(S::Semiconductor1D) = φfₚ.(S.points)

"""
    translate!(S::Semiconductor1D,pt::Number) Translates the semiconductor slab by pt
"""
function translate!(S::Semiconductor1D,pt::Number)
    v = S.points
    for i in 1:length(v)
        S.points[i] = SemiconductorPoint1D(S[i],S[i].p+pt)
    end
end

function vjoint(L1D::Semiconductor1D,R1D::Semiconductor1D,dx=1e-9)
    if R₀(R1D) < Rₗ(L1D)
        translate!(R1D,Rₗ(L1D)+dx)
    end
    #Create Log Spaced Points on Transition
    lastL = pop!(L1D.points)
    firstR = popfirst!(R1D.points)
    midle = 0.5*(r(L1D[end])+r(R1D[1]))
    r1 = collect(10.0.^LinRange(log10(r(L1D[end])),log10(midle),11))
    popfirst!(r1)
    pop!(r1)
    #Left Side
    for r in r1
        push!(L1D,SemiconductorPoint1D(lastL.Data,r,ψ(lastL)))
    end

    r2 = collect(10.0.^LinRange(log10(midle),log10(r(R1D[1])),11))
    pop!(r2)
    popfirst!(r2)
    modelR = R1D[1]
    for r in reverse(r2)
        pushfirst!(R1D.points,SemiconductorPoint1D(firstR.Data,r,ψ(firstR)))
    end

    ΔE = (Ef(L1D[1])-Ef(R1D[1]))
    applyPsiProfile(R1D,ΔE.*ones(length(R1D)))
    #Δχ(R1D,-ΔE.*ones(length(R1D)))
    v = vcat(L1D.points,R1D.points)
    Semiconductor1D(v)
end

function vjointLin(L1D::Semiconductor1D,R1D::Semiconductor1D,dx=1e-9)
    if R₀(R1D) < Rₗ(L1D)
        translate!(R1D,Rₗ(L1D))
    end
    #Create Log Spaced Points on Transition
    lastL = pop!(L1D.points)
    #firstR = popfirst!(R1D.points)

    ΔE = (Ef(L1D[1])-Ef(R1D[1]))
    Δχ(R1D,ΔE)
    applyPsiProfile(R1D,-ΔE.*ones(length(R1D)))

    v = vcat(L1D.points,R1D.points)
    Semiconductor1D(v)
end

"""
    BandDiagram(S::Semiconductor1D) Plots the Semiconductor1D BandDiagram
"""
function BandDiagram(S::Semiconductor1D)
    BandDiagram(S,Rₗ(S))
end

function BandDiagram(S::Semiconductor1D,LM)
    y = 1e4*r(S)
	if(maximum(y) < LM)
        Sp = SemiconductorPoint(S[end].Data,LM,0.0)
		push!(S,Sp)
	end
    ψₛ = ψ(S[1])
    #plot(r(S),Ec(S)+χ(S),label=L"E_0")
    plot(r(S),Ec(S),label=L"E_C")
    plot!(r(S),Ev(S),label=L"E_V")
    plot!(r(S),Ef(S)-ψ(S),label=L"\varPhi+E_F")
    plot!(r(S),Ei(S),label=L"E_i")

    LM = L(S)
    return plot!(label=L"E_F",xlims=(0,LM),xlabel=L"r\, [\mu m]",legend=:best)#,yticks=nothing)

end

function BandDiagramJ(S::Semiconductor1D)
    BandDiagramJ(S,Rₗ(S))
end

function BandDiagramJ(S::Semiconductor1D,LM)
    #plot(r(S),Ec(S)+χ(S),label=L"E_0")
    plot(1e4*r(S),Ec(S),label=L"E_C")
    plot!(1e4*r(S),Ev(S),label=L"E_V")
    #plot!(1e4*r(S),Ei(S),label=L"\varPhi+E_F")
    plot!(1e4*r(S),-ψ(S)+Ei(S),label=L"E_i")

    #plot!(r(S),Ei(S),label=L"E_i")

    LM = L(S)*1e4
    return plot!(label=L"E_F",xlims=(0,LM),xlabel=L"r\, [\mu m]",legend=:best)#,yticks=nothing)

end

function CarrierInjectionPlot(S::Semiconductor1D,LM = 5.0)
    y = r(prob)
    LM = maximum(y)
    Ny = N(prob)
    Py = P(prob)
    niy = ni(prob)
    plot(y,Ny./niy,label=L"N/n_i")
    plot!(y,Py./niy,label=L"P/n_i",xlims=(-0.1,LM),yaxis=:log10,xlabel=L"r\, [\mu m]")
end

"""
    SingleSlabBVProblem(ψₛ,S::SemiconductorMaterial,f::DistributionFunc,T::Number) Solves by integrating the analyticall 1/ξ(x) the uniform semiconductor SingleSlabBVProblem and return a SemiconductorPoint1D
"""
function SingleSlabBVProblem(ψₛ,S::SemiconductorMaterial,f::DistributionFunc,T::Number)
    ## Precalculate the Material Properties
    Sd = SemiconductorData(T,f,S)
    if ψₛ == 0
        y = collect(LinRange(0,150,101))
        psi = 0*y
        return Semiconductor1D(Sd,y,psi)
    end
    psi1 = LinRange(ψₛ,0.5*ψₛ,51)
    psi2 = 10.0.^LinRange(log10(abs(ψₛ*0.5)),log10(abs(ψₛ*1e-3)),101)
    if ψₛ < 0
        psi2 = -1.0 .* psi2
    end
    psi = vcat(psi1,psi2)
    SdV = [SemiconductorData(Sd,p*Eunit) for p in psi]
    integrand_y(ψ) = 1/ξ(Sd,ψ)
    Y(ψ) = quadgk(integrand_y,ψ,ψₛ)[1]
    L = abs.(Y.(psi))
    return Semiconductor1D(return Semiconductor1D(Sd,L,psi))
end

function relaxationCostFunction(S::Semiconductor1D)
    Eest = -∇ρ.(Ref(S),r(S))
    error = sqrt(sum((ξ(S).-Eest).^2))
end

"""
    applyPsiProfile(S::Semiconductor1D,ψ::T) where T <: AbstractArray

    Apply an electrostatic potential to the semiconductor 1D
"""
function applyPsiProfile(S::Semiconductor1D,ψ::T) where T <: AbstractArray
    if length(S) != length(ψ)
        throw(DomainError(ψ,"ψ must have the same number of points as spatial points as S"))
    end
    for i in 1:length(ψ)
        S[i] = applyPsi(S[i],ψ[i])
    end
end

"""
    applyEfeProfile(S::Semiconductor1D,φfₑ::AbstractArray)

    Apply an electron quasi-fermi level to the semiconductor 1D
"""
function applyEfeProfile(S::Semiconductor1D,φfₑ::AbstractArray)
    if length(S) != length(φfₑ)
        throw(DomainError(ψ,"ψ must have the same number of points as spatial points as S"))
    end
    for i in 1:length(φfₑ)
        S[i] = applyEfe(S[i],φfₑ[i])
    end
end
"""
    applyEfpProfile(S::Semiconductor1D,φfₚ::AbstractArray)

    Apply an holes quasi-fermi level to the semiconductor 1D
"""
function applyEfpProfile(S::Semiconductor1D,φfₚ::AbstractArray)
    if length(S) != length(φfₚ)
        throw(DomainError(φfₚ,"ψ must have the same number of points as spatial points as S"))
    end
    for i in 1:length(φfₚ)
        S[i] = applyEfp(S[i],φfₚ[i])
    end
end

function ApplyPsiAtPoint(S::Semiconductor1D,ψ,x)
    v = r(S)
    if x in v
        i = findlast(v .== x)
        S[i] = applyPsi(S[i],ψ)
    else
        m = findlast( v .< x )
        M = m + 1
        δr = v[M]-v[m]
        δM = v[M]-x
        δm = x-v[m]
        S[M] = applyPsi(S[M],δM/δr *ψ)
        S[m] = applyPsi(S[m],δm/δr *ψ)
    end
end

ext: jl


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, GeoDatasets
CairoMakie.activate!() # HIDE
let
    lon, lat, data = GeoDatasets.landseamask(;resolution='c',grid=10)
    fig = Figure(resolution = (1250,700), fontsize = 25)
    ax = Axis(fig, xlabel = "lon", ylabel = "lat", xtickalign = 1,
        ytickalign = 1, xticksize=15, yticksize=15, spinewidth = 0.75)
    cmap = cgrad(:Greys_4, 3, categorical = true, rev = false)
    pltobj = heatmap!(ax, lon, lat, data, colormap = cmap)
    cbar   = Colorbar(fig, pltobj, width = 11,ticklabelsize = 25,tickalign =1,
        spinewidth = 0.5, ticks = ([0.35,1,1.65],["Ocean", "Land", "Lakes"]),
        ticklabelrotation = Quaternionf0(0,0,0,0)) #
    limits!(ax, round(lon[1]), round(lon[end]), round(lat[1]), round(lat[end]))
    hidespines!(ax, :t, :r)
    ax.xticks = -180:60:180
    ax.yticks = -90:30:90
    ax.xtickformat = "{:d}ᵒ"
    ax.ytickformat = "{:d}ᵒ"
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "landSea.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["CairoMakie","GeoDatasets"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
include("clonetask.jl")
include("brokentask.jl")
include("tarray.jl")
include("tarray2.jl")
include("tarray3.jl")


text:
__precompile__(true)
## Julia Parallel Computing documentation: https://docs.julialang.org/en/v1/manual/parallel-computing/index.html
module ArettoClusterABMLAB
    using Distributed
    import Distributed.launch, Distributed.manage, Distributed.kill, Distributed.init_worker, Distributed.connect
    export ArettoManager
    export launch, manage

    #worker_arg = `--worker`

    #function __init__()
    #    global worker_arg
    #    worker_arg = `--worker=$(Distributed.cluster_cookie())`
    #end

    worker_arg() = `--worker=$(Distributed.init_multi(); cluster_cookie())`

    struct ArettoManager <: ClusterManager
        np::Integer
    end

    function launch(manager::ArettoManager, params::Dict, launched::Array, c::Condition)
        try
            exehome  = params[:dir]
            exename  = params[:exename]
            exeflags = params[:exeflags]
            # cleanup old files
            map(rm, filter(t -> occursin(r"job.*\.out", t), readdir(exehome)))

            ## get the default arguments passed to addprocs()
            ## filter the non-default keys from stdkeys
            ## the filtered keys are the ones that are sent to Slurm
            ## TODO if jobname supplied in params, don't supply default name.
            ## TODO workerarg can be passed in exeflags
            stdkeys  = keys(Distributed.default_addprocs_params())
            p = filter(x->!(x[1] in stdkeys), params)
            srunargs = []
            for k in keys(p)
                if length(string(k)) == 1
                    push!(srunargs, "-$k")
                    val = p[k]
                    if length(val) > 0
                        push!(srunargs, "$(p[k])")
                    end
                else
                    k2 = replace(string(k), "_"=>"-")
                    val = p[k]
                    if length(val) > 0
                        push!(srunargs, "--$(k2)=$(p[k])")
                    else
                        push!(srunargs, "--$(k2)")
                    end
                end
            end

            np = manager.np
            jobname = "julia-$(getpid())" ## use default jobname..

            srun_cmd = `srun -J $jobname -n $np -o "job%4t.out" -D $exehome $(srunargs) $exename $exeflags $(worker_arg())`
            srun_proc = open(srun_cmd)

            @info "passed in params: $params"
            @info "default keys: $stdkeys"
            @info "filtered keys is: $p"
            @info "constructed srun args: $srunargs"
            @info "$srun_cmd"
            @info "srun_process: $(dump(srun_proc))"

            ## the `srun` Slurm command is now running.
            ## When the job allocates, it will create number $np text files called "job$4t.out" which have the ip address of the allocated resource where the worker Julia has laucnhed
            ## if a Julia worker is successfully launched, it writes julia_worker:9040#172.16.1.27 to STDOUT
            ## the STDOUT is rerouted to the jobfile by srun which is why -J is important
            for i = 0:(np - 1)
                print("connecting to worker $(i + 1) out of $np\r")
                local w=[]
                fn = "$exehome/job$(lpad(i, 4, "0")).out"
                t0 = time()
                ## wait for a bit (atleast a min).
                while true
                    if time() > t0 + 60 + np
                        warn("dropping worker: file not created in $(60 + np) seconds")
                        break
                    end
                    sleep(0.001)
                    ## if a valid output file is made by Slurm, open the file and read the ip address/hostname
                    if isfile(fn) && filesize(fn) > 0
                        w = open(fn) do f
                            return split(split(readline(f), ":")[2], "#")
                        end
                        break
                    end
                end
                ## if we get a valid hostname/ipaddress add a WorkerConfig object
                if length(w) > 0
                    config = WorkerConfig()
                    config.port = parse(Int, w[1])
                    config.host = strip(w[2])
                    # Keep a reference to the proc, so it's properly closed once
                    # the last worker exits.
                    config.userdata = srun_proc
                    push!(launched, config)
                    notify(c)
                end
            end
        catch e
            println("Error launching Slurm job.")
            rethrow(e)
        end
    end

    function manage(manager::ArettoManager, id::Integer, config::WorkerConfig, op::Symbol)
        # This function needs to exist, but so far we don't do anything
        # manage(manager::FooManager, id::Integer, config::WorkerConfig, op::Symbol) is called at different times during the worker's lifetime with appropriate op values:
        # with :register/:deregister when a worker is added / removed from the Julia worker pool.
        # with :interrupt when interrupt(workers) is called. The ClusterManager should signal the appropriate worker with an interrupt signal.
        # with :finalize for cleanup purposes.
    end
    addprocs_slurm(np::Integer; kwargs...) = addprocs(ArettoManager(np); kwargs...)
end

#using .Aretto
#addprocs(a, partition="defq", N=16)
#include("C:\\Users\\affan\\Documents\\abm-lab-cluster\\Aretto.jl")

max_stars_repo_name: affans/ArettoClusterABMLAB

max_issues_repo_name: affans/ArettoClusterABMLAB

max_forks_repo_name: affans/ArettoClusterABMLAB


text:
d  = 2+rand(5)
dl = rand(4)
du = rand(4)
M = Tridiagonal(dl, d, du)
F = lufact(M)
src = rand(5,5,5)
for dim = 1:3
    dest1 = mapslices(x->A_ldiv_B!(F, x), copy(src), dim)
    dest2 = similar(src)
    AxisAlgorithms.A_ldiv_B_md!(dest2, F, src, dim)
    @test_approx_eq dest1 dest2
end


text:
#----------------------------------------------------------------------
#* Model of a  distillation column containing:
#*	- NumberOfTrays like tray;
#*	- a kettle reboiler;
#*	- dynamic condenser;
#*	- a splitter which separate reflux and distillate;
#*	- a pump in reflux stream;
#*---------------------------------------------------------------------
type User_Distillation_kettle_cond
	User_Distillation_kettle_cond()=begin
		new(
			User_Section_ColumnBasic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour Outlet in the section",
				:PosX=>1,
				:PosY=>0.46,
				:Protected=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid Outlet in the section",
				:PosX=>1,
				:PosY=>0.58,
				:Protected=>true
			)),
			condenser(),
			reboiler(),
			splitter(),
			pump(),
			DanaReal(),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Reboiler",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Condenser",
				:Hidden=>true
			)),
			sourceNoFlow(Dict{Symbol,Any}(
				:Brief=>"No Inlet Flow to Reboiler",
				:Hidden=>true
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream From Top Condenser",
				:PosX=>0.67,
				:PosY=>0
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Vapour outlet stream From Top Condenser",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Top Splitter",
				:PosX=>1,
				:PosY=>0.33
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Top Splitter",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Reboiler",
				:PosX=>1,
				:PosY=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Reboiler",
				:Hidden=>true
			)),
			[
				:(CondenserUnity.InletV.F*trays(1).vV = alfaTopo *_base_1.Ah * sqrt(2*(trays(1).OutletV.P - CondenserUnity.OutletL.P + 1e-8 * "atm") / (_base_1.alfa*trays(1).rhoV))),
				:(CondenserUnity.InletV.F = 0 * "mol/s"),
				:(ConnectorCondenserVout.T = VapourDistillate.T),
				:(ConnectorCondenserVout.P = VapourDistillate.P),
				:(ConnectorCondenserVout.F = VapourDistillate.F),
				:(ConnectorCondenserVout.z = VapourDistillate.z),
				:(ConnectorSplitterOut.T = LiquidDistillate.T),
				:(ConnectorSplitterOut.P = LiquidDistillate.P),
				:(ConnectorSplitterOut.F = LiquidDistillate.F),
				:(ConnectorSplitterOut.z = LiquidDistillate.z),
				:(ConnectorReboilerLout.T = BottomProduct.T),
				:(ConnectorReboilerLout.P = BottomProduct.P),
				:(ConnectorReboilerLout.F = BottomProduct.F),
				:(ConnectorReboilerLout.z = BottomProduct.z),
				:(VapourDrawOff.F*_base_1.VapSideTrayIndex= trays.VapourSideStream.F),
				:(VapourDrawOff.T = trays(_base_1.VapourSideStreamLocation).VapourSideStream.T),
				:(VapourDrawOff.P = trays(_base_1.VapourSideStreamLocation).VapourSideStream.P),
				:(VapourDrawOff.z = trays(_base_1.VapourSideStreamLocation).VapourSideStream.z),
				:(LiquidDrawOff.F*_base_1.LiqSideTrayIndex= trays.LiquidSideStream.F),
				:(LiquidDrawOff.T = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.T),
				:(LiquidDrawOff.P = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.P),
				:(LiquidDrawOff.z = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.z),
				:(VapourDrawOffFlow = VapourDrawOff.F),
				:(LiquidDrawOffFlow = LiquidDrawOff.F),
			],
			[
				"","","","","","","","","","","","","","","","","","","","","","","","",
			],
			[:CondenserVapourFlow,],
			[:VapourDrawOff,:LiquidDrawOff,:CondenserUnity,:ReboilerUnity,:SplitterTop,:PumpUnity,:alfaTopo,:HeatToReboiler,:HeatToCondenser,:RebNoFlow,:VapourDistillate,:ConnectorCondenserVout,:LiquidDistillate,:ConnectorSplitterOut,:BottomProduct,:ConnectorReboilerLout,]
		)
	end
	_base_1::User_Section_ColumnBasic
	CondenserVapourFlow::DanaSwitcher
	VapourDrawOff::vapour_stream
	LiquidDrawOff::liquid_stream
	CondenserUnity::condenser
	ReboilerUnity::reboiler
	SplitterTop::splitter
	PumpUnity::pump
	alfaTopo::DanaReal
	HeatToReboiler::energy_stream
	HeatToCondenser::energy_stream
	RebNoFlow::sourceNoFlow
	VapourDistillate::vapour_stream
	ConnectorCondenserVout::stream
	LiquidDistillate::liquid_stream
	ConnectorSplitterOut::stream
	BottomProduct::liquid_stream
	ConnectorReboilerLout::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export User_Distillation_kettle_cond
function setEquationFlow(in::User_Distillation_kettle_cond)
	let switch=CondenserVapourFlow
		if CondenserVapourFlow==CondenserUnity.InletV.F < 1e-6 * "kmol/h"
			set(switch,"off")
		end
		if CondenserVapourFlow==trays(1).OutletV.P > CondenserUnity.OutletL.P + 1e-1 * "atm"
			set(switch,"on")
		end
		if switch=="on"
			addEquation(1)
		elseif switch=="off"
			addEquation(2)
		end
	end
	# Condenser Connector Equations
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	# Splitter Connector Equations
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	# Reboiler Connector Equations
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
end
function atributes(in::User_Distillation_kettle_cond,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/DistillationKettleCond"
	fields[:Brief]="Model of a distillation column with dynamic condenser and dynamic reboiler."
	fields[:Info]="== Specify ==
* the feed stream of each tray (Inlet);
* the Murphree eficiency for each tray Emv;
* the pump pressure difference;
* the heat supllied in reboiler and condenser;
* the condenser vapor outlet flow (OutletV.F);
* the reboiler liquid outlet flow (OutletL.F);
* both splitter outlet flows OR one of the splitter outlet flows and the splitter frac.

== Initial Conditions ==
* the trays temperature (OutletL.T);
* the trays liquid level (Level) OR the trays liquid flow (OutletL.F);
* (NoComps - 1) OutletL (OR OutletV) compositions for each tray;

* the condenser temperature (OutletL.T);
* the condenser liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions;

* the reboiler temperature (OutletL.T);
* the reboiler liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions.
"
	drive!(fields,_)
	return fields
end
User_Distillation_kettle_cond(_::Dict{Symbol,Any})=begin
	newModel=User_Distillation_kettle_cond()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# EXCLUDE FROM TESTING
using VectorEngine

function kernel(v::Int64)
    @veprintf("v = %ld\n", v)
    return nothing
end

# Compile the kernel
vefunc = VectorEngine.vefunction(kernel, Tuple{Int64})

# Launch kernel
vefunc(123)

synchronize()


text:
const CPetscOptions = Ptr{Cvoid}

#TODO: should it be <: AbstractDict{String,String}?
abstract type AbstractOptions{T}
end

"""
    GlobalOptions{T}()

The PETSc global options database.
"""
struct GlobalOptions{T} <: AbstractOptions{T}
end
Base.cconvert(::Type{CPetscOptions}, obj::GlobalOptions) = C_NULL

mutable struct Options{T} <: AbstractOptions{T}
    ptr::CPetscOptions
end
Base.cconvert(::Type{CPetscOptions}, obj::Options) = obj.ptr
Base.unsafe_convert(::Type{Ptr{CPetscOptions}}, obj::Options) =
    convert(Ptr{CPetscOptions}, pointer_from_objref(obj))

scalartype(::Options{T}) where {T} = T

@for_libpetsc begin
    function Options{$PetscScalar}()
        initialize($PetscScalar)
        opts = Options{$PetscScalar}(C_NULL)
        @chk ccall((:PetscOptionsCreate, $libpetsc), PetscErrorCode, (Ptr{CPetscOptions},), opts)
        finalizer(destroy, opts)
        return opts
    end
    function destroy(opts::Options{$PetscScalar})
        finalized($PetscScalar) ||
        @chk ccall((:PetscOptionsDestroy, $libpetsc), PetscErrorCode, (Ptr{CPetscOptions},), opts)
        return nothing
    end

    function Base.push!(::GlobalOptions{$PetscScalar}, opts::Options{$PetscScalar})
        @chk ccall((:PetscOptionsPush, $libpetsc), PetscErrorCode, (CPetscOptions,), opts)
        return nothing
    end
    function Base.pop!(::GlobalOptions{$PetscScalar})
        @chk ccall((:PetscOptionsPop, $libpetsc), PetscErrorCode, ())
        return nothing
    end
    function Base.setindex!(opts::AbstractOptions{$PetscScalar}, val, key)
        @chk ccall((:PetscOptionsSetValue, $libpetsc), PetscErrorCode,
            (CPetscOptions, Cstring, Cstring),
            opts, string('-',key), val == true ? C_NULL : string(val))
    end

    function view(opts::AbstractOptions{$PetscScalar}, viewer::Viewer{$PetscScalar}=ViewerStdout{$PetscScalar}(MPI.COMM_SELF))
        @chk ccall((:PetscOptionsView, $libpetsc), PetscErrorCode,
                  (CPetscOptions, CPetscViewer),
                  opts, viewer);
        return nothing
    end
end

"""
    Options{T}(kw => arg, ...)

Create a new PETSc options database.
"""
function Options{T}(ps::Pair...) where {T}
    opts = Options{T}()
    for (k,v) in ps
        opts[k] = v
    end
    return opts
end

Base.show(io::IO, opts::AbstractOptions) = _show(io, opts)

"""
    with(f, opts::Options)

Call `f()` with the [`Options`](@ref) `opts` set temporarily (in addition to any global options).
"""
function with(f, opts::Options{T}) where {T}
  global_opts = GlobalOptions{T}()
  push!(global_opts, opts)
  try
    f()
  finally
    pop!(global_opts)
  end
end


text:
_getlocation(d::Distributions.LogNormal) = d.μ
_getlocation(d::Distributions.MvLogNormal) = Distributions.location(d)
_getscale(d::Distributions.Normal) = var(d)
_getscale(d::Distributions.MvNormal) = cov(d)
_getscale(d::Distributions.LogNormal) = d.σ^2
_getscale(d::Distributions.MvLogNormal) = Distributions.scale(d)

#testing of the normal distribution
for p in [(1.0,2.0,4.0,2.0),
          (ones(2),4.0*eye(2),4.0*eye(2),2.0*ones(2)),
          (ones(2),2.0*ones(2),4.0*eye(2),2.0*ones(2))]
    cp2 = copy(p[2])
    n = distribution(:normal,p[1],p[2])
    @test_approx_eq mean(n) p[1]
    n = recenter(n,p[4])
    @test_approx_eq mean(n) p[4]
    @test_approx_eq _getscale(n) p[3]
    n = rescale(n,2.0)
    @test_approx_eq _getscale(n) 2.0*2.0*p[3]
    n = update(n,p[1],cp2)
    @test_approx_eq mean(n) p[1]
    @test_approx_eq _getscale(n) p[3]
end

for p in [(eye(2),4.0*eye(2)),
          (ones(2),4.0*eye(2))]
    n = distribution(:normal,eye(2))
    @test mean(n) == zeros(2) && cov(n) == eye(2)
    @test_throws MethodError recenter(n,ones(2))
    n = rescale(n,2.0)
    @test_approx_eq cov(n) 4.0*eye(2)
    n = update(n,p[1])
    @test_approx_eq cov(n) eye(2)
end

#testing of the lognormal distribution
for p in [(1.0,2.0,4.0,2.0),
          (ones(2),4.0*eye(2),4.0*eye(2),2.0*ones(2)),
          (ones(2),2.0*ones(2),4.0*eye(2),2.0*ones(2))]
    cp2 = copy(p[2])
    l1 = distribution(:lognormal,p[1],p[2])
    @test_approx_eq _getlocation(l1) log(p[1])
    l1 = recenter(l1,p[4])
    @test_approx_eq _getlocation(l1) log(p[4])
    @test_approx_eq _getscale(l1) p[3]
    l1 = rescale(l1,2.0)
    @test_approx_eq _getscale(l1) 4.0*p[3]
    l1 = update(l1,p[1],cp2)
    @test_approx_eq _getlocation(l1) log(p[1])
    @test_approx_eq _getscale(l1) p[3]
end

#testing various uniform distributions
for s in [:uniform,:laplace,:triangular]
    l1 = distribution(s,1.0,0.5)
    @test mean(l1) == 1.0 && scale(l1) == 0.5
    l1 = recenter(l1,2.0)
    @test mean(l1) == 2.0 && scale(l1) == 0.5
    l1 = rescale(l1,3.0)
    @test mean(l1) == 2.0 && scale(l1) == 1.5
    l1 = update(l1,1.0,1.0)
    @test mean(l1) == 1.0 && scale(l1) == 1.0
end

#testing the bactrian distribution
for p in [(:normal,1.0),
          (:laplace,sqrt(2.0)),
          (:triangular,6.0)]
    bs1 = GeneralizedMetropolisHastings._bactrianscale(Val{p[1]},0.5,0.95)
    bs2 = GeneralizedMetropolisHastings._bactrianscale(Val{p[1]},1.0,0.95)
    @test_approx_eq bs1 0.5*sqrt(p[2]*(1-0.95*0.95))
    @test_approx_eq bs2 1.0*sqrt(p[2]*(1-0.95*0.95))
    m1 = GeneralizedMetropolisHastings._mixturemodel(p[1],2.0,0.5,0.95)
    @test m1.components[1] == distribution(p[1],2.0-0.5*0.95,bs1)
    @test m1.components[2] == distribution(p[1],2.0+0.5*0.95,bs1)
    b1 = distribution(:bactrian,2.0,0.5,p[1],0.95)
    @test mean(b1.mixture) == 2.0 && b1.μ == 2.0 && b1.σ == 0.5 && b1.m == 0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[1] 2.0-0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[2])[1] 2.0+0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[2] bs1
    @test_approx_eq Distributions.params(b1.mixture.components[2])[2] bs1
    b1 = recenter(b1,1.0)
    @test mean(b1.mixture) == 1.0 && b1.μ == 1.0 && b1.σ == 0.5 && b1.m == 0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[1] 1.0-0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[2])[1] 1.0+0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[2] bs1
    @test_approx_eq Distributions.params(b1.mixture.components[2])[2] bs1
    b1 = rescale(b1,2.0)
    @test mean(b1.mixture) == 1.0 && b1.μ == 1.0 && b1.σ == 1.0 && b1.m == 0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[1] 1.0-1.0*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[2])[1] 1.0+1.0*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[2] bs2
    @test_approx_eq Distributions.params(b1.mixture.components[2])[2] bs2
    b1 = update(b1,2.0,0.5)
    @test mean(b1.mixture) == 2.0 && b1.μ == 2.0 && b1.σ == 0.5 && b1.m == 0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[1] 2.0-0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[2])[1] 2.0+0.5*0.95
    @test_approx_eq Distributions.params(b1.mixture.components[1])[2] bs1
    @test_approx_eq Distributions.params(b1.mixture.components[2])[2] bs1
    @test_approx_eq Distributions.logpdf(b1,1.0) log(Distributions.pdf(b1.mixture,1.0))
    @test mean(b1) == mean(b1.mixture)
    @test var(b1) == var(b1.mixture)
    @test (srand(324) ; Distributions.rand(b1)) == (srand(324) ; Distributions.rand(b1.mixture))
    @test (srand(324) ; Distributions.rand!(b1,zeros(2))) == (srand(324) ; Distributions.rand!(b1.mixture,zeros(2)))
end

for p in [(:normal,Distributions.Normal),
          (:lognormal,Distributions.LogNormal),
          (:uniform,Distributions.Uniform),
          (:laplace,Distributions.Laplace),
          (:triangular,Distributions.SymTriangularDist)]
    d = distributions(p[1],[0.0,1.0],[1.0,2.0])
    @test eltype(d) == p[2]
end

#a test to ensure that all multivariate distributions copy their matrix arguments upon creation
#so that creating two from the same matrix does not introduce side-effects when recentering or rescaling
for p in [(:normal,ones(2),eye(2),2ones(2),mean,2ones(2),ones(2),cov,0.01*eye(2),eye(2)),
          (:lognormal,ones(2),eye(2),2ones(2),Distributions.location,log(2ones(2)),zeros(2),Distributions.scale,0.01*eye(2),eye(2))]
    println("Testing multivariate ",p[1]," copying behaviour")
    d1 = distribution(p[1],p[2],p[3])
    d2 = distribution(p[1],p[2],p[3])
    d1 = recenter(d1,p[4])
    @test_approx_eq p[5](d1) p[6]
    @test_approx_eq p[5](d2) p[7]
    d1 = rescale(d1,0.1)
    @test_approx_eq p[8](d1) p[9]
    @test_approx_eq p[8](d2) p[10]
end

for p in [(:normal,eye(2),cov,0.01eye(2),eye(2))]
    println("Testing multivariate ",p[1]," (zero mean) copying behaviour")
    d1 = distribution(p[1],p[2])
    d2 = distribution(p[1],p[2])
    d1 = rescale(d1,0.1)
    @test_approx_eq p[3](d1) p[4]
    @test_approx_eq p[3](d2) p[5]
end

#test of vectorized distribution factory function
d = distributions(:bactrian,[0.0,1.0],[1.0,2.0],:normal,0.95)
@test eltype(d) == GeneralizedMetropolisHastings.Bactrian{Float64}


text:
using ThreeBodyProblem
using DifferentialEquations
using Plots
using LinearAlgebra
# using ForwardDiff

# The first step is to define the system. ThreeBodyProblem.jl has built-in functions for common systems like Sun/Jupiter.
sys = sun_jupiter()

## There's even a recipe to plot the system for us!
plot(sys)

## We can see the Lagrange Points, but where are the Sun and Jupiter? Turns out that the distance between them is so much larger than the radius of either one that it makes it hard to see them on a truly scaled picture. Let's scale the Sun and Jupiter to make them visible.
plot(sys, scaled = true)

## That's better! Now let's compute the L1 and L2 Lagrange points of our system

L1, L2 = computeLpts(sys)
# Declare state vectors for L1 and L2 (with zero velocity)
rv1 = [L1; zeros(3)] # state at Lagrange point L1 (zero velocity added on)
rv2 = [L2; zeros(3)] # state at Lagrange point L2

## The Lagrange points are equilibrium points in our dynamics. This means that an object placed there perfectly will stay there forever. But, Lagrange points L1, L2, and L3 are unstable, meaning if you perturb the object slightly, it may fall away. This is how we will generate the invariant manifolds.
## We first need to linearize the dynamics about these Lagrange points to determine which directions are stable and unstable

Φ₀ = I(6) # The initial State Transition Matrix (STM)
w₀ = vcat(rv1,reshape(Φ₀,36,1)) # Augmented state with state and STM
# Set up the ODE problem
tspan = (0.,1.)
prob = ODEProblem(CR3BPstm!,w₀,tspan,sys)
sol = solve(prob, reltol=1e-6)
Φₜ = Matrix(reshape(sol.u[end][7:42],6,6))

# last time step should match very closely with rv1
rvₜ = sol.u[end][1:6]
err = norm(rv1 - rvₜ)

# Compute the eigenvectors of the Monodromy Matrix
Λ,V = eigen(Φₜ,sortby=isreal) # Λ is vector of eigenvalues, V is matrix of eigenvectors
Yw = real(V[:,findall(isreal, Λ)]) # Eigenvectors corresponding to real eigenvalues
Λ = Λ[findall(isreal, Λ)] # Purely real eigenvalues (have 0.0 imaginary component)
Yws = Yw[:,findmin(real(Λ))[2]] # Eigenvector associated with stable eigenvalue λ < 1
Ywu = Yw[:,findmax(real(Λ))[2]] # Eigenvector associated with unstable eigenvalue λ > 1

## We perturb our states in the stable and unstable directions
α = 1e-6

# L1 points
rv1up = rv1 + α*Ywu # unstable manifold + side
rv1un = rv1 - α*Ywu # unstable manifold - side
rv1sp = rv1 + α*Yws # stable manifold + side
rv1sn = rv1 - α*Yws # stable manifold - side

# Integrate trajectories for 2 Jupiter periods
tf = 2*2π
tspan_forward = (0.,tf) # make sure to add decimal points so tspan contains floating point values, not integers
tspan_back = (0.,-tf)
myreltol = 1e-12

# Set up the ODE Problems
prob1up = ODEProblem(CR3BPdynamics!,rv1up,tspan_forward,sys) # Unstable positive
prob1un = ODEProblem(CR3BPdynamics!,rv1un,tspan_forward,sys) # Unstable negative
prob1sp = ODEProblem(CR3BPdynamics!,rv1sp,tspan_back,sys) # Stable positive
prob1sn = ODEProblem(CR3BPdynamics!,rv1sn,tspan_back,sys) # Stable negative

# Solutions to the ODEs
sol1up = solve(prob1up, reltol=myreltol)
sol1un = solve(prob1un, reltol=myreltol)
sol1sp = solve(prob1sp, reltol=myreltol)
sol1sn = solve(prob1sn, reltol=myreltol)

plot(sys,scaled=true)
plot!(sol1up,vars=(1,2),label="Wu+",linecolor=:red)
plot!(sol1un,vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(sol1sp,vars=(1,2),label="Ws+",linecolor=:blue)
plot!(sol1sn,vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,ylims=[-1,1],xlims=[-0.6,1.3],legend=:topleft,flip=false)

## We can see two strands moving in towards the sun, but the rest is a jumbled mess. Let's zoom in to see what's going on. Also, I need to make sure things are scaled correctly, or Jupiter will look way too big.
plot(sys)
plot!(sol1up,vars=(1,2),label="Wu+",linecolor=:red)
plot!(sol1un,vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(sol1sp,vars=(1,2),label="Ws+",linecolor=:blue)
plot!(sol1sn,vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,ylims=[-0.04,0.04],xlims=[0.93,1.06],legend=:outerright,flip=false)

# L2 points
rv2up = rv2 + α*Ywu # unstable manifold + side
rv2un = rv2 - α*Ywu # unstable manifold - side
rv2sp = rv2 + α*Yws # stable manifold + side
rv2sn = rv2 - α*Yws # stable manifold - side

## Now we look at the L2 point manifolds
prob2up = ODEProblem(CR3BPdynamics!,rv2up,tspan_forward,sys) # Unstable positive
prob2un = ODEProblem(CR3BPdynamics!,rv2un,tspan_forward,sys) # Unstable negative
prob2sp = ODEProblem(CR3BPdynamics!,rv2sp,tspan_back,sys) # Stable positive
prob2sn = ODEProblem(CR3BPdynamics!,rv2sn,tspan_back,sys) # Stable negative

sol2up = solve(prob2up, reltol=1e-6)
sol2un = solve(prob2un, reltol=1e-6)
sol2sp = solve(prob2sp, reltol=1e-6)
sol2sn = solve(prob2sn, reltol=1e-6)

plot(sys, scaled=true)
plot!(sol2up,vars=(1,2),label="Wu+",linecolor=:red)
plot!(sol2un,vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(sol2sp,vars=(1,2),label="Ws+",linecolor=:blue)
plot!(sol2sn,vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,ylims=[-1.0,1.0],xlims=[-0.3,1.5],legend=:topright,flip=false)
## This time there are two strands moving out away from the Sun, along with another jumble around Jupiter. Let's zoom in again to take a look.
plot(sys)
plot!(sol2up,vars=(1,2),label="Wu+",linecolor=:red)
plot!(sol2un,vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(sol2sp,vars=(1,2),label="Ws+",linecolor=:blue)
plot!(sol2sn,vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,ylims=[-0.04,0.04],xlims=[0.93,1.07],legend=:outerright,flip=false)


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

module AirSpaceImpl

export
    initialize,
    update,

    updateObjectState,
    getAllObjectStates,
    getAll,
    updateAll,

    AirSpace,
    ASWMState

using AbstractWorldModelImpl
using AbstractWorldModelInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractWorldModelInterfaces.updateObjectState
import AbstractWorldModelInterfaces.getAllObjectStates
import AbstractWorldModelInterfaces.getAll
import AbstractWorldModelInterfaces.updateAll

type ASWMState

    t::Float64
    x::Float64
    y::Float64
    h::Float64
    vx::Float64
    vy::Float64
    vh::Float64
end

type AirSpace <: AbstractWorldModel

    number_of_aircraft::Int

    states::Vector{ASWMState}
    states_update::Vector{ASWMState}

    function AirSpace(number_of_aircraft::Int)

        obj = new()

        obj.number_of_aircraft = number_of_aircraft

        obj.states = Array(ASWMState, number_of_aircraft)
        obj.states_update = Array(ASWMState, number_of_aircraft)

        return obj
    end
end

function initialize(as::AirSpace, aircraft_number::Int, state)
  initialize(as,aircraft_number,convert(ASWMState, state))
end

function initialize(as::AirSpace, aircraft_number::Int, state::ASWMState)

    as.states[aircraft_number] = state
end

function updateObjectState(as::AirSpace, aircraft_number::Int, state::ASWMState)

    as.states_update[aircraft_number] = state
end

update(as::AirSpace, aircraft_number::Int, state) = update(as, aircraft_number, convert(ASWMState,state))

update(as::AirSpace, aircraft_number::Int, state::ASWMState) = updateObjectState(as, aircraft_number, state)

function getAllObjectStates(as::AirSpace)

    return as.states
end

getAll(as::AirSpace) = getAllObjectStates(as)

function updateAll(as::AirSpace)

    as.states = copy(as.states_update)
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
"""Gaussian filtering and smoothing"""

"""
    predict!(x_out, x_curr, Ah, Qh)

PREDICT step in Kalman filtering for linear dynamics models.
In-place implementation of [`predict`](@ref), saving the result in `x_out`.

```math
m_{n+1}^P = A(h)*m_n
P_{n+1}^P = A(h)*P_n*A(h) + Q(h)
```

See also: [`predict`](@ref)
"""
function predict!(x_out::Gaussian, x_curr::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)
    predict_mean!(x_out, x_curr, Ah, Qh)
    predict_cov!(x_out, x_curr, Ah, Qh)
    return x_out
end
function predict_mean!(x_out::Gaussian, x_curr::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)
    mul!(x_out.μ, Ah, x_curr.μ)
    return x_out.μ
end
function predict_cov!(x_out::Gaussian, x_curr::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)
    out_cov = X_A_Xt(x_curr.Σ, Ah) + Qh
    copy!(x_out.Σ, out_cov)
    return x_out.Σ
end

# PSDMatrix Version of this - But, before using QR try it with Cholesky!
function predict_cov!(x_out::PSDGaussian, x_curr::PSDGaussian, Ah::AbstractMatrix, Qh::PSDMatrix)
    _L = [Ah*x_curr.Σ.L Qh.L]
    out_cov = Symmetric(_L*_L')
    chol = cholesky!(out_cov, check=false)

    if issuccess(chol)
        PpL = chol.L
        copy!(x_out.Σ, PSDMatrix(PpL))
        return x_out.Σ
    else
        _, R = qr(_L')
        out_cov = PSDMatrix(LowerTriangular(collect(R')))
        copy!(x_out.Σ, out_cov)
        return x_out.Σ
    end
end

"""
    predict(x_curr::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)

See also: [`predict!`](@ref)
"""
function predict(x_curr::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)
    x_out = copy(x_curr)
    predict!(x_out, x_curr, Ah, Qh)
    return x_out
end

"""
    update!(x_out, x_pred, measurement, H, R=0)

UPDATE step in Kalman filtering for linear dynamics models, given a measurement `Z=N(z, S)`.
In-place implementation of [`update`](@ref), saving the result in `x_out`.

```math
K = P_{n+1}^P * H^T * S^{-1}
m_{n+1} = m_{n+1}^P + K * (0 - z)
P_{n+1} = P_{n+1}^P - K*S*K^T
```

Implemented in Joseph Form.

See also: [`predict`](@ref)
"""
function update!(x_out::Gaussian, x_pred::Gaussian, measurement::Gaussian,
                 H::AbstractMatrix, R=0)
    @assert iszero(R)
    z, S = measurement.μ, measurement.Σ
    m_p, P_p = x_pred.μ, x_pred.Σ

    S_inv = inv(S)
    K = P_p * H' * S_inv

    x_out.μ .= m_p .+ K * (0 .- z)
    copy!(x_out.Σ, X_A_Xt(P_p, (I-K*H)))
    return x_out
end
"""
    update(x_pred, measurement, H, R=0)

See also: [`update!`](@ref)
"""
function update(x_pred::Gaussian, measurement::Gaussian, H::AbstractMatrix, R=0)
    @assert iszero(R)
    x_out = copy(x_pred)
    update!(x_out, x_pred, measurement, H, R)
    return x_out
end

"""
    smooth(x_curr::Gaussian, x_next_smoothed::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)

SMOOTH step of the (extended) Kalman smoother, or (extended) Rauch-Tung-Striebel smoother.
It is implemented in Joseph Form:
```math
m_{n+1}^P = A(h)*m_n
P_{n+1}^P = A(h)*P_n*A(h) + Q(h)

G = P_n * A(h)^T * (P_{n+1}^P)^{-1}
m_n^S = m_n + G * (m_{n+1}^S - m_{n+1}^P)
P_n^S = (I - G*A(h)) P_n (I - G*A(h))^T + G * Q(h) * G + G * P_{n+1}^S * G
```
"""
function smooth(x_curr::Gaussian, x_next_smoothed::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)
    x_pred = predict(x_curr, Ah, Qh)

    P_p = x_pred.Σ
    P_p_inv = inv(P_p)

    G = x_curr.Σ * Ah' * P_p_inv

    smoothed_mean = x_curr.μ + G * (x_next_smoothed.μ - x_pred.μ)

    smoothed_cov = (
        X_A_Xt(x_curr.Σ, (I - G*Ah))
        + X_A_Xt(Qh, G)
        + X_A_Xt(x_next_smoothed.Σ, G)
    )

    assert_nonnegative_diagonal(smoothed_cov)

    x_curr_smoothed = Gaussian(smoothed_mean, smoothed_cov)
    return x_curr_smoothed, G
end


text:
using OrdinaryDiffEq, DiffEqDevTools, DiffEqBase, Test,
      ODEInterface, ODEInterfaceDiffEq

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_bigfloatlinear,
                                               prob_ode_linear,
                                               prob_ode_2Dlinear,
                                               prob_ode_bigfloat2Dlinear

probbig = prob_ode_bigfloat2Dlinear
probnum = prob_ode_linear
probnumbig = prob_ode_bigfloatlinear
prob = prob_ode_2Dlinear

dts = (1/2) .^ (7:-1:4)
testTol = .2
bools = Vector{Bool}(undef,0)

## DP5()
sim = test_convergence(dts,probnum,DP5())
@test abs.(sim.𝒪est[:l2]-5) < testTol
sim = test_convergence(dts,prob,DP5())
@test abs.(sim.𝒪est[:l2]-5) < testTol

tabalg = ExplicitRK()
sol1 =solve(probnum,DP5(),dt=1/2^6,adaptive=false,save_everystep=false)
sol2 =solve(probnum,tabalg,dt=1/2^6,adaptive=false,save_everystep=false)

@test sol1.u[end] - sol2.u[end] < 1e-10

sol1 =solve(prob,DP5(),dt=1/2^6,adaptive=false,save_everystep=false)
sol2 =solve(prob,tabalg,dt=1/2^6,adaptive=false,save_everystep=false)

@test minimum(sol1.u[end] - sol2.u[end] .< 3e-10)

sol1 =solve(probnum,DP5(),dt=1/2^6,beta2=0.04)
sol2 =solve(probnum,tabalg,dt=1/2^6,beta2=0.04)

# Should be identical
sol1 =solve(probnum,DP5())
sol2 =solve(probnum,tabalg,beta2=0.04,beta1=0.17)
sol3 =solve(probnum,dopri5())

@test sol1.t ≈ sol2.t
@test sol1.t ≈ sol3.t atol=1e-6

sol1 =solve(prob,DP5(),dt=1/8)
sol2 =solve(prob,tabalg,beta2=0.04,beta1=0.17,dt=1/8)
sol3 =solve(prob,dopri5(),dt=1/8)

@test sol1.t ≈ sol2.t
@test sol1.t ≈ sol3.t atol=1e-5

sol4 =solve(prob,DP5(),dt=1/8,calck=false)

@test sol1.t == sol4.t

### DP8()

dts = (1/2) .^ (3:-1:1)
sim = test_convergence(dts,probnumbig,DP8())
@test abs.(sim.𝒪est[:l2]-8) < testTol
sim = test_convergence(dts,probbig,DP8())
@test abs.(sim.𝒪est[:l2]-8) < testTol

sol1 =solve(probnum,DP8(),dt=1/2^6,adaptive=false,save_everystep=false)
sol2 =solve(probnum,DP8(),dt=1/2^6)

# Should be identical
sol1 =solve(probnum,DP8())
sol2 =solve(probnum,dop853())

@test sol1.t ≈ sol2.t

# Should be identical
sol1 =solve(probbig,DP8(),dt=1/2^6)
sol2 =solve(probbig,dop853(),dt=1/2^6)

@test sol1.t ≈ sol2.t


text:
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Test
using DataDrivenDiffEq.Optimize

include("./basis.jl")
include("./koopman.jl")
include("./sindy.jl")
include("./isindy.jl")
include("./utils.jl")
include("./optimize.jl")


text:
function get_component_to_pm(ix::Int, storage::S) where {S <: PSY.Storage}
    qlims = PSY.get_reactive_power_limits(storage)
    PM_storage = Dict{String, Any}(
        "x" => 0.0,
        "r" => 0.0,
        "ps" => 0.0,
        "qs" => 0.0,
        "p_loss" => 0.0,
        "q_loss" => 0.0,
        "energy" => PSY.get_initial_energy(storage),
        "energy_rating" => PSY.get_state_of_charge_limits(storage).max,
        "thermal_rating" => PSY.get_rating(storage),
        "status" => PSY.get_available(storage),
        "source_id" => Any[PSY.get_name(storage), ix],
        "discharge_rating" => PSY.get_output_active_power_limits(storage).max,
        "storage_bus" => PSY.get_number(PSY.get_bus(storage)),
        "name" => PSY.get_name(storage),
        "charge_efficiency" => PSY.get_efficiency(storage).in,
        "index" => ix,
        "qmax" => isnothing(qlims) ? 0.0 : qlims.max,
        "qmin" => isnothing(qlims) ? 0.0 : qlims.min,
        "charge_rating" => PSY.get_input_active_power_limits(storage).max,
        "discharge_efficiency" => PSY.get_efficiency(storage).out,
    )
    return PM_storage
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
module OffPolicyDecisionMaking

using DecisionMakingPolicies

using Random, Statistics, Distributions
using Zygote
using DecisionMakingEnvironments

export BanditExperience, Trajectory
export length, push!

export ope
export estimate_returns, estimate_returns!, estimate_returns_withentropy, estimate_returns_withentropy!

export AbstractImportanceSampling, UnweightedIS, WeightedIS, IS, PDIS, WIS, WPDIS

export andersons_ci, tdist_ci

export AbstractSeldonianProblem, SeldonianProblem
export solve

export AbstractSplitMethod, SplitLastK, SplitLastKKeepTest
export collect_and_split!
export HICOPI

include("importance_sampling/importancesampling.jl")
include("confidenceintervals.jl")
include("ope.jl")
include("seldonian.jl")
include("highconfidence.jl")

end


text:
#####
##### make sure there are no temporary files files (e.g., *.vtu) left around from the build in `GENERATED_DIR`
#####

file_extensions_to_remove = [".vtu", ".pvtu", ".csv", ".vtk", ".dat", ".nc"]
generated_files = [
    joinpath(root, f)
    for (root, dirs, files) in Base.Filesystem.walkdir(GENERATED_DIR)
    for f in files
]
println("Generated files: $(generated_files)")

filter!(
    x -> any([endswith(x, y) for y in file_extensions_to_remove]),
    generated_files,
)

println("Deleting files: $(generated_files)")
for f in generated_files
    rm(f)
end


text:
export latex

@api latex => LaTeX <: Tile begin
    arg(source::String)
    kwarg(block::Bool=false)
end

render(l::LaTeX) =
    Elem("ka-tex", source=l.source, block=l.block)


text:
module Test_sylvkr

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing Sylvester equation solvers based on Kronecker expansions" begin

n = 10; m = 7;
ar = rand(n,n)
br = rand(m,m)
cr = rand(n,m)
dr = rand(n,n)
er = rand(m,m)
fr = rand(n,m)
ac = ar+im*rand(n,n)
bc = br+im*rand(m,m)
cc = cr+im*rand(n,m)
dc = dr+im*rand(n,n)
ec = er+im*rand(m,m)
fc = fr+im*rand(n,m)
qr = cr*cr'
qc = cc*cc'
reltol = sqrt(eps(1.))

# solving Sylvester equations
@time x = sylvckr(ar,br,cr)
@test norm(ar*x+x*br-cr)/norm(x) < reltol

@time x = sylvckr(ac,bc,cc)
@test norm(ac*x+x*bc-cc)/norm(x) < reltol

@time x = sylvdkr(ar,br,cr)
@test norm(ar*x*br+x-cr)/norm(x) < reltol

@time x = sylvdkr(ac,bc,cc)
@test norm(ac*x*bc+x-cc)/norm(x) < reltol

# solving Lyapunov equations
@time x = sylvckr(ar,ar',-qr)
@test norm(ar*x+x*ar'+qr)/norm(x) < reltol

@time x = sylvckr(ac',ac,-qc)
@test norm(ac'*x+x*ac+qc)/norm(x) < reltol

@time x = sylvdkr(-ar,ar',qr)
@test norm(ar*x*ar'-x+qr)/norm(x) < reltol

@time x = sylvdkr(-ac',ac,qc)
@test norm(ac'*x*ac-x+qc)/norm(x) < reltol

# solving generalized Sylvester equations
@time x = gsylvkr(ar,br,dr,er,cr)
@test norm(ar*x*br+dr*x*er-cr)/norm(x) < reltol

@time x = gsylvkr(ac,bc,dc,ec,cc)
@test norm(ac*x*bc+dc*x*ec-cc)/norm(x) < reltol

@time x = gsylvkr(ar,br',-dr,er',cr)
@test norm(ar*x*br'-dr*x*er'-cr)/norm(x) < reltol

@time x = gsylvkr(ac',bc,-dc,ec',cc)
@test norm(ac'*x*bc-dc*x*ec'-cc)/norm(x) < reltol

# solving generalized Lyapunov equations
@time x = gsylvkr(ar,dr',dr,ar',-qr)
@test norm(ar*x*dr'+dr*x*ar'+qr)/norm(x) < reltol

@time x = gsylvkr(ac',dc,dc',ac,-qc)
@test norm(ac'*x*dc+dc'*x*ac+qc)/norm(x) < reltol

@time x = gsylvkr(ar,ar',-dr,dr',-qr)
@test norm(ar*x*ar'-dr*x*dr'+qr)/norm(x) < reltol

@time x = gsylvkr(ac',ac,-dc',dc,-qc)
@test norm(ac'*x*ac-dc'*x*dc+qc)/norm(x) < reltol

# solving Sylvester systems
@time x, y = sylvsyskr(ar,br,cr,dr,er,fr)
@test norm(ar*x+y*br-cr)/max(norm(x),norm(y)) < reltol &&
      norm(dr*x+y*er-fr)/max(norm(x),norm(y)) < reltol

@time x, y = sylvsyskr(ar,-br,cr,dr,-er,fr)
@test norm(ar*x-y*br-cr)/max(norm(x),norm(y)) < reltol &&
      norm(dr*x-y*er-fr)/max(norm(x),norm(y)) < reltol

@time x, y = sylvsyskr(ac,bc,cc,dc,ec,fc)
@test norm(ac*x+y*bc-cc)/max(norm(x),norm(y)) < reltol &&
      norm(dc*x+y*ec-fc)/max(norm(x),norm(y)) < reltol

@time x, y = sylvsyskr(ac,-bc,cc,dc,-ec,fc)
@test norm(ac*x-y*bc-cc)/max(norm(x),norm(y)) < reltol &&
      norm(dc*x-y*ec-fc)/max(norm(x),norm(y)) < reltol

@time x, y = dsylvsyskr(ar',br',cr,-dr',er',-fr)
@test norm(ar'*x-dr'*y-cr)/max(norm(x),norm(y)) < reltol &&
      norm(x*br'+y*er'+fr)/max(norm(x),norm(y)) < reltol

end

end


text:
using CompScienceMeshes, BEAST
using Test

fn = joinpath(dirname(@__FILE__),"assets","sphere35.in")
m = readmesh(fn)
t = Maxwell3D.singlelayer(wavenumber=1.0)
X = raviartthomas(m)
numfunctions(X)

##
X1 = subset(X,1:1)
numfunctions(X1)

T1 = assemble(t,X1,X)
T2 = BEAST.assemblerow(t,X1,X)
#
# T3 = assemble(t,X,X1)
# T4 = BEAST.assemblecol(t,X,X1)
#
# @test T1 == T2

# T2 = BEAST.assembleblock(t,X,X)
@test T1 == T2

I = [3,2,7]
X1 = subset(X,I)
T2 = zeros(scalartype(t,X1,X1),numfunctions(X1),numfunctions(X1))
store(v,m,n) = (T2[m,n] += v)

test_elements, test_assembly_data,
	trial_elements, trial_assembly_data,
	quadrature_data, zlocal = BEAST.assembleblock_primer(t,X,X)

BEAST.assembleblock_body!(t,
	X, I, test_elements,  test_assembly_data,
    X, I, trial_elements, trial_assembly_data,
    quadrature_data, zlocal, store)

T1 = assemble(t,X1,X1)
@test T1 == T2

# @time BEAST.assembleblock_body!(t,
# 	X, I, test_elements,  test_assembly_data,
#     X, I, trial_elements, trial_assembly_data,
#     quadrature_data, zlocal, store)
# @time assemble(t,X1,X1)

T3 = zeros(scalartype(t,X1,X1),numfunctions(X1),numfunctions(X1))
store3(v,m,n) = (T3[m,n] += v)

blkasm = BEAST.blockassembler(t,X,X)
blkasm(I,I,store3)

T4 = assemble(t,X,X)
@test T3 == T4[I,I]

# @time blkasm(I,I)
# @time assemble(t,X1,X1)


text:
@doc raw"""
    gaussian([T=Float64], point; x, y, fwhm, amp=1, theta=0, bkg=0)
    gaussian([T=Float64], px, py; x, y, fwhm, amp=1, theta=0, bkg=0)

An unnormalized bivariate Gaussian distribution. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. If `theta` is given, the PSF will be rotated by `theta` degrees counter-clockwise from the x-axis. If `bkg` is given it will be added as a scalar to the PSF.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that `theta` has no effect for isotropic distributions and is degenerate with the `fwhm` parameters (i.e., theta=90 is the same as reversing the `fwhm` tuple)

# Functional form
```
f(x | x̂, FWHM) = exp[-4ln(2) * ||x - x̂|| / FWHM^2]
```
where `x̂` and `x` are position vectors (indices) `||⋅||` represents the square-distance, and `FWHM` is the full width at half-maximum. If `FWHM` is a scalar, the Gaussian distribution will be isotropic. If `FWHM` is a vector or tuple, the weighting is applied along each axis (diagonal).
"""
gaussian(T, px, py; x, y, fwhm, amp=one(T), theta=0, bkg=0) = convert(T, _gaussian(px, py, x, y, fwhm, amp, theta, bkg))

"""
    normal

An alias for [`gaussian`](@ref)
"""
const normal = gaussian

# this is the factor to convert 1/(2σ²) to 1/(2fwhm²)
const GAUSS_PRE = -4 * log(2)

# isotropic
function _gaussian(px, py, x, y, fwhm, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    !iszero(theta) && @warn "isotropic gaussian is not affected by non-zero rotation angle $theta"
    # unnormalized gaussian likelihood
    sqmahab = dx^2 + dy^2
    sqmahab /= fwhm^2
    return amp * exp(GAUSS_PRE * sqmahab) + background
end

# bivariate
function _gaussian(px, py, x, y, fwhm::BivariateLike, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    if !iszero(theta)
        dx, dy = rotate_point(dx, dy, theta)
    end
    # unnormalized gaussian likelihood
    fwhmx, fwhmy = fwhm
    sqmahab = (dx / fwhmx)^2 + (dy / fwhmy)^2
    return amp * exp(GAUSS_PRE * sqmahab) + background
end


text:
using FRGRealTime
using Documenter

DocMeta.setdocmeta!(FRGRealTime, :DocTestSetup, :(using FRGRealTime); recursive=true)

makedocs(;
    modules=[FRGRealTime],
    authors="Yang-yang Tan, Yong-rui Chen, Wei-jie Fu",
    repo="https://github.com/Yangyang-Tan/FRGRealTime.jl/blob/{commit}{path}#{line}",
    sitename="FRGRealTime.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Yangyang-Tan.github.io/FRGRealTime.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Example" => "markdown/Example.md",
        "Numerical Tricks" => "markdown/NumericalTricks.md",
        "Check" => "markdown/Check.md"
    ],
)

deploydocs(;
    repo="github.com/Yangyang-Tan/FRGRealTime.jl",
)


text:
using SnoopCompile

SnoopCompile.@snoop "glv_compiles.csv" begin
    include("test/ExampleRunner.jl")
    using ExampleRunner
    config = ExampleRunner.RunnerConfig(
        number_of_frames = 2,
        interactive_time = 0.01,
        record=false,
        resolution = (100, 100)
    )
    ExampleRunner.run(config)
end

using GLVisualize

str = open("glv_compiles.csv") do io
    str = readstring(io)
    x = replace(str, r"[0-9]+\t\"<toplevel thunk> [\s\S]+?\)::Any\]\"\n", "")
end
open("glv_compiles.csv", "w") do io
    seekstart(io)
    print(io, str)
end
data = SnoopCompile.read("glv_compiles.csv")
blacklist = [
    "MIME", "Base", "Core",
    "ExampleRunner", "ImageMagick", "Contour", "MeshIO"
]
pc = SnoopCompile.format_userimg(data[end:-1:1,2], blacklist=blacklist)
SnoopCompile.write(GLVisualize.dir("precompile", "glv_userimg.jl"), pc)

#=
real	0m13.712s
user	0m15.872s
sys	0m0.176s
=#

include(joinpath(JULIA_HOME, "..", "..", "contrib", "build_sysimg.jl"))
build_sysimg(default_sysimg_path(), "native", nothing; force=true)

userimg_path = GLVisualize.dir("precompile", "userimg.jl")

build_sysimg(default_sysimg_path(), "native", userimg_path; force=true)


text:
module TestRegInfo

using MLStyle
using YaoLocations
using YaoCompiler
using YaoCompiler.Intrinsics
using YaoTargetQASM
using CompilerPluginTools
using Test

@device function test_basic(theta, phi)
    # syntax sugar
    1 => X
    @gate 2 => Z
    @ctrl 1 4 => Rx(theta)
    @ctrl 2 4 => Ry(phi)
    a = @measure 3
    # direct intrinsic
    apply(Y, 3)
    apply(X, 1, 4)
    c = measure((2, 3))
    return (a = a, b = c)
end

@device function test_pure_quantum()
    ret = @gate 1:4 => test_basic(1.0, 2.0)
    @ctrl 2 1 => Rx(2.2)
    return ret
end

@testset "RegInfo" begin
    op = test_pure_quantum()
    interp = YaoInterpreter()
    ci, type = code_typed(Intrinsics.apply, (AnyReg, typeof(op)); interp)[1]
    ci = CompilerPluginTools.rm_code_coverage_effect(ci)
    info = RegInfo(OpenQASMTarget(), ci)

    @test info.creg_ssa_map[:a] == 5
    @test info.creg_ssa_map[:b] == 8
    @test info.ssa_creg_map[5] === :a
    @test info.ssa_creg_map[8] === :b
    @test info.creg_size[:a] == 1
    @test info.creg_size[:b] == 2
    @test info.qreg_to_locs[1] == [1, 4]
    @test info.qreg_to_locs[2] == [2, 3]
    @test info.locs_to_qreg[1] == 1
    @test info.locs_to_qreg[2] == 2
    @test info.locs_to_qreg[3] == 2
    @test info.locs_to_qreg[4] == 1
    @test info.locs_to_addr[1] == 0
    @test info.locs_to_addr[2] == 0
    @test info.locs_to_addr[3] == 1
    @test info.locs_to_addr[4] == 1

    buf = IOBuffer()
    show(buf, MIME"text/plain"(), info)
    s = String(take!(buf))
    @test occursin("RegInfo\n", s)
    @test occursin("creg:\n", s)
    @test occursin("  a[1]\n", s)
    @test occursin("  b[2]\n", s)
    @test occursin("qreg:\n", s)
    @test occursin("  %2[2] => [2, 3]", s)
    @test occursin("  %1[2] => [1, 4]", s)
end

end


text:
using Test

ns = [14, 1, 17, 0, 3, 20]

function saidwhen(ns)
    return Dict([(n, i) for (i, n) in enumerate(ns)])
end

function play(sofar, target)
    ns = copy(sofar)
    said = saidwhen(ns)
    turn = length(sofar)

    while turn < target
        last = ns[turn]
        if !haskey(said, last)
            append!(ns, 0)
            said[last] = turn
        else
            append!(ns, turn - said[last])
            said[last] = turn
        end
        turn += 1
    end
    return ns
end

@test last(play([0, 3, 6], 10)) == 0
@test last(play([1, 3, 2], 2020)) == 1
@test last(play([2, 1, 3], 2020)) == 10
@test last(play([1, 2, 3], 2020)) == 27
@test last(play([2, 3, 1], 2020)) == 78
@test last(play([3, 2, 1], 2020)) == 438
@test last(play([3, 1, 2], 2020)) == 1836

p1 = last(play(ns, 2020))
p2 = last(play(ns, 30000000))

@assert(p1 == 387)
@assert(p2 == 6428)

println("-----------------------------------------------------------------------")
println("rambunctious recitation -- part one :: $p1")
println("rambunctious recitation -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
using Documenter, Winds

makedocs(
    modules = [Winds],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "Winds.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/yakir12/Winds.jl.git",
)


text:
"""
    solve!(
        KS::X,
        ctr::Y,
        face::Z,
        simTime,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{<:AbstractControlVolume,1},
        Z<:AbstractArray{<:AbstractInterface1D,1},
    }

    solve!(
        KS::X,
        ctr::Y,
        a1face::Z,
        a2face::Z,
        simTime,
    ) where {
        X<:AbstractSolverSet,
        Y<:AbstractArray{<:AbstractControlVolume,2},
        Z<:AbstractArray{<:AbstractInterface2D,2},
    }

Solution algorithm

- @args: solver setup
- @args: array of control volumes
- @args: array of interfaces
- @args & return: time instant

"""
function solve!(
    KS::X,
    ctr::Y,
    face::Z,
    simTime,
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{<:Union{ControlVolume1D,ControlVolume1D1F,ControlVolume1D2F},1},
    Z<:AbstractArray{<:AbstractInterface1D,1},
}

    #--- initial checkpoint ---#
    write_jld(KS, ctr, simTime)

    #--- setup ---#
    iter = 0
    t = deepcopy(simTime)
    dt = timestep(KS, ctr, simTime)
    nt = Int(floor(KS.set.maxTime / dt)) + 1
    res = zero(ctr[1].w)

    #--- main loop ---#
    #while true
    @showprogress for iter = 1:nt

        #dt = timestep(KS, ctr, simTime)
        reconstruct!(KS, ctr)
        evolve!(
            KS,
            ctr,
            face,
            dt;
            mode = symbolize(KS.set.flux),
            bc = symbolize(KS.set.boundary),
        )
        update!(
            KS,
            ctr,
            face,
            dt,
            res;
            coll = symbolize(KS.set.collision),
            bc = symbolize(KS.set.boundary),
        )

        #iter += 1
        t += dt

        if iter % 500 == 0
            println("iter: $(iter), time: $(t), dt: $(dt), res: $(res)")

            #if iter%1000 == 0
            #    write_jld(KS, ctr, iter)
            #end
        end

        if t > KS.set.maxTime || maximum(res) < 5.e-7
            break
        end

    end

    write_jld(KS, ctr, simTime)
    return t

end

function solve!(
    KS::X,
    ctr::Y,
    a1face::Z,
    a2face::Z,
    simTime,
) where {
    X<:AbstractSolverSet,
    Y<:AbstractArray{<:AbstractControlVolume,2},
    Z<:AbstractArray{<:AbstractInterface2D,2},
}

    #--- initial checkpoint ---#
    write_jld(KS, ctr, simTime)

    #--- setup ---#
    iter = 0
    t = deepcopy(simTime)
    dt = timestep(KS, ctr, simTime)
    nt = Int(floor(KS.set.maxTime / dt)) + 1
    res = zero(ctr[1].w)

    #--- main loop ---#
    @showprogress for iter = 1:nt
        reconstruct!(KS, ctr)
        evolve!(
            KS,
            ctr,
            a1face,
            a2face,
            dt;
            mode = symbolize(KS.set.flux),
            bc = symbolize(KS.set.boundary),
        )
        update!(
            KS,
            ctr,
            a1face,
            a2face,
            dt,
            res;
            coll = symbolize(KS.set.collision),
            bc = symbolize(KS.set.boundary),
        )

        t += dt

        if iter % 500 == 0
            println("iter: $(iter), time: $(t), dt: $(dt), res: $(res[1:end])")
        end

        if t > KS.set.maxTime || maximum(res) < 5.e-7
            break
        end
    end

    write_jld(KS, ctr, simTime)
    return t

end

"""
    timestep(KS, ctr, simTime)

Calculate timestep based on current solutions

"""
function timestep(
    KS::X,
    ctr::Y,
    simTime,
) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,1}}

    tmax = 0.0

    if ctr[1].w isa Number

        @inbounds Threads.@threads for i = 1:KS.pSpace.nx
            prim = ctr[i].prim
            vmax = abs(ctr[i].prim[2])
            tmax = max(tmax, vmax / KS.ps.dx[i])
        end

    elseif KS.set.nSpecies == 1

        @inbounds Threads.@threads for i = 1:KS.pSpace.nx
            prim = ctr[i].prim
            sos = sound_speed(prim, KS.gas.γ)
            vmax = begin
                if KS.vs isa Nothing
                    abs(prim[2]) + sos
                else
                    max(KS.vSpace.u1, abs(prim[2])) + sos
                end
            end
            tmax = max(tmax, vmax / KS.ps.dx[i])
        end

    elseif KS.set.nSpecies == 2

        @inbounds Threads.@threads for i = 1:KS.pSpace.nx
            prim = ctr[i].prim
            sos = sound_speed(prim, KS.gas.γ)
            vmax = begin
                if KS.vs isa Nothing
                    maximum(abs.(prim[2, :])) + sos
                else
                    max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :]))) + sos
                end
            end

            if KS.set.space[3:4] in ["3f", "4f"]
                tmax = max(tmax, vmax / KS.ps.dx[i], KS.gas.sol / KS.ps.dx[i])
            else
                tmax = max(tmax, vmax / KS.ps.dx[i])
            end
        end

    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt

end

function timestep(
    KS::X,
    ctr::Y,
    simTime,
) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,2}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    tmax = 0.0

    if KS.set.nSpecies == 1

        @inbounds Threads.@threads for j = 1:ny
            for i = 1:nx
                prim = ctr[i, j].prim
                sos = sound_speed(prim, KS.gas.γ)
                umax = ifelse(
                    KS.vs isa Nothing,
                    abs(prim[2]) + sos,
                    max(KS.vSpace.u1, abs(prim[2])) + sos,
                )
                vmax = ifelse(
                    KS.vs isa Nothing,
                    abs(prim[3]) + sos,
                    max(KS.vSpace.v1, abs(prim[3])) + sos,
                )
                tmax = max(tmax, umax / dx[i, j] + vmax / dy[i, j])
            end
        end

    elseif KS.set.nSpecies == 2

        @inbounds Threads.@threads for j = 1:ny
            for i = 1:nx
                prim = ctr[i, j].prim
                sos = sound_speed(prim, KS.gas.γ)
                umax = max(maximum(KS.vSpace.u1), maximum(abs.(prim[2, :]))) + sos
                vmax = max(maximum(KS.vSpace.v1), maximum(abs.(prim[3, :]))) + sos

                if KS.set.space[3:4] in ["3f", "4f"]
                    tmax = max(
                        tmax,
                        umax / dx[i, j] + vmax / dy[i, j],
                        KS.gas.sol / dx[i, j] + KS.gas.sol / dy[i, j],
                    )
                else
                    tmax = max(tmax, umax / dx[i, j] + vmax / dy[i, j])
                end
            end
        end

    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt

end

function timestep(
    KS::X,
    ctr::Y,
    simTime,
) where {X<:AbstractSolverSet,Y<:AbstractVector{<:AbstractUnstructControlVolume}}

    tmax = 0.0

    if KS.set.nSpecies == 1

        @inbounds Threads.@threads for i in eachindex(ctr)
            prim = ctr[i].prim
            sos = sound_speed(prim, KS.gas.γ)
            umax = KS.vSpace.u1 + sos
            vmax = KS.vSpace.v1 + sos
            pd =
                abs.([dot([umax, vmax], ctr[i].n[j]) for j in eachindex(ctr[i].n)]) ./ ctr[i].dx
            #tmax = max(tmax, maximum(pd))
            tmax = max(tmax, sqrt(umax^2 + vmax^2) / minimum(ctr[i].dx))
        end

    elseif KS.set.nSpecies == 2

    end

    dt = KS.set.cfl / tmax
    dt = ifelse(dt < (KS.set.maxTime - simTime), dt, KS.set.maxTime - simTime)

    return dt

end


text:
using Plots, Printf, LinearAlgebra

# enable plotting by default
if !@isdefined do_visu; do_visu = true end

@views function diffusion_1D(; do_visu=true)
    # Physics
    lx     = 10.0        # domain size
    D      = 1.0         # diffusion coefficient
    ttot   = 0.6         # total simulation time
    # Numerics
    nx     = 256         # numerical grid resolution
    # Derived numerics
    dx     = lx/nx       # grid size
    dt     = dx^2/D/2.1  # time step (obeys CFL condition)
    xc     = LinRange(dx/2, lx-dx/2, nx)
    # Array allocation
    qH     = zeros(nx-1) # on staggered grid
    dHdt   = zeros(nx-2) # normal grid, without boundary points
    # Initial condition
    H0     = exp.(-(xc.-lx/2).^2)
    H      = copy(H0)
    t = 0.0; it = 0
    # Physical time loop
    while t<ttot
        qH         .= -D*diff(H)/dx         # flux
        dHdt       .=  -diff(qH)/dx         # rate of change
        H[2:end-1] .= H[2:end-1] .+ dt*dHdt # update rule, sets the BC as H[1]=H[end]=0
        t += dt; it += 1
    end
    # Analytic solution
    Hana = 1/sqrt(4*(ttot+1/4)) * exp.(-(xc.-lx/2).^2 /(4*(ttot+1/4)))
    @printf("Total time = %1.2f, time steps = %d, error vs analytic = %1.2e \n", round(ttot, sigdigits=2), it, norm(H-Hana))
    # Visualize
    if do_visu
        plot(xc, H0, linewidth=3); display(plot!(xc, H, legend=false, framestyle=:box, linewidth=3, xlabel="lx", ylabel="H", title="explicit diffusion (nt=$it)"))
    end
    return xc, H0
end

diffusion_1D(; do_visu=do_visu);


text:
using DelimitedFiles, PlotlyJS
import Dates

# Important: Library Plot.jl shall _not_ be loaded!

timestamp = Dates.format(Dates.now(), "yyyy-mm-dd_HH_MM_SS");

# Files: file/folder names and paths
basepath        = "results/";
plotsfolder     = "plots/"
input_filename  = "input_history_scenario_";
output_filename = "output_history_scenario_";
temperature2D_filename = "temperature2D_scenario_";

const L = 300.0*10^(-3) # Length of rod
const H = 10.0*10^(-3)  # Height of rod

# Discretization: necessary to import the matrices correctly
const N₁ = 100; # 50;
const N₂ = 40; #20;
const dx₁ = L/(N₁ - 1);
const dx₂ = H/(N₂ - 1);
const x1span = dx₁/2 : dx₁ : L;
const x2span = (dx₂/2 : dx₂ : H) * 1000;

# Number of scenarios
const N_scenarios = 2;

# Temperature on topside
temperature_on_top = zeros(N₁, N_scenarios);

for i = 1 : N_scenarios

    # Input and Output temporal data
    input_filepath  = basepath*input_filename*string(i)*".txt"
    inputdata       = readdlm(input_filepath, '\t', Float64)
    num_inputsignals = size(inputdata)[2] - 1
    input_averaged   = sum(inputdata[:, 2:end], dims=2) / num_inputsignals;

    output_filepath = basepath*output_filename*string(i)*".txt"
    outputdata      = readdlm(output_filepath, '\t', Float64)
    num_outputsignals   = size(outputdata)[2] - 1
    output_averaged     = sum(outputdata[:, 2:end], dims=2) / num_outputsignals;
    output_variance     = sum((outputdata[:, 2:end] .- output_averaged).^2, dims=2)

    input_output_plotname = basepath * plotsfolder  * "input_output_"*string(i)*"_"*string(timestamp) * ".pdf"

    output_plot = Plot(outputdata[:,1], output_averaged, Layout( yaxis = attr(; title="Averaged y", titlefont=attr(; size=16))) )
    input_plot  = Plot(inputdata[:,1], input_averaged, Layout( yaxis = attr(; title="Averaged u", titlefont=attr(; size=16))) )
    output_var_plot = Plot(outputdata[:,1], output_variance, Layout( xaxis = attr(; title="Time [s]", titlefont=attr(; size=20) ) ,yaxis=attr(;title="Variance y", titlefont=attr(; size=16)) ))

    in_out_plot = [output_plot; input_plot; output_var_plot]
    in_out_plot.layout["showlegend"] = false;
    in_out_plot.layout["margin"][:l] = 70;

    PlotlyJS.savefig(in_out_plot, input_output_plotname)

    # Temperature distribution: 2 dimensional contour plot
    temperature2D_filepath  = basepath*temperature2D_filename*string(i)*".txt"
    temperature2D_data      = readdlm(temperature2D_filepath, '\t', Float64 )
    data    = contour(; x=x1span, y=x2span, z=temperature2D_data, colorbar=attr(;title="Temperature [K]", titleside="right", titlefont=attr(;size=20)))
    layout  = Layout(; xaxis = attr(; title="Length [m]", titlefont=attr(; size=20) ) ,yaxis=attr(;title="Heigth [mm]", titlefont=attr(; size=20)) )
    temperature2D_plot = Plot(data, layout)

    temperature2D_plotpath = basepath * plotsfolder * "temperature2D_contour_" * string(i) * "_" *string(timestamp) * ".pdf";
    PlotlyJS.savefig(temperature2D_plot, temperature2D_plotpath)

    temperature_on_top[:,i] = temperature2D_data[:, end];

end

#Plotting the temperature on top of the plate
trace1 = scatter(; x=x1span, y=temperature_on_top[:,1], mode="lines", name="Scenario 1")
trace2 = scatter(; x=x1span, y=temperature_on_top[:,2], mode="lines", name="Scenario 2")

layout = Layout(; xaxis = attr(; title="Length [m]", titlefont=attr(; size=20), range=[-0.005, 0.301] ) ,yaxis=attr(;title="Temperature [K]", titlefont=attr(; size=20), range=[388, 404]) )
topplot = Plot([trace1, trace2], layout)

topplot_path = basepath * plotsfolder * "top_" * string(timestamp) * ".pdf"

PlotlyJS.savefig(topplot, topplot_path)


text:
using .ForwardDiff: ForwardDiff, DiffResults, StaticArrays

struct ForwardDiffBackend <: AbstractBackend end

@primitive function pushforward_function(ba::ForwardDiffBackend, f, xs...)
    return function pushforward(vs)
        if length(xs) == 1
            v = vs isa Tuple ? only(vs) : vs
            return (ForwardDiff.derivative(h -> f(step_toward(xs[1], v, h)), 0),)
        else
            return ForwardDiff.derivative(h -> f(step_toward.(xs, vs, h)...), 0)
        end
    end
end

primal_value(x::ForwardDiff.Dual) = ForwardDiff.value(x)

# these implementations are more efficient than the fallbacks

gradient(::ForwardDiffBackend, f, x::AbstractArray) = (ForwardDiff.gradient(f, x),)

function jacobian(ba::ForwardDiffBackend, f, x::AbstractArray)
    return (ForwardDiff.jacobian(asarray ∘ f, x),)
end
jacobian(::ForwardDiffBackend, f, x::Number) = (ForwardDiff.derivative(f, x),)

hessian(::ForwardDiffBackend, f, x::AbstractArray) = (ForwardDiff.hessian(f, x),)

function value_and_gradient(::ForwardDiffBackend, f, x::AbstractArray)
    result = ForwardDiff.gradient!(DiffResults.GradientResult(x), f, x)
    return DiffResults.value(result), (DiffResults.derivative(result),)
end

function value_and_hessian(::ForwardDiffBackend, f, x)
    result = ForwardDiff.hessian!(DiffResults.HessianResult(x), f, x)
    return DiffResults.value(result), (DiffResults.hessian(result),)
end

@inline step_toward(x::Number, v::Number, h) = x + h * v
# support arrays and tuples
@noinline step_toward(x, v, h) = x .+ h .* v

@inline asarray(x) = [x]
@inline asarray(x::AbstractArray) = x


text:
mutable struct Bisect{V,E} <: State{V}
    value::V
    step::Symbol
    N::Int
    a::V
    b::V
    c::V
    d::V
    fa::E
    fb::E
    fc::E
end

Bisect(; unit, evalunit, _type, _...) = begin
    V = valuetype(_type, value(unit))
    E = valuetype(_type, value(evalunit))
    v = zero(V)
    e = zero(E)
    Bisect{V,E}(v, :z, 0, v, v, v, v, e, e, e)
end

constructortags(::Val{:Bisect}) = (:unit, :evalunit)

@generated evalunit(s::Bisect{V,E}) where {V,E} = unittype(E)

updatetags!(d, ::Val{:Bisect}; _...) = begin
    !haskey(d, :lower) && haskey(d, :min) && (d[:lower] = d[:min])
    !haskey(d, :upper) && haskey(d, :max) && (d[:upper] = d[:max])
    !haskey(d, :evalunit) && (d[:evalunit] = d[:unit])
    !haskey(d, :maxiter) && (d[:maxiter] = 150)
    !haskey(d, :tol) && (d[:tol] = 0.00001)
end

genvartype(v::VarInfo, ::Val{:Bisect}; N, V, _...) = begin
    EU = gettag(v, :evalunit)
    E = @q $C.valuetype($N, $EU)
    @q Bisect{$V,$E}
end

gendefault(v::VarInfo, ::Val{:Bisect}) = nothing

#HACK: needs update in case min/max variables changed during bisection loop
#TODO: other variables wanting to use min/max would require similar work
genupdate(v::VarInfo, ::Val{:Bisect}, ::PreStep; kw...) = begin
    @gensym s d
    @q let $s = $(symstate(v)),
           $d = $(genminmax(v, @q $C.value($s)))
        $C.store!($s, $d)
    end
end
genupdate(v::VarInfo, ::Val{:Bisect}, ::MainStep; kw...) = begin
    maxiter = gettag(v, :maxiter)
    tol = gettag(v, :tol)
    lstart = symlabel(v, PreStep())
    lrepeat = symlabel(v, MainStep(), :__repeat)
    lexit = symlabel(v, MainStep(), :__exit)
    @gensym s u Δ
    @q let $s = $(symstate(v))
        if $s.step == :z
            $s.N = 1
            $u = $C.value($(gettag(v, :unit)))
            $s.a = $(genminmax(v, @q $C.unitfy($C.value($(gettag(v, :lower))), $u)))
            $s.b = $(genminmax(v, @q $C.unitfy($C.value($(gettag(v, :upper))), $u)))
            $s.c = zero($s.c)
            $s.d = $s.b - $s.a
            $s.fa = zero($s.fa) * NaN
            $s.fb = zero($s.fb) * NaN
            $s.fc = zero($s.fc) * NaN
            $s.step = :a
            $C.store!($s, $s.a)
            @goto $lrepeat
        elseif $s.step == :a
            $s.fa = $(genbody(v))
            if $s.fa ≈ zero($s.fa)
                @debug "bisect: $($v.name) [$($s.N)] found! $($C.value($s))"
                @goto $lexit
            elseif isnan($s.fa)
                @warn "bisect: $($v.name) [$($s.N)] $($s.a) => $($s.fa)"
            else
                @debug "bisect: $($v.name) [$($s.N)] $($s.a) => $($s.fa)"
            end
            $s.step = :b
            $C.store!($s, $s.b)
            @goto $lrepeat
        elseif $s.step == :b
            $s.fb = $(genbody(v))
            if $s.fb ≈ zero($s.fb)
                @debug "bisect: $($v.name) [$($s.N)] found! $($C.value($s))"
                @goto $lexit
            elseif isnan($s.fb)
                @warn "bisect: $($v.name) [$($s.N)] $($s.b) => $($s.fb)"
            else
                @debug "bisect: $($v.name) [$($s.N)] $($s.b) => $($s.fb)"
            end
            if sign($s.fa) == sign($s.fb)
                #HACK: try expanding bracket
                #$s.N += round(Int, 0.1*$maxiter)
                $Δ = $s.d / 2
                #HACK: reduce redundant unitfy when generating min/max clipping
                #TODO: check no expansion case where Δ gets clipped by min/max
                $s.a = $(genminmax(v, @q $s.a - $Δ))
                $s.b = $(genminmax(v, @q $s.b + $Δ))
                @debug "bisect: $($v.name) [$($s.N)] $($s.a) <- a, b -> $($s.b) "
                $s.step = :a
                $C.store!($s, $s.a)
                @goto $lrepeat
            end
            $s.c = ($s.a + $s.b) / 2
            $C.store!($s, $s.c)
            $s.step = :c
            @goto $lrepeat
        elseif $s.step == :c
            $s.fc = $(genbody(v))
            if isnan($s.fc)
                @warn "bisect: $($v.name) [$($s.N)] $($s.c) => $($s.fc)"
            else
                @debug "bisect: $($v.name) [$($s.N)] $($s.c) => $($s.fc)"
            end
            if $s.fc ≈ zero($s.fc) || ($s.b - $s.a) / $s.d < $tol
                @debug "bisect: $($v.name) [$($s.N)] finished! $($C.value($s))"
                @goto $lexit
            end
            if sign($s.fc) == sign($s.fa)
                $s.a = $s.c
                $s.fa = $s.fc
                @debug "bisect: $($v.name) [$($s.N)] a <- $($s.c)"
            else
                $s.b = $s.c
                $s.fb = $s.fc
                @debug "bisect: $($v.name) [$($s.N)] b <- $($s.c)"
            end
            $s.c = ($s.a + $s.b) / 2
            $C.store!($s, $s.c)
            @goto $lrepeat
        end
        @label $lrepeat
        if $s.N <= $maxiter
            $s.N += 1
            @goto $lstart
        else
            @warn "bisect: $($v.name) [$($s.N)] convergence failed!" a=$s.a b=$s.b c=$s.c fa=$s.fa fb=$s.fb fc=$s.fc d=$s.d $(v.name)=$C.value($s)
        end
        @label $lexit
        $s.step = :z
        $C.value($s)
    end
end


text:
#=
# test a simple quad mesh from a file
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("quadtest");

@useLog("quadtestlog")

@domain(2)
@functionSpace(LEGENDRE, 1)

# First make a structured mesh
@mesh(QUADMESH, 10)

# write a msh file for it
@outputMesh("quadmesh.msh")

# read it back in as an unstructured mesh
@mesh("quadmesh.msh")

@variable(u)
@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)

@coefficient(f, "-2*pi*pi*sin(pi*x)*sin(pi*y)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

solve(u);

@finalize()

# exact solution is sin(pi*x)*sin(pi*y)
# check error
maxerr = 0;
exact(x,y) = sin(pi*x)*sin(pi*y);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    y = Femshop.grid_data.allnodes[2,i];
    err = abs(u.values[i] - exact(x,y));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# Plot the mesh to check
# m = Femshop.mesh_data;
# using Plots
# pyplot();
# quadx = zeros(4, m.nel);
# quady = zeros(4, m.nel);
# centers = zeros(2, m.nel);
# for i=1:m.nel
#     quadx[:,i] = m.nodes[1,m.elements[1:4,i]];
#     quady[:,i] = m.nodes[2,m.elements[1:4,i]];
#     centers[1,i] = sum(quadx[:,i])/4;
#     centers[2,i] = sum(quady[:,i])/4;
# end
# p1 = plot(quadx, quady, legend=false)
# for i=1:m.nel
#     annotate!(centers[1,i], centers[2,i], string(i));
# end
# display(p1)


text:
# Bogumił Kamiński, 2022

# Codes for chapter 9

# Code for section 9.1

using DataFrames
using CSV
using Plots
puzzles = CSV.read("puzzles.csv", DataFrame);

using Statistics
plays_lo = median(puzzles.NbPlays)
puzzles.NbPlays .> plays_lo

puzzles.NbPlays > plays_lo

rating_lo = 1500
rating_hi = quantile(puzzles.Rating, 0.99)
rating_lo .< puzzles.Rating .< rating_hi

row_selector = (puzzles.NbPlays .> plays_lo) .&&
               (rating_lo .< puzzles.Rating .< rating_hi)

sum(row_selector)
count(row_selector)

# Code for listing 9.1

good = puzzles[row_selector, ["Rating", "Popularity"]]

# Code for plotting histograms

plot(histogram(good.Rating; label="Rating"),
     histogram(good.Popularity; label="Popularity"))

# Code for column selectors

puzzles[1, "Rating"]

puzzles[:, "Rating"]

row1 = puzzles[1, ["Rating", "Popularity"]]

row1["Rating"]
row1[:Rating]
row1[1]
row1.Rating
row1."Rating"

good = puzzles[row_selector, ["Rating", "Popularity"]]

good[1, "Rating"]
good[1, :]
good[:, "Rating"]
good[:, :]

names(puzzles, ["Rating", "Popularity"])
names(puzzles, [:Rating, :Popularity])
names(puzzles, [4, 6])
names(puzzles, [false, false, false, true, false, true, false, false, false])
names(puzzles, r"Rating")
names(puzzles, Not([4, 6]))
names(puzzles, Not(r"Rating"))
names(puzzles, Between("Rating", "Popularity"))
names(puzzles, :)
names(puzzles, All())
names(puzzles, Cols(r"Rating", "NbPlays"))
names(puzzles, Cols(startswith("P")))

names(puzzles, startswith("P"))

names(puzzles, Real)

names(puzzles, AbstractString)

puzzles[:, names(puzzles, Real)]

# Code for row subsetting

df_small = DataFrame(id=1:4)
df_small[[1, 3], :]
df_small[[true, false, true, false], :]
df_small[Not([2, 4]), :]
df_small[Not([false, true, false, true]), :]

df1 = puzzles[:, ["Rating", "Popularity"]];
df2 = puzzles[!, ["Rating", "Popularity"]];

df1 == df2
df1 == puzzles
df2 == puzzles

df1.Rating === puzzles.Rating
df1.Popularity === puzzles.Popularity
df2.Rating === puzzles.Rating
df2.Popularity === puzzles.Popularity

using BenchmarkTools
@btime $puzzles[:, ["Rating", "Popularity"]];
@btime $puzzles[!, ["Rating", "Popularity"]];

puzzles[1, 1]
puzzles[[1], 1]
puzzles[1, [1]]
puzzles[[1], [1]]

# Code for making views

@view puzzles[1, 1]

@view puzzles[[1], 1]

@view puzzles[1, [1]]

@view puzzles[[1], [1]]

@btime $puzzles[$row_selector, ["Rating", "Popularity"]];
@btime @view $puzzles[$row_selector, ["Rating", "Popularity"]];

parentindices(@view puzzles[row_selector, ["Rating", "Popularity"]])

# Code for section 9.2

describe(good)

rating_mapping = Dict{Int, Vector{Int}}()
for (i, rating) in enumerate(good.Rating)
    if haskey(rating_mapping, rating)
        push!(rating_mapping[rating], i)
    else
        rating_mapping[rating] = [i]
    end
end
rating_mapping

good[rating_mapping[2108], :]

unique(good[rating_mapping[2108], :].Rating)

using Statistics
mean(good[rating_mapping[2108], "Popularity"])

ratings = unique(good.Rating)

mean_popularities = map(ratings) do rating
    indices = rating_mapping[rating]
    popularities = good[indices, "Popularity"]
    return mean(popularities)
end

scatter(ratings, mean_popularities;
        xlabel="rating", ylabel="mean popularity", legend=false)

using Loess
model = loess(ratings, mean_popularities);
ratings_predict = float.(sort(ratings))
popularity_predict = predict(model, ratings_predict)

plot!(ratings_predict, popularity_predict; width=5, color="black")

combine(groupby(good, :Rating), :Popularity => mean)

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis