text:
# Our developed libraries
using CombinatorialSpaces
using CombinatorialSpaces.ExteriorCalculus

using Catlab.Present
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using LinearAlgebra

using Decapodes.Simulations
using Decapodes.Examples
using Decapodes.Diagrams
using Decapodes.Schedules

# Julia community libraries
using MeshIO
using CairoMakie
using Decapodes.Debug
using DifferentialEquations
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

# Define used quantities
@present Flow2DQuantities(FreeExtCalc2D) begin
  X::Space
  C::Hom(munit(), Form0(X))     # concentration
  dC::Hom(munit(), Form0(X))     # concentration
  V::Hom(munit(), Form1(X))     # Flow field
  dV::Hom(munit(), Form1(X))    # Flow field time derivative
  P::Hom(munit(), DualForm1(X)) # Flow momentum
  p::Hom(munit(), Form0(X))     # pressure field
  dp::Hom(munit(), Form0(X))     # pressure field time derivative
  ϕ::Hom(munit(), DualForm1(X)) # negative diffusion flux
  k::Hom(Form1(X), Form1(X))    # diffusivity (usually scalar multiplication)
  kᵥ::Hom(Form1(X), Form1(X))    # viscosity (usually scalar multiplication)
  kₚ::Hom(Form0(X), Form0(X))    # compressibility (usually scalar multiplication)
  m⁻¹::Hom(DualForm1(X), DualForm1(X))    # diffusivity (usually scalar multiplication)
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
  L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
  ∂₀::Hom(Form0(X), Form0(X)) # all_wall boundary condition
  ∂₀₋::Hom(Form0(X), Form0(X)) # left/right wall boundary condition
  ∂₀₊::Hom(Form0(X), Form0(X)) # top/bottom wall boundary condition
  ∂₀ₛ::Hom(Form0(X), Form0(X)) # Sticking boundary condition
  ∂₁ₛ::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₗ₊::Hom(Form1(X), Form1(X)) # Sticking boundary condition
  ∂₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  dneg₁::Hom(DualForm1(X), DualForm1(X)) # In/Out edge flow boundary condition
  neg₁::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  dneg₀::Hom(DualForm0(X), DualForm0(X)) # In/Out edge flow boundary condition
  neg₀::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition
  mask₁ₑ::Hom(Form1(X), Form1(X)) # In/Out edge flow boundary condition
  mask₀ₗ::Hom(Form0(X), Form0(X)) # In/Out edge flow boundary condition

	# Boundary Conditions (currently does not impact actual BCs)
  bc₀::Hom(munit(), Form0(X))
  bc₁::Hom(munit(), Form1(X))
  bc₂::Hom(munit(), Form1(X))
  bc₃::Hom(munit(), Form0(X))
end

# Define Diffusion/Advection physics

@present DiffAdv2D <: Flow2DQuantities begin
  # Fick's first law
  ϕ == C ⋅ d₀(X) ⋅ k ⋅ ⋆₁(X)
  # Diffusion/advection equation
  dC == ϕ ⋅ dual_d₁(X) ⋅ ⋆₀⁻¹(X) + (V⊗(C⋅⋆₀(X))) ⋅ L₀ ⋅ ⋆₀⁻¹(X) ⋅ neg₀
  C ⋅ ∂ₜ(Form0(X)) == dC
  # Boundary condition
  C ⋅ ∂₀ == bc₀
end

# Extend to incompressible Navier Stokes

@present IncompressibleFlow2D <: DiffAdv2D begin
  P == V ⋅ ⋆₁(X)
  dV == (V ⊗ P) ⋅ L₁ ⋅ dneg₁ ⋅ ⋆₁⁻¹(X) ⋅ mask₁ₑ + V ⋅ Δ₁(X) ⋅ kᵥ + p ⋅ d₀(X) ⋅ neg₁
  V ⋅ ∂ₜ(Form1(X)) == dV
  dp == V ⋅ neg₁ ⋅ δ₁(X) ⋅ kₚ + (V ⊗ (p⋅⋆₀(X))) ⋅ L₀ ⋅ ⋆₀⁻¹(X) ⋅ mask₀ₗ ⋅ neg₀
  V ⋅ ∂₁ₛ == bc₁
  dC ⋅ ∂₀ₛ == bc₃
  dp ⋅ ∂₀₋ == bc₀
  dV ⋅ ∂₁ₗ₊ == bc₁
  p ⋅ ∂ₜ(Form0(X)) == dp
end

diag = eq_to_diagrams(IncompressibleFlow2D)
to_graphviz(diag ; edge_len="1.3")

dwd =diag2dwd(diag)
to_graphviz(dwd, orientation=LeftToRight)

exp_dwd = Examples.expand_dwd(dwd, gen_dec_rules())
to_graphviz(exp_dwd, orientation=LeftToRight)

# Load mesh

s = EmbeddedDeltaSet2D("../meshes/naca0012_8.stl")
sd = dual(s);

# Define non-default operators (multiplication by constants, boundary conditions)
k = 0.1
kₚ = 10.0
kᵥ = 1.0
funcs = sym2func(sd)
∂₀ = Examples.boundary_inds(Val{0}, s)
∂₁ = Examples.boundary_inds(Val{1}, s)

∂ₗ₀ = ∂₀[findall(p-> p[1] <= -50.0, s[∂₀, :point])]
∂ᵣ₀ = ∂₀[findall(p-> p[1] >= 50.0, s[∂₀, :point])]
∂ₜ₀ = ∂₀[findall(p-> p[2] >= 15.0, s[∂₀, :point])]
∂ᵦ₀ = ∂₀[findall(p-> p[2] <= -15.0, s[∂₀, :point])]
∂ₑ₀ = vcat(∂ₗ₀, ∂ᵣ₀, ∂ₜ₀, ∂ᵦ₀)
∂ₒ₀ = ∂₀[findall(p-> -49 <= p[1] <= 49 && -14 <= p[2] <= 14, s[∂₀, :point])]

∂ₗ₁ = Examples.bound_edges(s, ∂ₗ₀)
∂ᵣ₁ = Examples.bound_edges(s, ∂ᵣ₀)
∂ₗ₁₊ = Examples.adj_edges(s, ∂ₗ₀)
∂ᵣ₁₊ = Examples.adj_edges(s, ∂ᵣ₀)
∂ₒ₁ = Examples.bound_edges(s, ∂ₒ₀)

∂ₗ₀₊ = unique(vcat(s[∂ₗ₁₊, :tgt], s[∂ₗ₁₊, :src]))

funcs[:mcopy] = Dict(:operator => I(ne(sd)), :type => MatrixFunc())
funcs[:k] = Dict(:operator => k * I(ne(sd)), :type => MatrixFunc())
funcs[:half] = Dict(:operator => 0.5 * I(ne(sd)), :type => MatrixFunc())
funcs[:kₚ] = Dict(:operator => kₚ * I(nv(sd)), :type => MatrixFunc())
funcs[:kᵥ] = Dict(:operator => kᵥ * I(ne(sd)), :type => MatrixFunc())
funcs[:dneg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:neg₁] = Dict(:operator => -1 * I(ne(sd)), :type => MatrixFunc())
funcs[:neg₀] = Dict(:operator => -1 * I(nv(sd)), :type => MatrixFunc())
funcs[:∂₀] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₑ₀] .= 0), :type => InPlaceFunc())
funcs[:∂₀₋] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₀] .= 0; x′[∂ᵣ₀] .= 0), :type => InPlaceFunc())
funcs[:∂ₗ₀₊] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₀₊] .= 0), :type => InPlaceFunc())
funcs[:∂₁ₗ₊] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₁₊] .= 0), :type => InPlaceFunc())
funcs[:∂₀ₛ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₒ₀] .= 0), :type => InPlaceFunc())
funcs[:∂₁ₛ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₒ₁] .= 0), :type => InPlaceFunc())
funcs[:∂₁ₑ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₁₊] .= 0;  x′[∂ᵣ₁₊] .= 0), :type => InPlaceFunc())
funcs[:mask₁ₑ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₁₊] .= 0;  x′[∂ᵣ₁₊] .= 0), :type => InPlaceFunc())
funcs[:mask₀ₗ] = Dict(:operator => (x′,x) -> (x′ .= x; x′[∂ₗ₀₊] .= 0), :type => InPlaceFunc())

# Define initial conditions
c_objs = zeros(nv(s))
c_objs[∂ₒ₀] .= 1e-3
velocity(p) = [2.0 * exp(-(p[1]+50)/5),0.0,0.0]
v = ♭(sd, DualVectorField(velocity.(sd[triangle_center(sd),:dual_point]))).data;
p = [0.0 for p in s[:point]]
u0 = vcat(c_objs, v, p)
GC.gc()

fig, ax, ob = wireframe(s)
ax.aspect = AxisAspect(3)
fig
save("mesh.svg", fig)

# Compress wiring diagram over contiguous matrix multiplications
cont_dwd = deepcopy(exp_dwd)
Examples.contract_matrices!(cont_dwd, funcs)

# Generate simulation function
func, _ = gen_sim(cont_dwd, funcs, sd; autodiff=false);

# Solve problem
prob = ODEProblem(func, u0, (0, 30))
sol = solve(prob, Tsit5(), progress=true, progress_steps=100, p=v);

# Example of debugging simulation

# Key for debugging simulation
sim_key(exp_dwd, orientation=LeftToRight)

exp_func, _ = gen_sim(exp_dwd, funcs, sd; autodiff=false);

# Show the flow of concentration as arrows
fig, ax, ob = draw_wire(s, sd, exp_dwd, exp_func, sol(30.0), 52; axisaspect=1, xlim=(40, 50), ylim=(-5, 5), use_arrows=true, n_arrows=1000)
fig
save("debug_conc_flow.svg")

# Generate animations of simulation

times = range(0, sol.t[end], length=150)
colors = [sol(t)#=[(end-nv(s)+1):end]=# for t in times]

figure, axis, scatter_thing = mesh(s, color=colors[1],
                                   colorrange=(0,1e-3))#minimum(vcat(colors...)),maximum(vcat(colors...))))
axis.aspect = AxisAspect(100.0/30.0)
framerate = 30

record(figure, "conc_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
  scatter_thing.color = colors[i]
end

times = range(0, sol.t[end], length=150)
colors = [sol(t)[(end-nv(s)+1):end] for t in times]

figure, axis, scatter_thing = mesh(s, color=colors[1],
                                   colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))
axis.aspect = AxisAspect(100.0/30.0)
framerate = 30

record(figure, "press_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
  scatter_thing.color = colors[i]
end


text:
@testset "scaletransform" begin
    rng = MersenneTwister(123456)
    s = rand(rng) + 1e-3
    t = ScaleTransform(s)

    x = randn(rng, 7)
    XV = [randn(rng, 6) for _ in 1:4]
    XC = ColVecs(randn(rng, 10, 5))
    XR = RowVecs(randn(rng, 6, 11))

    @testset "$(typeof(x))" for x in [x, XV, XC, XR]
        x′ = map(t, x)
        @test all([t(x[n]) ≈ s .* x[n] for n in eachindex(x)])
        @test all([t(x[n]) ≈ x′[n] for n in eachindex(x)])
    end

    s2 = 2.0
    KernelFunctions.set!(t, s2)
    @test t.s == [s2]
    @test isequal(ScaleTransform(s), ScaleTransform(s))
    @test repr(t) == "Scale Transform (s = $(s2))"
    test_ADs(x -> SEKernel() ∘ ScaleTransform(exp(x[1])), randn(rng, 1))

    @testset "median heuristic" begin
        for x in (x, XV, XC, XR), dist in (Euclidean(), Cityblock())
            n = length(x)
            t = median_heuristic_transform(dist, x)
            @test t isa ScaleTransform
            @test first(t.s) ≈
                inv(median(dist(x[i], x[j]) for i in 1:n, j in 1:n if i != j))

            y = map(t, x)
            @test median(dist(y[i], y[j]) for i in 1:n, j in 1:n if i != j) ≈ 1
        end
    end
end

ext: jl


text:
using TimeSchemes
using Test

##
M = 1
B = 1
K = 1
v0 = [1.]
v1 = [0.]
t = convert(Array{Float64},range(0, stop=1, length=10001))'
F = zeros(length(t))'
v, cV = TimeSchemes.newmark(M, B, K, F, v0, v1, t)

sol = 1 / 3 * exp.(- t / 2) .* (sqrt(3) * sin.(sqrt(3) / 2 * t) + 3  * cos.(sqrt(3) / 2 * t))
sold = - 1 / 6 * exp.(- t / 2) .* (sqrt.(3) * sin.(sqrt(3) / 2 * t) + 3  * cos.(sqrt(3) / 2 * t)) +
    1 / 3 * exp.(- t / 2) .* (3 / 2 * cos.(sqrt(3) / 2 * t) - 3 * sqrt(3) / 2 * sin.(sqrt(3) / 2 * t))
@test (maximum(abs.(v[1, :]' - sol)) < 1e-9) && (maximum(abs.(v[2, :]' - sold)) < 1e-9)
println("All tests passed.")


text:
n = rand(1:10)
@test matrixdepot("fiedler", n) == matrixdepot("fiedler", [1:n;])

A = matrixdepot("fiedler", n)

@test issym(A)
for i = 1:n, j =1:n
    @test_approx_eq A[i,j] abs(i-j)
end
println("'fiedler' passed test...")


text:
# Helper functions for exploiting sparsity in the Jacobians of μ, ξ, and 𝒱

"""
```
compute_sparsity_pattern(f::Function, x::AbstractVector{<: Number}, nrow::Int;
                         sparsity::Union{AbstractArray, Nothing} = nothing,
                         sparsity_detection::Bool = false)
```
calculates the sparsity pattern of the Jacobian of the functions μ, ξ, and 𝒱.

### Inputs
- `f`: is the function to be differentiated, e.g. `z -> 𝒱(z, Ψ, (1, 2))`
- `x`: the vector at which differentiation occurs
- `nrow`: specifies the number of rows of the Jacobian

### Keywords
- `sparsity`: sparsity pattern of the Jacobian
- `sparsity_detection`: if true, use SparsityDetection.jl to determine the sparsity pattern.
    If false, then the sparsity pattern is determined by using automatic differentiation
    to calculate a Jacobian and assuming any zeros are supposed to be zero.
"""
function compute_sparsity_pattern(f::Function, x::AbstractVector{T}, nrow::Int;
                                  sparsity::Union{AbstractArray, Nothing} = nothing,
                                  sparsity_detection::Bool = false) where {T <: Number}

    if isnothing(sparsity)
        sparsity = if sparsity_detection
            convert(SparseMatrixCSC{eltype(x), Int}, jacobian_sparsity(f, similar(x, nrow), x))
        else
            sparse(ForwardDiff.jacobian(f, x))
        end
    end

    if isempty(nonzeros(sparsity))
        # default to differentiating a dense matrix if all zeros
        return sparse(ones(T, size(sparsity))), 1:length(x)
    else
        return sparsity, matrix_colors(sparsity)
    end
end

"""
```
update_sparsity_pattern!(m::RiskAdjustedLinearization, function_name::Union{Symbol, Vector{Symbol}};
                         z::AbstractVector{<: Number} = m.z,
                         y::AbstractVector{<: Number} = m.y,
                         Ψ::AbstractVector{<: Number} = m.Ψ,
                         sparsity::AbstractDict{Symbol, AbstractMatrix} = Dict{Symbol, AbstractMatrix}(),
                         colorvec::AbstractDict{Symbol, <: AbstractVector{Int}} = Dict{Symbol, Vector{Int}}(),
                         sparsity_detection::Bool = false)
```
updates the Jacobians of μ, ξ, and/or 𝒱 in `m` with a new sparsity pattern. The Jacobians
to be updated are specified by `function_name`, e.g. `function_name = [:μ, :ξ, :𝒱]`.

If the keyword `sparsity` is empty, then the function attempts to determine the new sparsity pattern by computing
the Jacobian via automatic differentiation and assuming any zeros are supposed to be zero.
Keywords provide guesses for the coefficients ``(z, y, \\Psi)`` that are required
to calculate the Jacobians.

### Keywords
- `z`: state coefficients at steady state
- `y`: jump coefficients at steady state
- `Ψ`: coefficients for mapping from states to jumps
- `sparsity`: key-value pairs can be used to specify new sparsity patterns for the Jacobian functions
    `μz`, `μy`, `ξz`, `ξy`, and `J𝒱 `.
- `colorvec`: key-value pairs can be used to specify new matrix coloring vectors for the Jacobian functions
    `μz`, `μy`, `ξz`, `ξy`, and `J𝒱 `.
- `sparsity_detection`: use SparsityDetection.jl to determine the sparsity pattern.
"""
function update_sparsity_pattern!(m::RiskAdjustedLinearization, function_name::Symbol;
                                  z::AbstractVector{<: Number} = m.z,
                                  y::AbstractVector{<: Number} = m.y,
                                  Ψ::AbstractMatrix{<: Number} = m.Ψ,
                                  sparsity::AbstractDict = Dict{Symbol, Matrix}(),
                                  colorvec::AbstractDict = Dict{Symbol, Vector{Int}}(),
                                  sparsity_detection::Bool = false)
    return update_sparsity_pattern!(m, [function_name]; z = z, y = y, Ψ = Ψ,
                                    sparsity = sparsity, colorvec = colorvec,
                                    sparsity_detection = sparsity_detection)
end

function update_sparsity_pattern!(m::RiskAdjustedLinearization, function_names::Vector{Symbol};
                                  z::AbstractVector{<: Number} = m.z,
                                  y::AbstractVector{<: Number} = m.y,
                                  Ψ::AbstractMatrix{<: Number} = m.Ψ,
                                  sparsity::AbstractDict = Dict{Symbol, Matrix}(),
                                  colorvec::AbstractDict = Dict{Symbol, Vector{Int}}(),
                                  sparsity_detection::Bool = false)

    if :μ in function_names
        μz, μy, μz_jac_cache, μy_jac_cache =
            construct_μ_jacobian_function(m.nonlinear.μ, z, y;
                                          sparsity_z = haskey(sparsity, :μz) ? sparsity[:μz] : nothing,
                                          sparsity_y = haskey(sparsity, :μy) ? sparsity[:μy] : nothing,
                                          colorvec_z = haskey(colorvec, :μz) ? colorvec[:μz] : nothing,
                                          colorvec_y = haskey(colorvec, :μy) ? colorvec[:μy] : nothing,
                                          sparsity_detection = sparsity_detection)

        m.linearization.μz = μz
        m.linearization.μy = μy
        m.linearization.sparse_jac_caches[:μz] = μz_jac_cache
        m.linearization.sparse_jac_caches[:μy] = μy_jac_cache
    end

    if :ξ in function_names
        ξz, ξy, ξz_jac_cache, ξy_jac_cache =
            construct_ξ_jacobian_function(m.nonlinear.ξ, z, y;
                                          sparsity_z = haskey(sparsity, :ξz) ? sparsity[:ξz] : nothing,
                                          sparsity_y = haskey(sparsity, :ξy) ? sparsity[:ξy] : nothing,
                                          colorvec_z = haskey(colorvec, :ξz) ? colorvec[:ξz] : nothing,
                                          colorvec_y = haskey(colorvec, :ξy) ? colorvec[:ξy] : nothing,
                                          sparsity_detection = sparsity_detection)

        m.linearization.ξz = ξz
        m.linearization.ξy = ξy
        m.linearization.sparse_jac_caches[:ξz] = ξz_jac_cache
        m.linearization.sparse_jac_caches[:ξy] = ξy_jac_cache
    end

    if :𝒱 in function_names
        J𝒱, J𝒱_jac_cache = if isa(m.nonlinear.𝒱, RALF2)
            construct_𝒱_jacobian_function(m.nonlinear.𝒱, m.nonlinear.ccgf, m.nonlinear.Λ, m.nonlinear.Σ,
                                          m[:Γ₅], m[:Γ₆], z, Ψ;
                                          sparsity = haskey(sparsity, :J𝒱) ? sparsity[:J𝒱] : nothing,
                                          colorvec = haskey(colorvec, :J𝒱) ? colorvec[:J𝒱] : nothing,
                                          sparsity_detection = sparsity_detection)
        else
            construct_𝒱_jacobian_function(m.nonlinear.𝒱, m.nonlinear.ccgf, m.nonlinear.Λ, m.nonlinear.Σ,
                                          m[:Γ₅], m[:Γ₆], z, y, Ψ; sparsity = haskey(sparsity, :J𝒱) ? sparsity[:J𝒱] : nothing,
                                          colorvec = haskey(colorvec, :J𝒱) ? colorvec[:J𝒱] : nothing,
                                          sparsity_detection = sparsity_detection)
        end

        m.linearization.J𝒱 = J𝒱
        m.linearization.sparse_jac_caches[:J𝒱] = J𝒱_jac_cache
    end

    m
end

## Helper functions for constructing the Jacobian functions of μ, ξ, and 𝒱
function construct_μ_jacobian_function(μ::RALF2, z::AbstractVector{T}, y::AbstractVector{T};
                                       sparsity_z::Union{AbstractArray, Nothing} = nothing,
                                       sparsity_y::Union{AbstractArray, Nothing} = nothing,
                                       colorvec_z = nothing, colorvec_y = nothing,
                                       sparsity_detection::Bool = false) where {T <: Number}

    # Define (temporary) objective functions
    _f_μz = z -> μ(z, y, (1, 2))
    _f_μy = y -> μ(z, y, (2, 3))

    # Infer sparsity patterns and matrix coloring vector
    Nz = length(z)
    if isnothing(sparsity_z)
        sparsity_z, colorvec_z = compute_sparsity_pattern(_f_μz, z, Nz; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec_z)
        colorvec_z = matrix_colors(sparsity_z)
    end
    if isnothing(sparsity_y)
        sparsity_y, colorvec_y = compute_sparsity_pattern(_f_μy, y, Nz; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec_y)
        colorvec_y = matrix_colors(sparsity_y)
    end

    #=        # Create caches for the sparse Jacobian methods # This code is left here for when
    jac_cache_μz = ForwardColorJacCache(_f_μz, z, min(m.Nz, m.Ny); # Jacobians of μ and ξ are refactored
    sparsity = sparsity_μz, colorvec = colorvec_μz)
    jac_cache_μy = ForwardColorJacCache(_f_μy, y, min(m.Nz, m.Ny);
    sparsity = sparsity_μy, colorvec = colorvec_μy)=#

    # Create RALF2 objects. Note that we cannot pre-allocate the caches for
    # forwarddiff_color_jacobian! by using ForwardColorJacCache b/c the objective function
    # changes as z and y change. If Jacobians of μ and ξ are refactored to be done once,
    # then it'll be possible to cache.
    μ_dz = similar(z)
    μ_dy = similar(z)

    μz = RALF2((F, z, y) -> forwarddiff_color_jacobian!(F, (F0, x) -> μ.f0(F0, x, y), z, dx = μ_dz,
                                                        colorvec = colorvec_z, sparsity = sparsity_z),
               z, y, deepcopy(sparsity_z))
    μy = RALF2((F, z, y) -> forwarddiff_color_jacobian!(F, (F0, x) -> μ.f0(F0, z, x), y, dx = μ_dy,
                                                        colorvec = colorvec_y, sparsity = sparsity_y),
               z, y, deepcopy(sparsity_y))

    # Create mini-version of the Jacobian cache
    μz_jac_cache = (dx = μ_dz, sparsity = sparsity_z, colorvec = colorvec_z)
    μy_jac_cache = (dx = μ_dy, sparsity = sparsity_y, colorvec = colorvec_y)

    return μz, μy, μz_jac_cache, μy_jac_cache
end

function construct_ξ_jacobian_function(ξ::RALF2, z::AbstractVector{T}, y::AbstractVector{T};
                                       sparsity_z::Union{AbstractArray, Nothing} = nothing,
                                       sparsity_y::Union{AbstractArray, Nothing} = nothing,
                                       colorvec_z = nothing, colorvec_y = nothing,
                                       sparsity_detection::Bool = false) where {T <: Number}

    # Define (temporary) objective functions
    _f_ξz = z -> ξ(z, y, (1, 2))
    _f_ξy = y -> ξ(z, y, (2, 3))

    # Infer sparsity patterns and matrix coloring vector
    Ny = length(y)
    if isnothing(sparsity_z)
        sparsity_z, colorvec_z = compute_sparsity_pattern(_f_ξz, z, Ny; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec_z)
        colorvec_z = matrix_colors(sparsity_z)
    end
    if isnothing(sparsity_y)
        sparsity_y, colorvec_y = compute_sparsity_pattern(_f_ξy, y, Ny; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec_y)
        colorvec_y = matrix_colors(sparsity_y)
    end

    #=        # Create caches for the sparse Jacobian methods # This code is left here for when
    jac_cache_ξz = ForwardColorJacCache(_f_ξz, z, min(m.Nz, m.Ny);
    sparsity = sparsity_ξz, colorvec = colorvec_ξz)
    jac_cache_ξy = ForwardColorJacCache(_f_ξy, y, min(m.Nz, m.Ny);
    sparsity = sparsity_ξy, colorvec = colorvec_ξy)=#

    # Create RALF2 objects. Note that we cannot pre-allocate the caches for
    # forwarddiff_color_jacobian! by using ForwardColorJacCache b/c the objective function
    # changes as z and y change. If Jacobians of μ and ξ are refactored to be done once,
    # then it'll be possible to cache.
    ξ_dz = similar(y)
    ξ_dy = similar(y)

    ξz = RALF2((F, z, y) -> forwarddiff_color_jacobian!(F, (F0, x) -> ξ.f0(F0, x, y), z, dx = ξ_dz,
                                                        colorvec = colorvec_z, sparsity = sparsity_z),
               z, y, deepcopy(sparsity_z))
    ξy = RALF2((F, z, y) -> forwarddiff_color_jacobian!(F, (F0, x) -> ξ.f0(F0, z, x), y, dx = ξ_dy,
                                                        colorvec = colorvec_y, sparsity = sparsity_y),
               z, y, deepcopy(sparsity_y))

    # Create mini-version of the Jacobian cache
    ξz_jac_cache = (dx = ξ_dz, sparsity = sparsity_z, colorvec = colorvec_z)
    ξy_jac_cache = (dx = ξ_dy, sparsity = sparsity_y, colorvec = colorvec_y)

    return ξz, ξy, ξz_jac_cache, ξy_jac_cache
end

function construct_𝒱_jacobian_function(𝒱::RALF2, ccgf::Function, Λ::RALF1, Σ::RALF1{LC}, Γ₅::AbstractArray{<: Number},
                                       Γ₆::AbstractArray{<: Number}, z::AbstractVector{T}, Ψ::AbstractMatrix{T};
                                       sparsity::Union{AbstractArray, Nothing} = nothing,
                                       colorvec = nothing, sparsity_detection::Bool = false) where {T <: Number, LC}

    # Define (temporary) objective functions
    _f_𝒱z = z -> 𝒱(z, Ψ, (1, 2))

    # Need to grab some dimensions
    Ny, Nz = size(Ψ)
    Nε     = size(LC <: AbstractArray ? Σ.cache : Σ.cache.du, 2)

    # Infer sparsity patterns and matrix coloring vector
    if isnothing(sparsity)
        sparsity, colorvec = compute_sparsity_pattern(_f_𝒱z, z, Ny; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec)
        colorvec = matrix_colors(sparsity)
    end

    # Create RALF2 object. Note that we cannot pre-allocate the caches for
    # forwarddiff_color_jacobian! by using ForwardColorJacCache b/c the objective function
    # changes as the coefficients (z, y, Ψ) change.
    #
    # For 𝒱 specifically, to avoid problems reinterpreting arrays to make autodiff work,
    # we redefine the 𝒱 function to use Λ.f0 and Σ.f0 rather than using the RALF objects
    # (as we do when applying autodiff for dense Jacobians). If we use Λ(z) and Σ(z)
    # directly, then the reinterpret step may either fail (cannot reinterpret the array
    # to the new desired chunk size), or the reinterpreted array
    # will have the wrong dimensions.

    ## Infer whether Λ and Σ are in place
    FΛ0 = similar(z, Nz, Ny)
    FΣ0 = similar(z, Nz, Nε)
    Λ0 = if applicable(Λ.f0, FΛ0, z)
        Λ.f0
    else
        function _Λ_op(F, z)
            F .= Λ.f0(z)
        end
    end
    Σ0 = if applicable(Σ.f0, FΛ0, z)
        Σ.f0
    else
        function _Σ_op(F, z)
            F .= Σ.f0(z)
        end
    end

    _𝒱_sparse = if applicable(ccgf, Γ₅, z)
        function _𝒱_sparse_op(F, z, Ψ)
            FΛ = similar(F, Nz, Ny)
            FΣ = similar(F, Nz, Nε)
            Λ0(FΛ, z)
            Σ0(FΣ, z)
            F .= ccgf((Γ₅ + Γ₆ * Ψ) * ((I - (FΛ * Ψ)) \ FΣ), z)
        end
    else
        function _𝒱_sparse_ip(F, z, Ψ)
            FΛ = similar(F, Nz, Ny)
            FΣ = similar(F, Nz, Nε)
            Λ0(FΛ, z)
            Σ0(FΣ, z)
            ccgf(F, (Γ₅ + Γ₆ * Ψ) * ((I - (FΛ * Ψ)) \ FΣ), z)
        end
    end
    𝒱_dz = similar(z, Ny)

    J𝒱 = RALF2((F, z, Ψ) -> forwarddiff_color_jacobian!(F, (F0, x) -> _𝒱_sparse(F0, x, Ψ), z, dx = 𝒱_dz,
                                                        colorvec = colorvec, sparsity = sparsity),
               z, Ψ, deepcopy(sparsity))

    # Create mini-version of the Jacobian cache
    J𝒱_jac_cache = (dx = 𝒱_dz, sparsity = sparsity, colorvec = colorvec)

    return J𝒱, J𝒱_jac_cache
end

function construct_𝒱_jacobian_function(𝒱::RALF4, ccgf::Function, Λ::RALF2, Σ::RALF2{LC}, Γ₅::AbstractArray{<: Number},
                                       Γ₆::AbstractArray{<: Number}, z::AbstractVector{T}, y::AbstractVector{T},
                                       Ψ::AbstractMatrix{T};
                                       sparsity::Union{AbstractArray, Nothing} = nothing,
                                       colorvec = nothing, sparsity_detection::Bool = false) where {T <: Number, LC}

    # Define (temporary) objective functions
    _f_𝒱z = zₜ -> 𝒱(zₜ, y, Ψ, zₜ, (4, 2))

    # Need to grab some dimensions
    Ny, Nz = size(Ψ)
    Nε     = size(LC <: AbstractArray ? Σ.cache : Σ.cache.du, 2)

    # Infer sparsity patterns and matrix coloring vector
    if isnothing(sparsity)
        sparsity, colorvec = compute_sparsity_pattern(_f_𝒱z, z, Ny; sparsity_detection = sparsity_detection)
    elseif isnothing(colorvec)
        colorvec = matrix_colors(sparsity)
    end

    # Create RALF2 objects. Note that we cannot pre-allocate the caches for
    # forwarddiff_color_jacobian! by using ForwardColorJacCache b/c the objective function
    # changes as z and y change. If Jacobians of μ and ξ are refactored to be done once,
    # then it'll be possible to cache.
    #
    # See the previous version of construct_𝒱_jacobian_function for comments on
    # why we re-implement 𝒱 as done below.

    FΛ0 = similar(z, Nz, Ny)
    FΣ0 = similar(z, Nz, Nε)
    Λ0 = if applicable(Λ.f0, FΛ0, z, y)
        Λ.f0
    else
        function _Λ_op(F, z, y)
            F .= Λ.f0(z, y)
        end
    end
    Σ0 = if applicable(Σ.f0, FΛ0, z, y)
        Σ.f0
    else
        function _Σ_op(F, z, y)
            F .= Σ.f0(z, y)
        end
    end

    _𝒱_sparse = if applicable(ccgf, Γ₅, z)
        function _𝒱_sparse_op(F, z, y, Ψ, zₜ)
            FΛ = similar(F, Nz, Ny)
            FΣ = similar(F, Nz, Nε)
            yₜ = y + Ψ * (zₜ - z)
            Λ0(FΛ, zₜ, yₜ)
            Σ0(FΣ, zₜ, yₜ)
            F .= ccgf((Γ₅ + Γ₆ * Ψ) * ((I - (FΛ * Ψ)) \ FΣ), zₜ)
        end
    else
        function _𝒱_sparse_ip(F, z, y, Ψ, zₜ)
            FΛ = similar(F, Nz, Ny)
            FΣ = similar(F, Nz, Nε)
            yₜ = y + Ψ * (zₜ - z)
            Λ0(FΛ, zₜ, yₜ)
            Σ0(FΣ, zₜ, yₜ)
            ccgf(F, (Γ₅ + Γ₆ * Ψ) * ((I - (FΛ * Ψ)) \ FΣ), zₜ)
        end
    end

    𝒱_dz = similar(y, Ny)
    J𝒱 = RALF3((F, z, y, Ψ) -> forwarddiff_color_jacobian!(F, (F0, zₜ) -> _𝒱_sparse(F0, z, y, Ψ, zₜ), z, dx = 𝒱_dz,
                                                           colorvec = colorvec, sparsity = sparsity),
               z, y, Ψ, deepcopy(sparsity))

    # Create mini-version of the Jacobian cache
    J𝒱_jac_cache = (dx = 𝒱_dz, sparsity = sparsity, colorvec = colorvec)

    return J𝒱, J𝒱_jac_cache
end

# Helper functions for exploiting sparsity in calls to nlsolve

## Helper function for compute_sparsity_pattern
function infer_objective_function(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1)

    f = if algorithm == :deterministic
        (F, x) -> _deterministic_equations(F, x, m)
    elseif algorithm == :relaxation
        (F, x) -> _relaxation_equations(F, x, m, m.Ψ, m[:𝒱_sss])
    elseif algorithm == :homotopy
        if Λ_eltype(m.nonlinear) <: RALF1 && Σ_eltype(m.nonlinear) <: RALF1
            (F, x) -> _homotopy_equations1(F, x, m, q)
        else
            (F, x) -> _homotopy_equations2(F, x, m, q)
        end
    end

    return f
end

"""
```
compute_sparsity_pattern(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                         sparsity::Union{AbstractArray, Nothing} = nothing,
                         sparsity_detection::Bool = false)
```
calculates the sparsity pattern and matrix coloring vector of the Jacobian
of the nonlinear system of equations for either the deterministic or
stochastic steady state, depending on which `algorithm` is called.

### Keywords
- `q`: step size for homotopy. Should satisfy `0 < q < 1` and is only required to ensure
    that the sparsity pattern is correctly determined when `algorithm = :homotopy`
    and thus the dependence of the entropy `𝒱` on the coefficients `(z, y, Ψ)` matters.
- `sparsity`: sparsity pattern of the Jacobian of the nonlinear system of equations
- `sparsity_detection`: if true, use SparsityDetection.jl to determine the sparsity pattern.
    If false, then the sparsity pattern is determined by using finite differences
    to calculate a Jacobian and assuming any zeros are supposed to be zero.
"""
function compute_sparsity_pattern(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                                  sparsity::Union{AbstractArray, Nothing} = nothing,
                                  sparsity_detection::Bool = false)
    @assert algorithm in [:deterministic, :relaxation, :homotopy] "The algorithm must be :deterministic, :relaxation, or :homotopy"
    @assert 1 > q > 0 "The step size q must satisfy 0 < q < 1."

    f = infer_objective_function(m, algorithm; q = q)

    input = algorithm == :homotopy ? vcat(m.z, m.y, vec(m.Ψ)) : vcat(m.z, m.y)
    if isnothing(sparsity)
        sparsity = if sparsity_detection
            convert(SparseMatrixCSC{eltype(x), Int}, jacobian_sparsity(f, similar(input), input))
        else
            jac = similar(input, length(input), length(input))
            FiniteDiff.finite_difference_jacobian!(jac, f, input)
            sparse(jac)
        end
    end
    colorvec = matrix_colors(sparsity)

    return sparsity, colorvec
end

"""
```
preallocate_jac_cache(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                      sparsity::Union{AbstractArray, Nothing} = nothing,
                      sparsity_detection::Bool = false)
```
pre-allocates the cache for the Jacobian of the nonlinear system of equations
for either the deterministic or stochastic steady state, depending on which
`algorithm` is called.

### Keywords
- `q`: step size for homotopy. Should satisfy `0 < q < 1` and is only required to ensure
    that the sparsity pattern is correctly determined when `algorithm = :homotopy`
    and thus the dependence of the entropy `𝒱` on the coefficients `(z, y, Ψ)` matters.
- `sparsity`: the sparsity pattern of the Jacobian of the nonlinear system of equations
- `sparsity_detection`: if true, use SparsityDetection.jl to determine the sparsity pattern.
    If false, then the sparsity pattern is determined by using finite differences
    to calculate a Jacobian and assuming any zeros are supposed to be zero.
"""
function preallocate_jac_cache(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                               sparsity::Union{AbstractArray, Nothing} = nothing,
                               sparsity_detection::Bool = false)

    sparsity, colorvec = compute_sparsity_pattern(m, algorithm; q = q,
                                                  sparsity = sparsity, sparsity_detection = sparsity_detection)
    input = algorithm == :homotopy ? vcat(m.z, m.y, vec(m.Ψ)) : vcat(m.z, m.y)

    return FiniteDiff.JacobianCache(input, colorvec = colorvec, sparsity = sparsity)
end

function construct_sparse_jacobian_function(m::RiskAdjustedLinearization, f::Function,
                                            algorithm::Symbol, autodiff::Symbol;
                                            sparsity::Union{AbstractArray, Nothing} = nothing,
                                            colorvec = nothing, jac_cache = nothing,
                                            sparsity_detection::Bool = false)

    if isnothing(jac_cache)
        # Create Jacobian function that does not assume the existence of a cache

        if isnothing(sparsity) # No sparsity pattern provided, so need to make one
            sparsity, colorvec = compute_sparsity_pattern(m, algorithm; sparsity = sparsity,
                                                          sparsity_detection = sparsity_detection)
        elseif isnothing(colorvec) # Sparsity pattern, but no colorvec, so apply matrix_colors
            colorvec = matrix_colors(sparsity)
        end

        nlsolve_jacobian! = if autodiff == :forward
            (F, x) -> forwarddiff_color_jacobian!(F, f, x, # homotopy doesn't work with autodiff, so assuming
                                                  ForwardColorJacCache(f, x, min(m.Nz, m.Ny); # only using deterministic/relaxation,
                                                                       colorvec = colorvec, sparsity = sparsity)) # hence the chunk size
        else
            (F, x) -> FiniteDiff.finite_difference_jacobian!(F, f, x; colorvec = colorvec,
                                                             sparsity = sparsity)
        end

        return nlsolve_jacobian!, sparsity
    else
        # Create Jacobian function that assumes the existence of a cache

        nlsolve_jacobian! = if autodiff == :forward
            (F, x) -> forwarddiff_color_jacobian!(F, f, x, jac_cache)
        else
            (F, x) -> FiniteDiff.finite_difference_jacobian!(F, f, x, jac_cache)
        end

        return nlsolve_jacobian!, jac_cache.sparsity
    end
end


text:
module GalacticCMAEvolutionStrategy

using CMAEvolutionStrategy, GalacticOptim, GalacticOptim.SciMLBase

export CMAEvolutionStrategyOpt

struct CMAEvolutionStrategyOpt end

function __map_optimizer_args(prob::OptimizationProblem, opt::CMAEvolutionStrategyOpt;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    mapped_args = (; lower=prob.lb,
        upper=prob.ub,
        kwargs...)

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., maxiter=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., maxtime=maxtime)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., ftol=abstol)
    end

    return mapped_args
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::CMAEvolutionStrategyOpt, data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)
    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = callback(decompose_trace(trace).metadata["x"], trace.value...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    _loss = function (θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = CMAEvolutionStrategy.minimize(_loss, prob.u0, 0.1; opt_args...)
    t1 = time()

    opt_ret = opt_res.stop.reason

    SciMLBase.build_solution(prob, opt, opt_res.logger.xbest[end], opt_res.logger.fbest[end]; original=opt_res, retcode=opt_ret)
end

end


text:
using UUIDs

# Refs
#
# * https://arxiv.org/html/0901.4016
# * https://github.com/dsw/proquint/blob/master/proquint.c

_split_16(u::UInt8) = (u,) .% UInt16
_split_16(u::UInt16) = (u,) .% UInt16
_split_16(u::UInt32) = (u >> 16,
                        u & 0xffff) .% UInt16
_split_16(u::UInt64) = (u >> 48,
                        (u >> 32) & 0xffff,
                        (u >> 16) & 0xffff,
                        u & 0xffff) .% UInt16
_split_16(u::UInt128) = ((u >> 112),
                         (u >> 96) & 0xffff,
                         (u >> 80) & 0xffff,
                         (u >> 64) & 0xffff,
                         (u >> 48) & 0xffff,
                         (u >> 32) & 0xffff,
                         (u >> 16) & 0xffff,
                         u & 0xffff) .% UInt16

const consonants = ('b','d','f','g','h','j','k','l','m','n','p','r','s','t','v','z')
const vowels = ('a','i','o','u')

# Encode 16 bit segment as a "quint" of 5 characters
function encode_quint(segment)
    # From the spec
    #
    #    0 1 2 3 4 5 6 7 8 9 A B C D E F
    #   +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
    #   |con    |vo |con    |vo |con    |
    #   +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

    @inbounds begin
        (consonants[(segment >> 12) & 0xf + 1],
             vowels[(segment >> 10) & 0x3 + 1],
         consonants[(segment >> 6)  & 0xf + 1],
             vowels[(segment >> 4)  & 0x3 + 1],
         consonants[segment         & 0xf + 1])
    end
end

# Decode a "quint" of 5 characters as a 16 bit number
function decode_quint(str, index)
    u = UInt16(0)
    for i=index:index+4
        c = str[i]
        if c == 'b'      u <<= 4 ; u |= 0x0;
        elseif c == 'd'  u <<= 4 ; u |= 0x1;
        elseif c == 'f'  u <<= 4 ; u |= 0x2;
        elseif c == 'g'  u <<= 4 ; u |= 0x3;
        elseif c == 'h'  u <<= 4 ; u |= 0x4;
        elseif c == 'j'  u <<= 4 ; u |= 0x5;
        elseif c == 'k'  u <<= 4 ; u |= 0x6;
        elseif c == 'l'  u <<= 4 ; u |= 0x7;
        elseif c == 'm'  u <<= 4 ; u |= 0x8;
        elseif c == 'n'  u <<= 4 ; u |= 0x9;
        elseif c == 'p'  u <<= 4 ; u |= 0xa;
        elseif c == 'r'  u <<= 4 ; u |= 0xb;
        elseif c == 's'  u <<= 4 ; u |= 0xc;
        elseif c == 't'  u <<= 4 ; u |= 0xd;
        elseif c == 'v'  u <<= 4 ; u |= 0xe;
        elseif c == 'z'  u <<= 4 ; u |= 0xf;
        elseif c == 'a'  u <<= 2 ; u |= 0x0;
        elseif c == 'i'  u <<= 2 ; u |= 0x1;
        elseif c == 'o'  u <<= 2 ; u |= 0x2;
        elseif c == 'u'  u <<= 2 ; u |= 0x3;
        else
            throw(ArgumentError("Unexpected character '$c' in proquint"))
        end
    end
    u
end

function encode_proquint(u::Unsigned, separator::Char='-')
    segments = encode_quint.(_split_16(u))
    buf = Vector{UInt8}(undef, 6*length(segments) - 1)
    j = 1
    @inbounds for i in 1:length(segments)
        for c in segments[i]
            buf[j] = c
            j += 1
        end
        if i != length(segments)
            buf[j] = separator
            j += 1
        end
    end
    String(buf)
end

encode_proquint(u::UUID, separator::Char='-') = encode_proquint(UInt128(u), separator)

function decode_proquint(T::Type{<:Unsigned}, str)
    val = zero(T)
    i = 1
    while i <= length(str)
        val <<= 16
        val |= decode_quint(str, i)
        if i+5 <= length(str)
            c = str[i+5]
            if c != '-'
                throw(ArgumentError("Unexpected proquint separator character '$c'"))
            end
        end
        i += 6
    end
    val
end

decode_proquint(T::Type{UInt8}, str) = UInt8(decode_proquint(UInt16, str))
decode_proquint(::Type{UUID}, str) = UUID(decode_proquint(UInt128, str))


text:
################################################################################
#  Copyright 2021, Tom Van Acker                                               #
################################################################################
# StochasticPowerModels.jl                                                     #
# An extention package of PowerModels.jl for Stochastic (Optimal) Power Flow   #
# See http://github.com/timmyfaraday/StochasticPowerModels.jl                  #
################################################################################

"expected cost of active power generation"
function objective_min_expected_generation_cost(pm::AbstractPowerModel; kwargs...)
    gen_cost = Dict()

    T2 = pm.data["T2"]

    for (g, gen) in _PM.ref(pm, :gen, nw=1)
        pg = Dict(nw => _PM.var(pm, nw, :pg, g) for nw in _PM.nw_ids(pm))

        if length(gen["cost"]) == 1
            gen_cost[g] = gen["cost"][1]
        elseif length(gen["cost"]) == 2
            gen_cost[g] = gen["cost"][1]*pg[1] +
                          gen["cost"][2]
        elseif length(gen["cost"]) == 3
            gen_cost[g] = gen["cost"][1]*sum(T2.get([n-1,n-1]) * pg[n]^2 for n in _PM.nw_ids(pm)) +
                          gen["cost"][2]*pg[1] +
                          gen["cost"][3]
        else
            gen_cost[g] = 0.0
        end
    end

    return JuMP.@objective(pm.model, Min,
            sum(gen_cost[g] for g in _PM.ids(pm, :gen, nw=1))
    )
end


text:
function checkInputs()
#Check various inputs

 println("Let's check inputs...!")

# Check that the hard coded income and PDF have the correct dimension
if (isUncertainty == 1) && (uncertaintyMethod == 0)
    if (size(hcIncome, 1)!= numPointsY) || (size(hcIncPDF, 1)!= numPointsY)
        error("hcIncome and must have dimension numPointsY * 1 and numPointsY * 1 respectively")
    end
end

#
# Check if uncertainty is off that numPointsY is set equal to 1
if (numPointsY!=1) && (isUncertainty == 0)
    error("There is no uncertainty but numPointsY was not set to 1.")
end

#
# Check that the value of rho is not greater than 1. Warn if close to 1
if (rho>0.999) || (rho <-0.999)
    if (rho>1)  || (rho < 1)
    error("rho is greater than 1. This code solves a stationary income process. rho greater than 1 implies a non-stationary process")
    else
    warning("rho is greater than 0.99. This code solves a stationary income process. As rho gets closer to 1 the process becomes nonstationary - possibility of numerical instability")
    end
end

#
# Check that the standard deviation is not too small
if (sigma<1e-10) && (isUncertainty == 1)
    if (sigma<=0)
        error("sigma is less than or equal to zero")
    else
        warning("sigma is very small and close to zero - possibility of numerical instability. Consider turning uncertainty off")
    end
end

# Check a number of inputs that need to be either 0 or 1
if (linearise != 0) && (linearise !=1)
    error("linearise should be either 0 or 1")
end

if (uncertaintyMethod != 0) && (uncertaintyMethod !=1)
    error("uncertaintyMethod should be either 0 or 1")
end

if (borrowingAllowed != 0) && (borrowingAllowed !=1)
    error("borrowingAllowed should be either 0 or 1")
end

if (isUncertainty != 0) && (isUncertainty !=1)
    error("isUncertainty should be either 0 or 1")
end

println("Inputs are fine :)")

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
module NeXLCore

using Requires
using Reexport
using LinearAlgebra
using DataFrames
using DataDeps

@reexport using PeriodicTable
@reexport using NeXLUncertainties
@reexport using PhysicalConstants.CODATA2018

# Abstract model types
abstract type NeXLAlgorithm end
export NeXLAlgorithm

include("ffast.jl") # Algorithms implemented in FFAST.jl
include("dtsa_mac.jl") # Implement's Heinrich's IXCOM 11 MACs
include("botesalvat.jl") # Algorithms implemented in BoteSalvat.jl
include("cullen_eadl.jl") # Data from Cullen's Evaluated Atomic Data Library
include("element.jl") # Element data from PeriodicTable.jl
include("shell.jl") # Atomic shell methods
include("transition.jl") # Atomic transition methods
include("characteristic.jl") # Characteristic X-ray methods
include("parse.jl") # Parse elements, shells, transitions from strings
include("algorithm.jl") # Default implementation of algorithms

export element # Construct Element structs
export eachelement # Elements for which there is a complete set of data :-(
export Shell # K, L, M, N etc
export KShell, LShell, MShell, NShell # The first few Shell(s)
export SubShell # K, L1, L2, L3, M1, ...
export K, L1, L2, L3, M1, M2, M3, M4, M5, N1, N2, N3, N4, N5, N6, N7 # SubShell(s)
export n, l, j # Quantum numbers
export allsubshells, ksubshells, lsubshells, msubshells, nsubshells, osubshells, pshell # SubShell type lists
export subshell # Construct SubShell structs from a string
export firstsubshell, lastsubshell # Given a shell find the lowest/highest subshell in that shell i.e. Shell[M]=>M1/M5
export AtomicSubShell # A SubShell plus an Element
export atomicsubshell # Construct AtomicSubShell structs from a string
export capacity # The total shell capacity
export jumpratio # The jump ratio for the specified shell
export meanfluorescenceyield # The mean shell-based fluorescence yield
export configuration # A string containing the electronic configuration for an Element
export Transition # An X-ray transition
export transition # Constructs Transition from SubShell objects or a string
export alltransitions, ktransitions, ltransitions, mtransitions, ntransitions, otransitions
export kalpha, kbeta, kother # K-L?, K->M? and K->!L
export lalpha, lbeta, malpha, mbeta
export CharXRay # A characteristic X-ray
export characteristic # Constructs CharXRay
export inner, outer  # Returns AtomiShell for inner and outer CharXRay
export a  # Atomic weight
export z # Atomic number
export atomic_weight # Atomic weights with uncertainties / intervals according to https://ciaaw.org/atomic-weights.htm
export symbol # Atomic symbol ("H", "He" etc)
export name # Full English name
export energy # Returns CharXRay and AtomicSubShell eneries
export density # Returns Element or Mateial data
export λ, ν, ω, wavenumber # wavelength, frequency, angular frequency and wavenumber of X-ray
export edgeenergy # Ionization edge energy for an X-ray
export NormalizeBySubShell, NormalizeByShell, NormalizeToUnity, RawYield
export weight # Returns CharXRay weights as scaled by NormalizeBySubShell, NormalizeByShell, NormalizeToUnity
export has # Element has a specific Transition
export FFASTDB # Chantler's FFAST database
export DTSA   # Heinrich's IXCOM 11 MACs
export mac # Calculates the MAC using the default or a specified algorithm
export macU # Calculates the MAC using the default or a specified algorithm
export set_user_mac! # Specify a custom user mac
export delete_user_mac! # Delete a specific user mac
export clear_user_macs! # Clear all user macs
export shell # The shell (Shell(K),Shell(L),Shell(M),...) for an AtomicSubShell, Transition, CharXRay etc.
export transitionsbyshell # Dictionary mapping transition Shell to lists of Transition(s)
export transitionsbysubshell # Dictionary mapping transition SubShell to lists of Transition(s)
export atomicsubshells # Gets an iterator of AtomicSubShell for the specified element
export brightest # Returns the brightest characteristic X-ray from a set of transitions for an element
export splitbyshell # Creates a Dict{AtomicSubShell,CharXRay} from a collection of CharXRay and the associated inner AtomicSubShell.
export relativeionizationcrosssection # Computes a number proportional to the ionization crosssection
export ionizationcrosssection # Computes the absolute ionization crosssection
export exists # Does a transition occur (according to our list...)
export @n_str # Parses a string into an Element, SubShell, AtomicSubShell, Transition or CharXRay
export @enx_str # Energy of a atomic sub-shell or characteristic X-ray in string form
export Bote2009

include("compton.jl")
export comptonAngular # Computes the angular distribution of Compton
export comptonShift # Computes the fractional Compton shift
export comptonEnergy # Computes the resulting Compton-shifted X-ray energy
export comptonDifferential # Computes the differential Compton cross-section

include("material.jl")
include("matparser.jl")
export material # Construct a Material struct
export pure # Construct a pure element material
export massfraction # Returns the composition as mass fraction
export atomicfraction # Returns the composition as atom fraction
export normalizedmassfraction # Returns the composition as normalized mass fraction
export normalized # A single element normalized
export asnormalized # Returns a normalized copy of a Material
export analyticaltotal # Returns the analytical total
export labeled # Transform a data item into a Dict of (Label, value)
export Material # Material struct
# export Base.keys # Element keys into Material
export name # Material name
export compare # Compare compositions as a DataFrame
export elms # Use elms instead of elements since elements taken by PeriodicTable
export nonneg # Returns the mass fraction as a Float64 >= 0.0
export atoms_per_cm³ # Returns the number of atoms in one cm² of the Material
export atoms_per_g # Returns the number of atoms in one gram of the Material
export @mat_str
export delete # Create a new Material with certain elements removed

include("film.jl")
export Film # A thin film of Material
export transmission # Transmission fraction through a Film
export thickness # Of Film in cm
export massthickness # Of film in g/cm²

# For DTSA-II interop
export todtsa2comp    # Write a mass fraction (+opt density) to a string parsable by parsedtsa2comp
export parsedtsa2comp # Parse a mass fraction (+opt density) as written by todtsa2comp

include("stoichiometry.jl")
export asoxide # Compute the standard oxide
export valence # a table of elemental valences (valence[z(n"O")] = -2)
export obystoichiometry # computes the mass fraction of O using valence rules.

export NeXLPalette
export plotXrayEnergies # A Gadfly plot of X-ray energies for a set of transitions
export plotXrayWeights # Plot weights of lines
export plotEdgeEnergies # Plot edge energies

include("bremsstrahlung.jl")
export NeXLBremsstrahlung # <: NeXLAlgorithm
export Kramers1923 # <: NeXLBremsstrahlung
export Lifshin1974 # <: NeXLBremsstrahlung
export Reed1975 # <: NeXLBremsstrahlung
export Smith1975 # <: NeXLBremsstrahlung
export Small1987 # <: NeXLBremsstrahlung
export Trincavelli1997 # <: NeXLBremsstrahlung
export Castellano2004a # <: NeXLBremsstrahlung
export Castellano2004b # <: NeXLBremsstrahlung
export bremsstrahlung # Bremsstrahlung model

include("meanionizationpotential.jl")
export NeXLMeanIonizationPotential
export Bloch1933 # <: NeXLMeanIonizationPotential
export Jensen1937 # <: NeXLMeanIonizationPotential
export Wilson1941 # <: NeXLMeanIonizationPotential
export Sternheimer1964 # <: NeXLMeanIonizationPotential
export Springer1967 # <: NeXLMeanIonizationPotential
export Zeller1973 # <: NeXLMeanIonizationPotential
export Brizuela1990 # <: NeXLMeanIonizationPotential
export Berger1982 # <: NeXLMeanIonizationPotential
export J # Mean ionization potential

include("eta.jl")
export LoveScott1978η
export Tomlin1963
export August1989η
export Pouchou1991η
export η
export zbar  # Calculated using Donovan /

include("bethe.jl")
export BetheEnergyLoss # <: BetheEnergyLoss
export Bethe # <: BetheEnergyLoss
export JoyLuo # <: BetheEnergyLoss
export dEds
export Kanaya1972

include("electron.jl")
export λₑ # Wavelength of an electron
export kₑ # Wave number of an electron
export mₑ # mass of an electron in MeV

include("misc.jl")
export Poehn1985 # A jump ratio model
export jumpratio # jump ratio algorithm
export klinewidths # K shell linewidths from Bambynek'1974 errata to Bambynek 1972
export Burhop1965, Sogut2002, Krause1979, Kahoul2012, Reed1975ω # Fluorescence yield models

export CullenEADL
export fluorescenceyield # fluorescence yield models

include("matu.jl")
export MaterialLabel
export MassFractionLabel
export NormMassFractionLabel
export AtomicFractionLabel
export AtomicWeightLabel
export MaterialFractionLabel
export MatStatTypes
export MeanAz, MeanZ, AnalyticalTotal
export μoρElementLabel
export μoρLabel
export μoρMaterial # MeasurementModel
export MeanZ # MeasurementModel
export MeanAz # MeasurementModel
export MFtoAF # MeasurementModel
export MFtoNMF # MeasurementModel
export AFtoMF # MeasurementModel
export MatStats # MeasurementModel
export mf2comp
export materiallabels

include("scattering.jl")
export ElasticScatteringCrossSection # abstract
export ScreenedRutherfordType # abstract
export ScreenedRutherford
export Liljequist1989
export Browning1991
export Browning1994
export λ, λₜᵣ  # Mean free path and transport mean free path
export σₜ # total scattering cross section
export δσδΩ # The differential cross section
export Rₐ # Classical atomic scattering radius

include("kratio.jl")
export KRatio # Represents a measured intensity ratio
export KRatios # The array equivalent of KRatio
export nonnegk # Returns the k-ratio value truncated to non-negative.
export elms  # Returns a list of the elements in a `List{KRatio}`
export colorize # Colorize a X-ray
export normalizek # Normalize a Vector{KRatio[s]}
#export strip # Removes one or more elements from a Vector{KRatio}

include("properties.jl")
# Methods for checking what properties are required for an algorithm.
export minproperties # A list of the minimum required properties
export hasminrequired # Checks whether a spectrum has necessary properties
export requiredbutmissing # Lists missing properties

include("custommac.jl")
export CustomMAC  # Tied to "data\specialmacs.csv"
export getcustommacs # Retrieve a set of custom MAC values (:Henke1974, :Henke1982, :Bastin19XX, etc. (see specialmacs.csv))

include("materials.jl")
export loadsmithsoniandata, parsedsmithsoniandata
export loadmineraldata
export srm470_k411, srm470_k412

include("palettes.jl")
export Log3BandC, Log3Band, LogScale, LinearScale # Converts [0.0,1.0] to a colored image.
export loadlegend # Load a legend for one of the above scales

include("mc.jl")
include("mchelpers.jl")
export Position, previous # Base.position
export Particle, Electron
export RectangularShape, SphericalShape # aliases for 3D GeometryBasics shapes
export Region # A shape and a material
export scatter # The Particle transport function
export gun # A source of starter Particle objects
export trajectory # Calculates Point and Region pairs as the Particle traverses the sample
export intersection # Compute how far along a ray, the ray intersects a shape.
export chamber, particle, bulk, thin_film, coated_particle
export colorize # Maps Material to Color for all Material in a Region

include("staging.jl")
export StageMapping, DefaultStageMapping
export stage2image, image2stage
export compute_tilt # Computes the tilt and orientation from three focus points

include("standardize.jl")
export isstandard # Does a k-ratio have the necessary properties to be a standard
export standardize # Apply similar standards to a KRatio or KRatios

function __init__()
    @require Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004" include("gadflyplot.jl")
    @require MeshCat = "283c5d60-a78f-5afe-a0af-af636b173e11" include("meshcat.jl")
    @require SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" include("materialdb.jl")

    register(
        DataDep(
            "RUFFDatabase",
            """
            Dataset: RUFF Mineral Database
            Website: https://rruff.info/ima/
            License: Creative Commons Attribution-ShareAlike 3.0 Unported License (http://creativecommons.org/licenses/by-sa/3.0/)
            """,
            "https://drive.google.com/uc?export=download&id=1Ackbz0YtaliQNCdZmPfj9uWTwhVNBzy8",
            "c87ce711c6fa9d8a012e29f949c075e176ea9a7c753e9e79270e4e9fd988e613",
            post_fetch_method = unpack
        )
    )
end

end


text:
abstract type MeanFunction end

"""
    ZeroMean{T<:Real} <: MeanFunction

Returns `zero(T)` everywhere.
"""
struct ZeroMean{T<:Real} <: MeanFunction end

_map(::ZeroMean{T}, x::AbstractVector) where {T} = zeros(T, length(x))

function ChainRulesCore.rrule(::typeof(_map), m::ZeroMean, x::AbstractVector)
    map_ZeroMean_pullback(Δ) = (NoTangent(), NoTangent(), ZeroTangent())
    return _map(m, x), map_ZeroMean_pullback
end

ZeroMean() = ZeroMean{Float64}()

"""
    ConstMean{T<:Real} <: MeanFunction

Returns `c` everywhere.
"""
struct ConstMean{T<:Real} <: MeanFunction
    c::T
end

_map(m::ConstMean, x::AbstractVector) = fill(m.c, length(x))

"""
    CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy. Must be able to be mapped over an
`AbstractVector` of inputs.
"""
struct CustomMean{Tf} <: MeanFunction
    f::Tf
end

_map(f::CustomMean, x::AbstractVector) = map(f.f, x)


text:
module Photodynamics

using Reexport
@reexport using NbodyGradient

using StaticArrays, LinearAlgebra, Roots, SpecialFunctions, Limbdark

import Limbdark: transit_poly_g, transit_poly_g!, Transit_Struct
import NbodyGradient: InitialConditions, Derivatives, ElementsIC
import NbodyGradient: check_step, set_state!

export dot, compute_lightcurve!
export TransitSeries, Lightcurve

include("simpson.jl")
include("impact.jl")
include("TransitSeries.jl")
include("Lightcurve.jl")

"""Constant coefficients for the series expansion."""
struct Coefficients{T<:Real}
    vd::StaticVector{3,T}
    ad::StaticVector{4,T}
    jd::StaticVector{3,T}
    sd::StaticVector{4,T}

    function Coefficients(::Type{T}=Float64) where T<:Real
        vd = SVector{3,T}(1.0/60, 9/60, 45/60)
        ad = SVector{4,T}(1.0/90, -3/20, 3/2, -49/18)
        jd = SVector{3,T}(1.0/8, 1, 13/8)
        sd = SVector{4,T}(-1.0/6, 2, -13/2, 28/3)
        return new{T}(vd,ad,jd,sd)
    end
end

const COEFF = Coefficients();

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ._ProtoBuf_Top_.opentelemetry

const Status_StatusCode = (;
    [
        Symbol("STATUS_CODE_UNSET") => Int32(0),
        Symbol("STATUS_CODE_OK") => Int32(1),
        Symbol("STATUS_CODE_ERROR") => Int32(2),
    ]...
)

mutable struct Status <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Status(; kwargs...)
        obj = new(meta(Status), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Status
const __meta_Status = Ref{ProtoMeta}()
function meta(::Type{Status})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Status)
            __meta_Status[] = target = ProtoMeta(Status)
            fnum = Int[2, 3]
            allflds =
                Pair{Symbol,Union{Type,String}}[:message=>AbstractString, :code=>Int32]
            meta(
                target,
                Status,
                allflds,
                ProtoBuf.DEF_REQ,
                fnum,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_Status[]
    end
end
function Base.getproperty(obj::Status, name::Symbol)
    if name === :message
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :code
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

const Span_SpanKind = (;
    [
        Symbol("SPAN_KIND_UNSPECIFIED") => Int32(0),
        Symbol("SPAN_KIND_INTERNAL") => Int32(1),
        Symbol("SPAN_KIND_SERVER") => Int32(2),
        Symbol("SPAN_KIND_CLIENT") => Int32(3),
        Symbol("SPAN_KIND_PRODUCER") => Int32(4),
        Symbol("SPAN_KIND_CONSUMER") => Int32(5),
    ]...
)

mutable struct Span_Event <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Span_Event(; kwargs...)
        obj = new(meta(Span_Event), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Span_Event
const __meta_Span_Event = Ref{ProtoMeta}()
function meta(::Type{Span_Event})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Span_Event)
            __meta_Span_Event[] = target = ProtoMeta(Span_Event)
            wtype = Dict(:time_unix_nano => :fixed64)
            allflds = Pair{Symbol,Union{Type,String}}[
                :time_unix_nano=>UInt64,
                :name=>AbstractString,
                :attributes=>Base.Vector{opentelemetry.proto.common.v1.KeyValue},
                :dropped_attributes_count=>UInt32,
            ]
            meta(
                target,
                Span_Event,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                wtype,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_Span_Event[]
    end
end
function Base.getproperty(obj::Span_Event, name::Symbol)
    if name === :time_unix_nano
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.common.v1.KeyValue,
        }
    elseif name === :dropped_attributes_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    else
        getfield(obj, name)
    end
end

mutable struct Span_Link <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Span_Link(; kwargs...)
        obj = new(meta(Span_Link), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Span_Link
const __meta_Span_Link = Ref{ProtoMeta}()
function meta(::Type{Span_Link})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Span_Link)
            __meta_Span_Link[] = target = ProtoMeta(Span_Link)
            allflds = Pair{Symbol,Union{Type,String}}[
                :trace_id=>Vector{UInt8},
                :span_id=>Vector{UInt8},
                :trace_state=>AbstractString,
                :attributes=>Base.Vector{opentelemetry.proto.common.v1.KeyValue},
                :dropped_attributes_count=>UInt32,
            ]
            meta(
                target,
                Span_Link,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_Span_Link[]
    end
end
function Base.getproperty(obj::Span_Link, name::Symbol)
    if name === :trace_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :span_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :trace_state
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.common.v1.KeyValue,
        }
    elseif name === :dropped_attributes_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    else
        getfield(obj, name)
    end
end

mutable struct Span <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Span(; kwargs...)
        obj = new(meta(Span), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Span
const __meta_Span = Ref{ProtoMeta}()
function meta(::Type{Span})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Span)
            __meta_Span[] = target = ProtoMeta(Span)
            wtype = Dict(:start_time_unix_nano => :fixed64, :end_time_unix_nano => :fixed64)
            allflds = Pair{Symbol,Union{Type,String}}[
                :trace_id=>Vector{UInt8},
                :span_id=>Vector{UInt8},
                :trace_state=>AbstractString,
                :parent_span_id=>Vector{UInt8},
                :name=>AbstractString,
                :kind=>Int32,
                :start_time_unix_nano=>UInt64,
                :end_time_unix_nano=>UInt64,
                :attributes=>Base.Vector{opentelemetry.proto.common.v1.KeyValue},
                :dropped_attributes_count=>UInt32,
                :events=>Base.Vector{Span_Event},
                :dropped_events_count=>UInt32,
                :links=>Base.Vector{Span_Link},
                :dropped_links_count=>UInt32,
                :status=>Status,
            ]
            meta(
                target,
                Span,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                wtype,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_Span[]
    end
end
function Base.getproperty(obj::Span, name::Symbol)
    if name === :trace_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :span_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :trace_state
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :parent_span_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :kind
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :start_time_unix_nano
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :end_time_unix_nano
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.common.v1.KeyValue,
        }
    elseif name === :dropped_attributes_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :events
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Span_Event}
    elseif name === :dropped_events_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :links
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Span_Link}
    elseif name === :dropped_links_count
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :status
        return (obj.__protobuf_jl_internal_values[name])::Status
    else
        getfield(obj, name)
    end
end

mutable struct InstrumentationLibrarySpans <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function InstrumentationLibrarySpans(; kwargs...)
        obj = new(meta(InstrumentationLibrarySpans), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct InstrumentationLibrarySpans
const __meta_InstrumentationLibrarySpans = Ref{ProtoMeta}()
function meta(::Type{InstrumentationLibrarySpans})
    ProtoBuf.metalock() do
        if !isassigned(__meta_InstrumentationLibrarySpans)
            __meta_InstrumentationLibrarySpans[] =
                target = ProtoMeta(InstrumentationLibrarySpans)
            allflds = Pair{Symbol,Union{Type,String}}[
                :instrumentation_library=>opentelemetry.proto.common.v1.InstrumentationLibrary,
                :spans=>Base.Vector{Span},
                :schema_url=>AbstractString,
            ]
            meta(
                target,
                InstrumentationLibrarySpans,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_InstrumentationLibrarySpans[]
    end
end
function Base.getproperty(obj::InstrumentationLibrarySpans, name::Symbol)
    if name === :instrumentation_library
        return (obj.__protobuf_jl_internal_values[name])::opentelemetry.proto.common.v1.InstrumentationLibrary
    elseif name === :spans
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Span}
    elseif name === :schema_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct ResourceSpans <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ResourceSpans(; kwargs...)
        obj = new(meta(ResourceSpans), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ResourceSpans
const __meta_ResourceSpans = Ref{ProtoMeta}()
function meta(::Type{ResourceSpans})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ResourceSpans)
            __meta_ResourceSpans[] = target = ProtoMeta(ResourceSpans)
            allflds = Pair{Symbol,Union{Type,String}}[
                :resource=>opentelemetry.proto.resource.v1.Resource,
                :instrumentation_library_spans=>Base.Vector{InstrumentationLibrarySpans},
                :schema_url=>AbstractString,
            ]
            meta(
                target,
                ResourceSpans,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ResourceSpans[]
    end
end
function Base.getproperty(obj::ResourceSpans, name::Symbol)
    if name === :resource
        return (obj.__protobuf_jl_internal_values[name])::opentelemetry.proto.resource.v1.Resource
    elseif name === :instrumentation_library_spans
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            InstrumentationLibrarySpans,
        }
    elseif name === :schema_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TracesData <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TracesData(; kwargs...)
        obj = new(meta(TracesData), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct TracesData
const __meta_TracesData = Ref{ProtoMeta}()
function meta(::Type{TracesData})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TracesData)
            __meta_TracesData[] = target = ProtoMeta(TracesData)
            allflds =
                Pair{Symbol,Union{Type,String}}[:resource_spans=>Base.Vector{ResourceSpans}]
            meta(
                target,
                TracesData,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_TracesData[]
    end
end
function Base.getproperty(obj::TracesData, name::Symbol)
    if name === :resource_spans
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ResourceSpans}
    else
        getfield(obj, name)
    end
end

export TracesData,
    ResourceSpans,
    InstrumentationLibrarySpans,
    Span_SpanKind,
    Span_Event,
    Span_Link,
    Span,
    Status_StatusCode,
    Status


text:
"""Module for the modeling of the dispersion of pollutants in aquifers."""
module Aquifers

using ArchGDAL
using Dates
using SpecialFunctions

include(".\\Utils\\Functions.jl")
include(".\\Utils\\FunctionsDB.jl")

export run_leaching

const agd = ArchGDAL

mutable struct Leach
 # Parameters
    henry_const::Float64                # Henry's constant
    volumetric_water_content::Float64   # Volumetric water content (Θw)
    volumetric_air_content::Float64     # Volumetric air content   (Θa)
    soil_adsorption::Float64            # Soil adsorption coefficient (kw)
    effective_infiltration::Float64     # Effective infiltration (ief)
    soil_density::Float64               # Soil density (ρs)
    source_thickness::Float64           # Source thickness
    aquifer_depth::Float64              # Aquifer depth
    darcy_velocity::Float64             # Darcy velocity
    mixing_zone_depth::Float64          # Depth of mixing zone
    orthogonal_width::Float64           # Source extend
 # Computational results
    kw::Float64                         # Soil-Water partition coefficient
    koc::Float64                        # Organic carbon partition co-efficient
    ldf::Float64                        # Diluition factor
    sam::Float64                        # Soil attenuation coefficient
    leaching_factor::Float64            # Leaching factor

    Leach(henry_const,volumetric_water_content,volumetric_air_content,soil_adsorption,effective_infiltration,soil_density,source_thickness,aquifer_depth,darcy_velocity,mixing_zone_depth,orthogonal_width) = new(henry_const,volumetric_water_content,volumetric_air_content,soil_adsorption,effective_infiltration,soil_density,source_thickness,aquifer_depth,darcy_velocity,mixing_zone_depth,orthogonal_width,0.0,0.0,0.0,0.0,0.0)
end

mutable struct DAF <: Functions.AbstractAnalysisObject
 # Parameters
    secondary_source_concentration::Float64   # Concentration of pollutant at the secondary source
    x::Float64                                # X coordinate of target point
    y::Float64                                # Y coordinate of target point
    α_x::Float64                              #
    α_y::Float64                              #
    decay_coeff::Float64                      # Decay coefficient
    darcy_velocity::Float64                   # Darcy velocity
    soil_adsorption::Float64                  # Soil adsorption coefficient
    soil_density::Float64                     # Soil density
    tera_e::Float64                           #
    orthogonal_width::Float64                 #
    time::Int64                               # Time for the analysis
    direction::Int64                          # Direction of the flow
    algorithm::Symbol                         # Type of alorithm to use
    option::Symbol                            # Version of the algorithm to use
 # Computational results
    R::Float64
    DAF::Float64
    DAF_tot::Float64

    DAF(secondary_source_concentration,x,y,α_x,α_y,decay_coeff,darcy_velocity,soil_adsorption,soil_density,tera_e,orthogonal_width,time,direction,algorithm,option) = new(secondary_source_concentration,x,y,α_x,α_y,decay_coeff,darcy_velocity,soil_adsorption,soil_density,tera_e,orthogonal_width,time,direction,algorithm,option)
end

# =============================================== Leaching functions ================================================================================

function calc_kw!( l::Leach )
    l.koc = l.soil_adsorption
    l.soil_adsorption = 0.01l.koc
    l.kw = l.soil_density / ( l.volumetric_water_content + ( l.soil_adsorption * l.soil_density ) + ( l.henry_const * l.volumetric_air_content ) )
    return l.kw
end

function calc_ldf!( l::Leach )
    darcy = l.darcy_velocity * 100.0 * 86400.0 * 365.0
    l.ldf = 1.0 + ( darcy * ( l.mixing_zone_depth / ( l.effective_infiltration * l.orthogonal_width ) ) )
    return l.ldf
end

function calc_sam!( l::Leach )
    l.sam = l.source_thickness / l.aquifer_depth
    return l.sam
end

function calc_LF!( l::Leach )
    l.leaching_factor = ( l.kw * l.sam ) / l.ldf
    return l.leaching_factor
end

# ================================================ DAF functions ====================================================================================

function calc_R!( c::DAF )
    c.R = 1.0 + ( c.soil_adsorption * ( c.soil_density / c.tera_e ) )
    return c.R
end

function calc_DAF_ispra!( c::DAF )
    ######################## modello di domenico ###########################
    # vedere appendice C pagina 2 del documento Criteri metodologici per l'applicazione dell'analisi assoluta di rischio ai siti contaminati
    # la formula originale prevede la produttoria delle 3 componenti x,y,z moltiplicata per 1/4
    # eliminando la terza componente dell'asse z è necessario moltplicare per 1/2 (quindi 0.5)
    # per verifica vedere Domenico P.A. e Schwartz F.W. (1998), Physical and Chemical Hydrogeology, John Wiley and Sons, New York.
    # da pagina 642 a pag 644
    if c.α_x == 0.0
      c.α_x = 0.1c.x
    end
    if c.α_y == 0.0
      c.α_y = c.α_x / 3.0
    end
    R = 1.0 + ( c.soil_adsorption * ( c.soil_density / c.tera_e ) )
    daf1 = 0.50ℯ^( ( c.x / 2.0c.α_x ) * ( 1 - √( 1.0 + ( ( 4.0c.decay_coeff * c.α_x * R ) / c.darcy_velocity ) ) ) )
    daf21 = erf( ( c.y + 0.5c.orthogonal_width ) / ( 2.0√( c.α_y * c.x ) ) )
    daf22 = erf( ( c.y - 0.5c.orthogonal_width ) / ( 2.0√( c.α_y * c.x ) ) )
    daf3 = daf21 - daf22
    DAF_tot = daf1 * daf3
    return DAF_tot
end

function calc_DAF_ispra2!( c::DAF )
    if c.α_x == 0.0
      c.α_x = 0.1c.x
    end
    if c.α_y == 0.0
      c.α_y = c.α_x / 3.0
    end
    #daf1 = ( c.x / 2c.α_x ) * ( 1 - √( 1 + ( ( 4c.decay_coeff * c.α_x * c.R ) / c.darcy_velocity ) ) )
    daf1 = exp( c.x / ( 2.0c.α_x ) * 0 )
    #daf1e = exp(daf1)
    daf2 = erf( c.orthogonal_width / ( 4.0√( c.α_y * c.x ) ) )
    c.DAF = daf1 * daf2
    return c.DAF
end

function calc_DAF!( c::DAF )
    if c.α_x == 0
      c.α_x = 0.1( c.x / 100.0 )
    end
    if c.α_y == 0
      c.α_y = c.α_x / 3.0
    end
    dx = c.α_x * c.darcy_velocity
    dy = c.α_y * c.darcy_velocity
    daf_a = c.secondary_source_concentration / ( 4.0 * c.tera_e * π * c.time * √(dx * dy) )
    daf_b = ℯ^( -( ( (( c.x - (c.darcy_velocity * c.time) )^2.0) / (4.0 * dx * c.time) ) + ((c.y^2.0) / (4.0 * dy *c.time)) ) )
    c.DAF = daf_a * daf_b
    return c.DAF
end

function calc_DAF_uni!( c::DAF )
    if c.α_x == 0
      c.α_x = 0.1c.x
    end
    dx = c.α_x * c.darcy_velocity
    daf_a = c.secondary_source_concentration / ( 2.0 * c.tera_e * √( 4.0 * dx * π * c.time ) )
    daf_b = ℯ^( -(( ((c.x - (c.darcy_velocity * c.time))^2.0) / (4.0 * dx * c.time) )) )
    c.DAF = daf_a * daf_b
    return c.DAF
end

function calc_DAF_c!( c::DAF )
    #continuous
    if c.α_x == 0
      c.α_x = 0.1c.x
    end
    if c.α_y == 0
      c.α_y = c.α_x / 3.0
    end
    dx = c.α_x * c.darcy_velocity
    dy = c.α_x * c.darcy_velocity
    r = √( (c.x^2.0) + ( (c.y^2.0) * (dx / dy) ) )
    daf_a = c.secondary_source_concentration / ( 4.0 * c.tera_e * √(π) * √(c.darcy_velocity * r) * √(dy) )
    daf_b = ℯ^( ( (c.x - r) * c.darcy_velocity ) / (2.0dx) )
    c.DAF = daf_a * daf_b
    return c.DAF
end

function Functions.compute_concentration!( d::DAF )
    concentration = 0.0
    if d.x > 0
     #=
        if d.algorithm == :fickian
            if d.option == :pulse
                concentration = calc_DAF!(d)
            else
                concentration = calc_DAF_c!(d)
            end
        else
            concentration = d.secondary_source_concentration * calc_DAF_ispra!(d)
        end
     =#
        concentration = d.algorithm == :fickian ? d.option == :pulse ? calc_DAF!(d) : calc_DAF_c!(d) : d.secondary_source_concentration * calc_DAF_ispra!(d)
    end
    return concentration
end

"""
    run_leaching(; dem_file::String, source_file::String, area_file::String="", contaminantCASNum::String, concentration::Float64, tollerance::Int64=2,
                   aquifer_depth::Float64, aquifer_flow_direction::Int64, mean_rainfall::Float64, texture::String, resolution::Int64, time::Int64=1,
                   orthogonal_width::Float64=10000.0, soil_density::Float64=1.70, source_thickness::Float64=1.0, darcy_velocity::Float64=0.000025,
                   mixing_zone_depth::Float64=1.0, decay_coeff::Float64=0.0, algorithm::Symbol=:fickian, option::Symbol=:continuous,
                   output_path::String=".\\aquifer_output_model.tiff" )

Run the simulation of leaching and dispersion of contaminants in an aquifer, returning a map of the possible worst case spreading of the contaminants.

# Arguments
- `dtm_file::String`: path to the raster of terrain.
- `source_file::String`: path to the shapefile containing source point of the contaminants.
- `area_file::String=""`: path to the shapefile containing the poligon delimiting the area for the analysis.
- `contaminantCASNum::String`: CAS number identifier of a substance.
- `concentration::Float64`: concentration of the contaminants at the source.
- `tollerance::Int64=2`: value used to determine wether the concentration of pollutant in a cell is relevant.
    Specifically, a concentration value is considered relevant if its value is within "tollerance" orders of magnitute from the concentration on other cells.
- `aquifer_depth::Float64`: depth of the aquifer in meters.
- `aquifer_flow_direction::Int64`: direction of the water flow within the aquifer as an angle in degrees.
- `mean_rainfall::Float64`: average rainfall volume.
- `texture::String`: type of terrain at the source.
- `resolution::Float64`: dimension of a cell for the analysis.
- `time::Int64=1`: starting time.
- `orthogonal_width::Float64=10000.0`
- `soil_density::Float64=1.70`: density of the terrain.
- `source_thickness::Float64::Float64=1.0`: thickness of the terrain layer at the source.
- `darcy_velocity::Float64=0.000025`
- `mixing_zone_depth::Float64=1.0`
- `decay_coeff::Float64=0.0`
- `algorithm::Symbol=:fickian`: type of algorithm to be used.
- `option::Symbol=:continuous`: second option to define the kind o algorithm to use.
- `output_path::String=".\\aquifer_output_model.tiff": output file path.
"""
function run_leaching(; dem_file::String, source_file::String, area_file::String="", contaminantCASNum::String, concentration::Float64, tollerance::Int64=2,
                        aquifer_depth::Float64, aquifer_flow_direction::Int64, mean_rainfall::Float64, texture::String, resolution::Float64, time::Int64=1,
                        orthogonal_width::Float64=10000.0, soil_density::Float64=1.70, source_thickness::Float64=1.0, darcy_velocity::Float64=0.000025,
                        mixing_zone_depth::Float64=1.0, decay_coeff::Float64=0.0, algorithm::Symbol=:fickian, option::Symbol=:continuous,
                        output_path::String=".\\aquifer_output_model.tiff" )

    if algorithm ∉ [:fickian, :domenico]
        throw(DomainError(algorithm, "`algorithm` must either be `:fickian` or `:domenico`"))
    end

    if option ∉ [:pulse, :continuous]
        throw(DomainError(option, "`option` must either be `:continuous` or `:pulse`"))
    end

    henry_const, soil_adsorption = FunctionsDB.substance_extract(contaminantCASNum, ["c_henry", "koc_kd"])[1, :]

    if isempty(henry_const) && isempty(soil_adsorption)
        throw(DomainError(contaminantCASNum, "Analysis error, check input parameters"))
    end

    volumetric_air_content, volumetric_water_content, effective_infiltration, tera_e, grain = FunctionsDB.texture_extract(texture, ["tot_por", "c_water_avg", "ief", "por_eff", "grain"])[1, :]

    if any(isempty.([volumetric_air_content, volumetric_water_content, effective_infiltration, tera_e, grain]))
        throw(DomainError(texture, "Analysis error, check input parameters"))
    end

    src_geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(src_geom) != 0
        throw(DomainError(source, "`source` must be a point"))
    end

    dem = agd.read(dem_file)
    refsys = agd.getproj(dem)

    if refsys != agd.toWKT(agd.getspatialref(src_geom))
        throw(DomainError("The reference systems are not uniform."))
    end

    if !isempty(area_file)
        trg_geom = agd.getgeom(collect(agd.getlayer(agd.read(area_file), 0))[1])

        agd.geomdim(trg_geom) != 2 && throw(DomainError(area_file, "`area` must be a polygon."))
        agd.toWKT(agd.getspatialref(trg_geom)) != refsys && throw(DomainError("The reference systems are not uniform."))
        !agd.contains(trg_geom, src_geom) && throw(DomainError(area_file, "`area` polygon must contain the source."))
    end

    effective_infiltration *= (mean_rainfall / 10.0)^2.0

    x_source = agd.getx(src_geom, 0)
    y_source = agd.gety(src_geom, 0)
    r_source, c_source = Functions.toIndexes(dem, x_source, y_source)

    element = Leach(henry_const, volumetric_water_content, volumetric_air_content, soil_adsorption, effective_infiltration, soil_density, source_thickness, aquifer_depth, darcy_velocity, mixing_zone_depth, orthogonal_width)
    calc_kw!(element)
    calc_ldf!(element)
    calc_sam!(element)
    secondary_source_concentration = concentration * calc_LF!(element)

    daf = DAF(secondary_source_concentration, x_source, y_source, 0.0, 0.0, decay_coeff, darcy_velocity, soil_adsorption, soil_density, tera_e, orthogonal_width, time, aquifer_flow_direction, algorithm, option)
    points = !isempty(area_file) ? Functions.expand(r_source, c_source, concentration, tollerance, dem, trg_geom, daf) :
        Functions.expand(r_source, c_source, concentration, tollerance, dem, daf)

    maxR = maximum( point -> point[1], points )
    minR = minimum( point -> point[1], points )
    maxC = maximum( point -> point[2], points )
    minC = minimum( point -> point[2], points )

    geotransform = agd.getgeotransform(dem)
    geotransform[[1, 4]] = Functions.toCoords(dem, minR, minC)
    noData = Float32(agd.getnodatavalue(agd.getband(dem, 1)))
    data = fill(noData, maxR-minR+1, maxC-minC+1)
    for r in minR:maxR, c in minC:maxC
        match = findfirst( p -> p[1] == r && p[2] == c, points )
        if !isnothing(match)
            data[r - minR + 1, c - minC + 1] = points[match][3]
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, refsys, noData, output_path)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
function twodAdvectionDiffusion(x,eConn,innerNodes,outerNodes,velocity, κ = 1.0,ω=-10.0)
#
#  Solves the advection-diffusion equation in 2D
#     - ∇⋅(κ∇θ) + v⋅∇θ = q,
#  Ω is the domain between a Bspline disk and an
#  outer circle.  We prescribe zero Dirichlet conditions to the outer boundary.
#
#  As a unit test, choose a value of q such that the solution θ when
#  v = (-ωy,ωx) becomes cos(2πr)-1 (which is zero at the outer boundary and has
#  zero flux on the inner boundary).
#---------------------------------------------------------------------------78--

#  include("twodQuadratureRule.jl")
#  include("twodShape.jl")
#  include("twodBilinear.jl")
#  include("twodLinForm.jl")

  rule  = 7;   # points in quadrature formula

  function q(x::Array{Float64,2},κ,ω)
    d = size(x,1)
    ans = zeros(Float64,d,1)
    for i=1:d
      rad = sqrt(x[i,1]^2+x[i,2]^2)
      ans[i] = 4*κ*π^2*cos(2*π*rad) + 2*κ*π*sin(2*π*rad)/rad
    end
    return ans
  end

  #  Get problem dimensions
  #-----------------------------------------------------------------------------
  nNodes = size(x,1)
  nElements = size(eConn,1)

  nInnerNodes = length(innerNodes)
  nOuterNodes = length(outerNodes)
#  nDirichlet = nInnerNodes+nOuterNodes   # use for Dirichlet BCs
  nDirichlet = nOuterNodes

  #  Set the index into equation numbers
  #-----------------------------------------------------------------------------
  ideθ = zeros(Int64,nNodes,1)

  global nUnk = 0
  for i=1:nNodes
    global nUnk = nUnk+1
    ideθ[i,1] = nUnk
  end

  #  Integrate system matrices, element-by-element
  #-----------------------------------------------------------------------------
  r,s,w = twodQuadratureRule(rule)
  κ_g   = κ*ones(rule)
  one   = ones(rule)

  nElDOF   = size(eConn,2)
  nElDOFsq = nElDOF*nElDOF
  nEntries = nElements*nElDOFsq
#  II  = SharedArray{Int64}(nEntries)
#  JJ  = SharedArray{Int64}(nEntries)
#  AA  = SharedArray{Float64}(nEntries)
#  b   = SharedArray{Float64}(nNodes)
  II  = Array{Int64,1}(undef, nEntries)
  JJ  = Array{Int64,1}(undef, nEntries)
  AA  = Array{Float64,1}(undef, nEntries)
  b   = zeros(Float64,nNodes,1)

#  @sync @distributed for k=1:nElements
  for k=1:nElements
    xg  = zeros(Float64,rule,2)
    wg  = zeros(Float64,rule)
    ϕ = zeros(Float64,rule,nElDOF)
    ϕ_x = zeros(Float64,rule,nElDOF)
    ϕ_y = zeros(Float64,rule,nElDOF)

    nLocal = eConn[k,:][:]
    xLocal = x[nLocal,:]

    xg,wg,ϕ,ϕ_x,ϕ_y = twodShape( xLocal, r, s, w )

    q_g   = q(xg,κ,ω)           # forcing function evaluated at quadrature points
    u_g   = ϕ*velocity[nLocal,1]
    v_g   = ϕ*velocity[nLocal,2]

    ALoc = twodBilinear( u_g, ϕ_x, ϕ  , wg) + twodBilinear( v_g, ϕ_y, ϕ  , wg) +
           twodBilinear( κ_g, ϕ_x, ϕ_x, wg) + twodBilinear( κ_g, ϕ_y, ϕ_y, wg)
    FLoc = twodLinForm(  q_g, ϕ  ,      wg)

    index = (k-1)*nElDOFsq  # compute base index (since k could be in any order)
    lDOF = ideθ[nLocal,1][:]

    for nt = 1:nElDOF
      nTest = ideθ[nLocal[nt],1]
      for nu = 1:nElDOF
        nUnkT = ideθ[nLocal[nu],1]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkT
        AA[index] = ALoc[nt,nu]
      end

      b[nTest] = b[nTest] + FLoc[nt]
    end
  end

  A = sparse(II,JJ,AA)

  interiorNodes = convert(Array{Int64,1},1:nNodes)
#  interiorNodes = setdiff(interiorNodes,innerNodes)
  outerN = convert(Array{Int64,1},outerNodes)
  setdiff!(interiorNodes,outerN)

#  knownIndexθ = [innerNodes;outerNodes]
  knownIndexθ = outerNodes

  dirichletθ = Array{Float64,1}(undef,nDirichlet)
  for i=1:nDirichlet
    nO = knownIndexθ[i]
    dirichletθ[i] = 0.0
#    dirichletθ[i] = x[nO,1]
  end

  nUnknowns = length(interiorNodes)
  unKnownIndex = Array{UInt64,1}(undef,nUnknowns)

  global nUnknownCounter = 0
  for i ∈ interiorNodes
    global nUnknownCounter = nUnknownCounter+1
    unKnownIndex[nUnknownCounter] = i
  end

  θunk = A[unKnownIndex,unKnownIndex]\(b[unKnownIndex]-A[unKnownIndex,knownIndexθ[:]]*dirichletθ)

  θ = zeros(nNodes,1)

  global nUnknownCounter = 0
  for i ∈ interiorNodes
    global nUnknownCounter = nUnknownCounter+1
    θ[i] = θunk[nUnknownCounter]
  end

#  for i=1:nInnerNodes
#    θ[innerNodes[i]] = dirichletθ[i]
#  end

#  for i=1:nOuterNodes
#    θ[outerNodes[i]] = dirichletθ[nInnerNodes+i]
#  end
  for i=1:nOuterNodes
    θ[outerNodes[i]] = dirichletθ[i]
  end

  θexact = zeros(nNodes,1)
  for i=1:nNodes
    ri = sqrt(x[i,1]^2+x[i,2]^2)
    θexact[i] = cos(2*π*ri)-1.0
  end

  return θ,θexact
end


text:
#!/usr/bin/env julia6
using ReadGlobal

function model4(r1,r2,r3,d11,d12,d13,d22,d23,d33,v1,v2,v3)
    Threads.@threads for i in 1:length(d11)
        @inbounds v1[i] = d11[i]*r1[i] + d12[i]*r2[i] + d13[i]*r3[i]
        @inbounds v2[i] = d12[i]*r1[i] + d22[i]*r2[i] + d23[i]*r3[i]
        @inbounds v3[i] = d13[i]*r1[i] + d23[i]*r2[i] + d33[i]*r3[i]
    end
end

function model3(r1,r2,r3,d11,d12,d13,d22,d23,d33,w12,w13,w23,v1,v2,v3)
    Threads.@threads for i in 1:length(d11)
        @inbounds v1[i] = d11[i]*r1[i] + (d12[i]-w12[i])*r2[i] + (d13[i]-w13[i])*r3[i]
        @inbounds v2[i] = (d12[i]+w12[i])*r1[i] + d22[i]*r2[i] + (d23[i]-w23[i])*r3[i]
        @inbounds v3[i] = (d13[i]+w13[i])*r1[i] + (d23[i]+w23[i])*r2[i] + d33[i]*r3[i]
    end
end

function model5(r1,r2,r3,w12,w13,w23,v1,v2,v3)
    Threads.@threads for i in 1:length(r1)
        @inbounds v1[i] = -w12[i]*r2[i] - w13[i]*r3[i]
        @inbounds v2[i] = w12[i]*r1[i] - w23[i]*r3[i]
        @inbounds v3[i] = w13[i]*r1[i] + w23[i]*r2[i]
    end
end

function fluxmodels2()

    N, Fil = getnfilter()

    d11 = "$(Fil)D11_N$N"
    d22 = "$(Fil)D22_N$N"
    d33 = "$(Fil)D33_N$N"
    d12 = "$(Fil)D12_N$N"
    d13 = "$(Fil)D13_N$N"
    d23 = "$(Fil)D23_N$N"

    w12 = "$(Fil)W12_N$N"
    w13 = "$(Fil)W13_N$N"
    w23 = "$(Fil)W23_N$N"

    r1 = "$(Fil)drhodx_N$N"
    r2 = "$(Fil)drhody_N$N"
    r3 = "$(Fil)drhodz_N$N"

    doinchunks(model4,input=(r1,r2,r3,d11,d12,d13,d22,d23,d33),output=("4modelflux1","4modelflux2","4modelflux3"))

    doinchunks(model3,input=(r1,r2,r3,d11,d12,d13,d22,d23,d33,w12,w13,w23),output=("3modelflux1","3modelflux2","3modelflux3"))

    doinchunks(model5,input=(r1,r2,r3,w12,w13,w23),output=("5modelflux1","5modelflux2","5modelflux3"))

end

fluxmodels2()