text:
const ScalarPotential{T, N, S} = Union{ElectricPotential{T, N, S}, WeightingPotential{T, N, S}, PointTypes{T, N, S}, ChargeDensity{T, N, S}}

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T}) where {T}
    @assert int.right != 0 "Right boundary of φ interval is not allowed to be 0"
    Potential::Type = typeof(wp)
    l::Int = length( axφ )
    Δφ::T = int.right - int.left
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(0, 2π)
    n::Int = Int(round(T(2π) / Δφ, sigdigits = 6))
    new_ticks::Vector{T} = Vector{T}(undef, l * n)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), l * n, size(wp, 3))
    for idx_n in 1:n
        for il in 1:l
            idx::Int = il + (idx_n - 1) * l
            new_ticks[idx] = (idx_n - 1) * Δφ + axφ.ticks[il]
            new_pot[:, idx, :] = wp[:, il, :]
        end
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return new_wp
end

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T}) where {T}
    @assert int.right != 0 "Right boundary of φ interval is not allowed to be 0"
    Potential::Type = typeof(wp)
    l::Int = 2 * length( axφ ) - 2
    Δφ::T = int.right - int.left
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(int.left, int.right + Δφ)
    new_ticks::Vector{T} = Vector{T}(undef, l)
    new_ticks[1:length( axφ )] = collect(axφ.ticks)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), l, size(wp, 3))
    new_pot[:, 1:length( axφ ), :] = wp.data[:, :, :]
    for i in 1:(length(axφ) - 2)
        new_ticks[length(axφ) + i] = 2 * int.right - axφ.ticks[length(axφ) - i]
        new_pot[:, length(axφ) + i, :] = wp[:, length(axφ) - i, :]
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return get_2π_potential(new_wp, new_axφ, new_int)
    # return new_wp
end

function extend_2D_to_3D_by_n_points(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
        n_points_in_φ::Int ) where {T}
    Potential::Type = typeof(wp)
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(0, 2π)
    new_ticks::Vector{T} = Vector{T}(undef, n_points_in_φ)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), n_points_in_φ, size(wp, 3))
    Δφ::T = 2π / n_points_in_φ
    for i in 1:n_points_in_φ
        new_ticks[i] = (i - 1) * Δφ
        new_pot[:, i, :] = wp[:, 1, :]
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return new_wp
end

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}; n_points_in_φ::Union{Missing, Int} = missing) where {T}
    axφ::DiscreteAxis{T} = wp.grid[2]
    int::Interval = axφ.interval
    if int.right == 0 && length(axφ) == 1 # 2D
        if ismissing(n_points_in_φ)
            error(ArgumentError, ": First Argument is a 2D potential (only 1 point in φ). User has to set the keyword `n_points_in_φ::Int` (e.g. `n_points_in_φ = 18`) in order to get a 3D potential.")
        else
            return extend_2D_to_3D_by_n_points(wp, axφ, int, n_points_in_φ)
        end
    else
        return get_2π_potential(wp, axφ, int)
    end
end


text:
module Test2D

using Test
using Trixi2Img

include("test_trixi2img.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "2D" begin
  run_trixi(joinpath("2d", "elixir_advection_extended.jl"), maxiters=1)

  @testset "uniform mesh as PNG" begin
    test_trixi2img("solution_000000.h5", outdir)
        # To be able to compare hashes, we would need to make sure that exactly the same versions
        # of all required libraries are installed everywhere, which does not seem to be a good option
        # right now.
        # hashes=[("solution_000000_scalar.png", "f7da610fc8ab16aa691b1b0ecd956d67d6f68e06")])
  end

  # PDF tests do not work on Windows
  if !Sys.iswindows()
    @testset "uniform mesh as PDF with grid lines" begin
      test_trixi2img("solution_000000.h5", outdir,
          format=:pdf, grid_lines=true)
          # To be able to compare hashes, we would need to make sure that exactly the same versions
          # of all required libraries are installed everywhere, which does not seem to be a good option
          # right now.
          # hashes=[("solution_000000_scalar.pdf", "e62024c18ae010748f8c7b840bf7475ae398799a")])
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end

ext: jl


text:
@testset "RadialMapComponent (Nx,p) = (1,0)" begin

    C = RadialMapComponent(1,0)

    @test C.Nx == 1
    @test C.p == 0
    @test size(C.ξ, 1)==1
    @test size(C.ξ[1], 1)==0

    @test size(C.σ, 1)==1
    @test size(C.σ[1], 1)==0

    @test size(C.coeff,1)==1
    @test size(C.coeff[1],1)==2

    C.coeff[1]=[2.0;1.0]
    @test C(1.0)==3.0
end

@testset "RadialMapComponent (Nx,p) = (1,1)" begin

    C = RadialMapComponent(1,1)

    @test C.Nx == 1
    @test C.p == 1

    @test size(C.ξ,1) ==1
    @test size(C.ξ[1],1) ==1+2

    @test size(C.σ,1) ==1
    @test size(C.σ[1],1) ==1+2

    @test size(C.ξ,1) ==1
    @test size(C.ξ[1],1) ==1+2

    @test size(C.coeff,1) ==1
    @test size(C.coeff[1],1) ==1+3

    C.coeff[1] .= [1.0; 2.0; 3.0; 1.0]

    x = randn()
    @test abs(C(x) - (1.0 + 2.0*ψ₀(0.0,1.0)(x) +
                     3.0*ψj(0.0,1.0)(x) + 1.0*ψpp1(0.0,1.0)(x)))< 1e-10

    @test abs(C([x]) - (1.0 + 2.0*ψ₀(0.0,1.0)(x) +
                      3.0*ψj(0.0,1.0)(x) + 1.0*ψpp1(0.0,1.0)(x)))< 1e-10

end

@testset "RadialMapComponent (Nx,p) = (Nx,0) with Nx>1" begin

    C = RadialMapComponent(3, 0)
    x1 = randn()
    x2 = randn()
    x3 = randn()

    @test C.Nx == 3
    @test C.p == 0

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==0
    @test size(C.σ[1],1)==0
    @test size(C.coeff[1],1)==1

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==1

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==2

    @test C(1.0) == 0.0

    C.coeff[1][1] = 1.5

    @test isapprox(C([x1;0.0;0.0]), 1.5*x1, atol = 1e-10)

    C.coeff[2][1] = 2.0
    @test isapprox(C([x1;x2;0.0]), 1.5*x1 + 2.0*x2, atol = 1e-10)

    C.coeff[3][1] = 3.0
    @test isapprox(C([x1; x2; x3]), 1.5*x1 + 2.0*x2 + 3.0, atol = 1e-10)
    @test isapprox(C(x1), 1.5*x1+2.0*x1+3.0, atol = 1e-10)

    C.coeff[3][2] = 4.0
    @test isapprox(C([x1; x2; x3]), 1.5*x1 + 2.0*x2 + 3.0 + 4.0*x3, atol = 1e-10)
    @test isapprox(C(x1), 1.5*x1+2.0*x1+3.0 + 4.0*x1, atol = 1e-10)

end

@testset "RadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p>1" begin

    C = RadialMapComponent(3, 2)
    x1 = randn()
    x2 = randn()
    x3 = randn()

    @test C.Nx == 3
    @test C.p == 2

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==2
    @test size(C.σ[1],1)==2
    @test size(C.coeff[1],1)==2+1

    @test size(C.ξ[2],1)==2
    @test size(C.σ[2],1)==2
    @test size(C.coeff[2],1)==2+1

    @test size(C.ξ[3],1)==2+2
    @test size(C.σ[3],1)==2+2
    @test size(C.coeff[3],1)==2+3

    @test C.ξ[1]==zeros(2)
    @test C.ξ[2]==zeros(2)
    @test C.ξ[3]==zeros(2+2)

    @test C.σ[1]==ones(2)
    @test C.σ[2]==ones(2)
    @test C.σ[3]==ones(2+2)

    @test C.coeff[1]==zeros(2+1)
    @test C.coeff[2]==zeros(2+1)
    @test C.coeff[3]==zeros(2+3)

    C.coeff[1] .= ones(3)

    @test isapprox(C(x1), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x1), atol = 1e-10)

    @test isapprox(C([x1; 0.0; 0.0]), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x1), atol = 1e-10)

    C.coeff[2] .= 2*ones(3)

    @test isapprox(C([x1; x2; 0.0]), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x1) + ui(2,C.ξ[2], C.σ[2], C.coeff[2])(x2), atol = 1e-10)

    C.coeff[3] .= 3*ones(5)

    @test isapprox(C([x1; x2; x3]), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x1) +
                                 ui(2,C.ξ[2], C.σ[2], C.coeff[2])(x2) +
                                 uk(2,C.ξ[3], C.σ[3], C.coeff[3])(x3), atol = 1e-10)

    @test isapprox(C(x1),   ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x1) +
                            ui(2,C.ξ[2], C.σ[2], C.coeff[2])(x1) +
                            uk(2,C.ξ[3], C.σ[3], C.coeff[3])(x1), atol = 1e-10)
end

@testset "Component function for RadialMapComponent type" begin

    C = RadialMapComponent(1, 2)

    @test C.Nx == 1
    @test C.p == 2

    @test  TransportBasedInference.component(C, 1).ξk == zeros(2+2)
    @test  TransportBasedInference.component(C, 1).σk == ones(2+2)
    @test  TransportBasedInference.component(C, 1).coeffk == zeros(2+3)

    C = RadialMapComponent(5, 2)

    @test C.Nx == 5
    @test C.p == 2

    @test  TransportBasedInference.component(C, 5).ξk == zeros(2+2)
    @test  TransportBasedInference.component(C, 5).σk == ones(2+2)
    @test  TransportBasedInference.component(C, 5).coeffk == zeros(2+3)

    @test  TransportBasedInference.component(C, 2).ξi == zeros(2)
    @test  TransportBasedInference.component(C, 2).σi == ones(2)
    @test  TransportBasedInference.component(C, 2).coeffi == zeros(2+1)
end

@testset "Verify off_diagonal function" begin

    # Nx=1  & p=0
    C = RadialMapComponent(1,0)
    C.coeff[1] .= randn(2)

    @test off_diagonal(C, randn()) == 0.0

    # Nx=1 & p= 3
    C = RadialMapComponent(1, 3)
    C.coeff[1] .= rand(6)

    @test off_diagonal(C, randn()) == 0.0

    # Nx=3 & p = 0
    C = RadialMapComponent(3, 0)
    a1 =randn()
    a2 = randn()
    a3 = rand(2)
    C.coeff[1] .= a1
    C.coeff[2] .= a2
    C.coeff[3] .= a3

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10

    # Nx=3 & p = 3
    C = RadialMapComponent(3, 3)
    for i=1:2
    C.ξ[i] .= randn(3)
    C.σ[i] .= rand(3)
    end

    C.ξ[3] .= randn(5)
    C.σ[3] .= rand(5)

    a1 =randn(4)
    a2 = randn(4)
    a3 = rand(6)
    C.coeff[1] .= a1
    C.coeff[2] .= a2
    C.coeff[3] .= a3

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10
end

@testset "extract and modify coefficients of RadialMapComponent" begin
    # Nx=1 and p=0
    C = RadialMapComponent(1, 0)

    @test C.Nx == 1
    @test C.p == 0

    modifycoeff!(C, [1.0;2.0])
    @test extractcoeff(C) == [1.0; 2.0]

    C = RadialMapComponent(1, 3)
    modifycoeff!(C, [1.0;2.0; 3.0; 4.0; 5.0; 6.0])
    @test extractcoeff(C) == [1.0;2.0; 3.0; 4.0; 5.0; 6.0]

    # Nx =3 and p=2
    A = collect(1.0:1.0:11.0)
    C = RadialMapComponent(3, 2)

    modifycoeff!(C, A)
    @test  extractcoeff(C) == A

    C = RadialMapComponent(3, 2)
    C.coeff[1] = [1.0; 2.0; 3.0]
    C.coeff[2] = [4.0; 5.0; 6.0]
    C.coeff[3] = [7.0; 8.0; 9.0; 10.0; 11.0]

    @test extractcoeff(C) ==A

    # Nx = 3 and p = 0

    C = RadialMapComponent(3, 0)

    @test C.Nx == 3
    @test C.p == 0

    C.coeff[1] .= [1.0]
    C.coeff[2] .= [2.0]
    C.coeff[3] .= [3.0; 4.0]

    @test extractcoeff(C)==collect(1.0:4.0)

    modifycoeff!(C, collect(5.0:8.0))

    @test extractcoeff(C)==collect(5.0:8.0)
end


text:
mutable struct DiscreteSystem
    probability::Vector{Vector{Int64}}
    table::Vector{Array{Float64, N} where N}
    incidence::SparseMatrixCSC{Float64,Int64}
    incidenceTranspose::SparseMatrixCSC{Float64,Int64}
    data::String
end

mutable struct DiscreteTreeGraph
    Nvariable::Int64
    Nfactor::Int64
    Nlink::Int64
    root::Int64
    virtualMessage::Float64
    state::Array{Int64,1}
    messageDirect::Vector{Vector{Float64}}
    rowForward::Vector{Vector{Int64}}
    rowBackward::Vector{Vector{Int64}}
    colForward::Vector{Vector{Int64}}
    colBackward::Vector{Vector{Int64}}
    incomingToFactor::Vector{Vector{Int64}}
    incomingToVariable::Vector{Vector{Int64}}
    iterateFactor::Array{Int64,1}
    iterateVariable::Array{Int64,1}
    passFactorVariable::Int64
    passVariableFactor::Int64
    forward::Bool
    backward::Bool
end

mutable struct DiscreteInference
    fromFactor::Array{Int64,1}
    toVariable::Array{Int64,1}
    messageFactorVariable::Vector{Vector{Float64}}
    fromVariable::Array{Int64,1}
    toFactor::Array{Int64,1}
    messageVariableFactor::Vector{Vector{Float64}}
    marginal::Vector{Vector{Float64}}
end

struct DiscreteTreeModel
    graph::DiscreteTreeGraph
    inference::DiscreteInference
    system::DiscreteSystem
end

########## Form a graph with continuous variables and initialize messages ##########
function discreteTreeModel(
    args...;
    virtual::Float64 = 1.0,
    root::Int64 = 1)

    if checkFileOrArguments(args)
        system = readDiscreteFile(args)
    else
        system = readDiscreteArguments(args)
    end

    graph, inference = makeDiscreteTreeGraph(system, virtual, root)

    return DiscreteTreeModel(graph, inference, system)
end

########## Read in-Julia discrete system model ##########
function readDiscreteArguments(args)
    if isa(args[1], Dict)
        probability = []; table = []
        @inbounds for (key, value) in args[1]
            if key == "probability"
                for k in value
                    push!(probability, k)
                end
            end
            if key == "table"
                for k in value
                    push!(table, k)
                end
            end
        end
    else
        probability = args[1]
        table = args[2]
    end

    incidence, incidenceTranspose = incidenceMatrix(probability)

    return DiscreteSystem(probability, table, incidence, incidenceTranspose, "noname")
end

########## Load from HDF5 or XLSX files ##########
function readDiscreteFile(args)
    fullpath, extension, dataname = checkImportFile(args)

    #### Read from HDF5 or XLSX file
    if extension == ".h5"
        list = h5open(fullpath, "r")
        probability = []; table = Any[]
        probabilityH5 = h5read(fullpath, "/probability")
        tableH5 = h5read(fullpath, "/table")
        for key in sort(collect(keys(probabilityH5)))
            push!(probability, vec(probabilityH5[key]))
            if size(tableH5[key], 2) == 1
                push!(table, vec(tableH5[key]))
            else
                push!(table, tableH5[key])
            end
        end
    elseif extension == ".xlsx"
        xf = XLSX.openxlsx(fullpath, mode = "r")
        if "graphical model" in XLSX.sheetnames(xf)
            list = xf["graphical model"][:][1:end, :]

            probability = Any[]; table = Any[]; cpt = []
            Nlist = size(list, 1);  flag = false; last = 0
            @inbounds for i = 1:Nlist - 1
                if !ismissing(list[i, 1])
                    if list[i, 1] == "conditional independence"
                        for (q, k) in enumerate(list[i, 2:end])
                            if ismissing(k)
                                last = q
                                break
                            end
                        end
                        push!(probability, list[i, 2:last])
                    end

                    if list[i, 1] == "discrete variable states"
                        cpt = zeros(tuple(list[i, 2:last]...))
                    end
                    if list[i, 1] == "conditional probability table"
                        flag = true
                    end
                end
                if flag && !ismissing(list[i, 2])
                    idx = CartesianIndex(tuple(list[i, 2:last]...))
                    cpt[idx] = list[i, last + 1]
                end
                if !ismissing(list[i + 1, 1]) && list[i + 1, 1] == "conditional independence"
                    push!(table, cpt)
                    flag = false
                end
                if i + 1 == Nlist && !ismissing(list[i + 1, 2])
                    idx = CartesianIndex(tuple(list[i + 1, 2:last]...))
                    cpt[idx] = list[i + 1, last + 1]
                    push!(table, cpt)
                end
            end
        else
            throw(ErrorException("Error opening sheet: graphical model."))
        end
    else
        error("The input data is not a valid format.")
    end

    incidence, incidenceTranspose = incidenceMatrix(probability)

    return DiscreteSystem(probability, table, incidence, incidenceTranspose, dataname)
end

########## Incidence matrix ##########
function incidenceMatrix(probabilities)
    row = Int64[]; col = Int64[]

    @inbounds for (k, probability) in enumerate(probabilities)
        for variable in probability
            push!(row, k)
            push!(col, variable)
        end
    end
    incidence = sparse(row, col, 1)
    incidenceTranspose = sparse(col, row, 1)

    return incidence, incidenceTranspose
end

########## Produce the graphical model ##########
function makeDiscreteTreeGraph(system, virtualMessage, root)
    ### Number of factor and variable nodes
    Nfactor, Nvariable = size(system.incidence)

    ### Pass through probabilities
    Nlink = 0
    states = fill(0, Nvariable)
    messageDirect = [Float64[] for i = 1:Nvariable]
    rowForward =  [Int64[] for i = 1:Nfactor]
    colForward = [Int64[] for i = 1:Nvariable]
    @inbounds for (k, probability) in enumerate(system.probability)
        for (kk, variable) in enumerate(probability)
            if states[variable] == 0
                states[variable] = size(system.table[k], kk)
                messageDirect[variable] = fill(virtualMessage, states[variable])
            elseif states[variable] != size(system.table[k], kk)
                error("The possible states of the variable $variable do not agree through tables.")
            end
            if length(probability) != 1
                push!(rowForward[k], variable)
                push!(colForward[variable], k)
                Nlink += 1
            else
                messageDirect[variable] .*= system.table[k]
            end
        end
    end
    rowBackward = deepcopy(rowForward)
    colBackward = deepcopy(colForward)

    ### Initialize the first forward step
    iterateVariable = fill(0, Nvariable)
    counter = 0
    @inbounds for col = 1:Nvariable
        if length(colForward[col]) == 1 && col != root
            counter += 1
            iterateVariable[counter] = col
        end
    end
    resize!(iterateVariable, counter)

    ### Initialize data
    fromVariable = fill(0, Nlink)
    toFactor = fill(0, Nlink)
    messageVariableFactor = [Float64[] for i = 1:Nlink]

    fromFactor = fill(0, Nlink)
    toVariable = fill(0, Nlink)
    messageFactorVariable = [Float64[] for i = 1:Nlink]

    marginal = [Float64[] for i = 1:Nvariable]
    incomingToVariable = [Int64[] for i = 1:Nvariable]
    incomingToFactor = [zeros(length(system.probability[i])) for i = 1:Nfactor]
    passFactorVariable = 0
    passVariableFactor = 0
    iterateFactor = Int64[]

    return DiscreteTreeGraph(Nvariable, Nfactor, Nlink, root, virtualMessage, states, messageDirect,
        rowForward, rowBackward, colForward, colBackward, incomingToFactor, incomingToVariable,
        iterateFactor, iterateVariable, passFactorVariable, passVariableFactor, true, true),
           DiscreteInference(fromFactor, toVariable, messageFactorVariable, fromVariable, toFactor, messageVariableFactor, marginal)
end


text:
mutable struct Remember{V} <: State{V}
    value::V
    done::Bool
end

Remember(; unit, _value, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    v = _value
    Remember{V}(v, false)
end

constructortags(::Val{:Remember}) = (:unit,)

genvartype(v::VarInfo, ::Val{:Remember}; V, _...) = @q Remember{$V}

gendefault(v::VarInfo, ::Val{:Remember}) = begin
    N = gettag(v, :_type)
    U = gettag(v, :unit)
    i = gettag(v, :init)
    if isnothing(i)
        @q zero($C.valuetype($N, $U))
    else
        @q $C.unitfy($C.value($i), $C.value($U))
    end
end

genupdate(v::VarInfo, ::Val{:Remember}, ::MainStep; kw...) = begin
    w = gettag(v, :when)
    @gensym s
    @q let $s = $(symstate(v))
        if !($s.done) && $C.value($w)
            $(genstore(v; when=false))
            $s.done = true
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
begin
    using PlutoUI
	using ExtendableGrids
	using SimplexGridFactory
	using PlutoVista
	using GridVisualize
        using Triangulate
        isdefined(Main,:PlutoRunner) ? 	default_plotter!(PlutoVista) : default_plotter!(nothing)
end

# ╔═╡ 07194f25-5735-452f-9bed-cf791958d44d
md"""
## Grid with quadratic boundary subregions at electrode
"""

# ╔═╡ d0f3483a-2bf4-4e2d-80b0-6c869b45cda8
function qxygrid(;nref=0,xymax=1.0,hxymin=0.1,hxymax=0.1)

Xp=geomspace(0,xymax,hxymin,hxymax)
Xm=-reverse(Xp)
X=glue(Xm,Xp)
gridxy=simplexgrid(X,X)
fill!(gridxy[BFaceRegions],5)
gridxy
end

# ╔═╡ fd8f6c48-0030-4f6f-9fbf-9f2340decbf2
g=qxygrid(hxymin=0.05, hxymax=0.2)

# ╔═╡ bbeeebc0-ae16-4022-975d-35588e62166a
gridplot(g)

# ╔═╡ e6c7e47d-6159-4bb6-a71d-7c633e2fefbd
function qxyzgrid(gridxy ;zmax=1.0,hzmin=0.05,hzmax=0.2)
	Z=geomspace(0,zmax,hzmin,hzmax)
	gridxyz=simplexgrid(gridxy,Z,top_offset=1)
	xymax=gridxy[Coordinates][1,end]

	bfacemask!(gridxyz,[-xymax,-xymax,0],[0,0,0],1)
    bfacemask!(gridxyz,[0,0,0],[xymax,xymax,0],2)
    bfacemask!(gridxyz,[0,-xymax,0],[xymax,0,0],3)
    bfacemask!(gridxyz,[-xymax,0,0],[0,xymax,0],4)
    gridxyz
end

# ╔═╡ 9800c97e-c25b-4d08-a014-cccd746f7d71
gxyz=qxyzgrid(g)

# ╔═╡ 9a219111-0275-4cd9-b97e-648d3fcfcbb9
vis=GridVisualizer(dim=3,resolution=(400,300));vis

# ╔═╡ 85be1677-87ff-49dc-af9a-557e575bc55f
@bind z Slider(0:0.01:1,show_value=true,default=0.5)

# ╔═╡ 33ecc78c-16ac-46ac-8d83-cbb26288a6fb
gridplot!(vis,gxyz,show=true,zplanes=[z],xplanes=[1],yplanes=[1],outlinealpha=0.3)

# ╔═╡ 9e640256-61a4-4fb9-a449-d5f948fb2d26
md"""
## Grid with inner circle at electrode
"""

# ╔═╡ bf0004f4-b3c8-4385-b87d-df2ef1420408
md"""
### Plain triangulation
"""

# ╔═╡ 1b70e44b-cbe8-4812-b267-ef9bce20a5b7
function cxygrid(;maxvol=0.01,nref=0,xyzmax=1,rad=0.5)
	builder=SimplexGridBuilder(Generator=Triangulate)

	regionpoint!(builder,(xyzmax-0.1,xyzmax-0.1))
	cellregion!(builder,1)
	rect2d!(builder,[-xyzmax,-xyzmax],[xyzmax,xyzmax],facetregions=1)

	facetregion!(builder,3)
	cellregion!(builder,2)
	regionpoint!(builder, (0,0))
	circle!(builder,(0,0), rad,n=20)
	simplexgrid(builder,maxvolume=maxvol)
end

# ╔═╡ 10ebaeab-3ece-480f-9c4d-a676814e2f7f
gcxy=cxygrid()

# ╔═╡ b711885e-cc61-4e08-9367-437731323863
gridplot(gcxy)

# ╔═╡ e4029d2c-18b7-4592-b47d-294e4af93497
gcxy.components

# ╔═╡ 0c3f5b1c-7bcf-491e-8318-73f2c880e0b5
function cxyzgrid(gcxy,zmax=1.0,hzmin=0.05,hzmax=0.2)
	Z=geomspace(0,zmax,hzmin,hzmax)
	gridxyz=simplexgrid(gcxy,Z,bot_offset=0)
	xymax=maximum(gcxy[Coordinates][1,:])

	bfacemask!(gridxyz,[-xymax,-xymax,zmax],[xymax,xymax,zmax],5,allow_new=false)
	gridxyz

end

# ╔═╡ f9130b5d-b4ea-4ebd-b73f-82a3c1bc1307
gcxyz=cxyzgrid(gcxy)

# ╔═╡ 73847699-40c0-4043-a16f-37883def6858
gcxyz.components

# ╔═╡ bff0a3f1-c545-4fcc-8be8-ef460f8479bd
visc=GridVisualizer(dim=3,resolution=(400,300))

# ╔═╡ 66ce2939-2957-43e0-b7ea-cdd10d9fc17e
@bind zc Slider(0:0.01:1,show_value=true,default=0.5)

# ╔═╡ a310027d-3f9e-4de5-bc9e-5457cb19eef5
gridplot!(visc,gcxyz,show=true,zplanes=[zc],xplanes=[1],yplanes=[1],outlinealpha=0.3)

# ╔═╡ 752bdb8f-de5a-4863-b2d1-53d69aff7dcb
md"""
### Grid with anisotropic local refinement at electrode
"""

# ╔═╡ 5aef30a0-8712-4c6c-9465-25da0624b408
function grxy(;nang=40,nxy=15, r=0.5, dr=0.1,hrmin=0.01,hrmax=0.05,xymax=1.0)
		rad1=geomspace(r-dr,r,hrmax,hrmin)
		rad2=geomspace(r,r+dr,hrmin,hrmax)
		ang=range(0,2π,length=nang)
		δr=2π*r/nang
	    maxvol=0.5*δr^2
		ring1=ringsector(rad1,ang)
		ring1[CellRegions].=1
		ring2=ringsector(rad2,ang)
		ring1[CellRegions].=2
		ring=glue(ring1,ring2,breg=3)

		binner=SimplexGridBuilder(Generator=Triangulate)
		regionpoint!(binner,(0,0))
		cellregion!(binner,2)

		facetregion!(binner,2)
		bregions!(binner,ring,[1])
		ginner=simplexgrid(binner,maxvolume=maxvol,nosteiner=true)
		ginner[CellRegions].=2

		bouter=SimplexGridBuilder(Generator=Triangulate)

		holepoint!(bouter,(0,0))

		facetregion!(bouter,3)
		bregions!(bouter,ring,[2])

	regionpoint!(bouter,(xymax-0.1,xymax-0.1))
	cellregion!(bouter,3)
	rect2d!(bouter,[-xymax,-xymax],[xymax,xymax],facetregions=1,nx=nxy,ny=nxy)

	gouter=simplexgrid(bouter;maxvolume=maxvol*2,nosteiner=true)

	glue(glue(gouter,ring),ginner)

end

# ╔═╡ 4debce38-96c5-4661-a965-1ca723fa8a36
rxygrid=grxy()

# ╔═╡ 52f4c693-0218-4f56-ac95-abd6d99a0b67
gridplot(rxygrid)

# ╔═╡ f43bd8ff-a2ed-4740-a61e-4f5d2c615c10
grxyz=cxyzgrid(rxygrid)

# ╔═╡ b9ced73f-e7a0-4573-9413-c2a231d21c78
  visr=GridVisualizer(dim=3,resolution=(400,300))

# ╔═╡ 2c285e4f-9f10-474b-9482-83a5c4cbfe09
@bind zr Slider(0:0.01:1,show_value=true,default=0.5)

# ╔═╡ 97d8db91-16e6-4070-81a0-7bcaff6cc9f6
gridplot!(visr,grxyz,show=true,zplanes=[zr],xplanes=[1],yplanes=[1],outlinealpha=0.3)

# ╔═╡ 78dba5b2-52c9-40cd-bc1d-d5c343271f97
html"""<hr> """

# ╔═╡ b9cc0359-7286-4c02-ba10-35303da26a50
TableOfContents()

# ╔═╡ 605c914d-607b-4d2e-80c5-d14cb6918e32
md"""
    begin
       using Pkg
       Pkg.activate(mktempdir())
       Pkg.add(["PlutoUI","Revise","Triangulate"])
       using Revise
	   Pkg.develop(["ExtendableGrids","SimplexGridFactory",
           "GridVisualize","PlutoVista"])
    end
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PlutoVista = "646e1f28-b900-46d7-9d87-d554eb38a413"
SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"

[compat]
ExtendableGrids = "~0.8.7"
GridVisualize = "~0.3.9"
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.6"
SimplexGridFactory = "~0.5.9"
Triangulate = "~2.1.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[ElasticArrays]]
deps = ["Adapt"]
git-tree-sha1 = "a0fcc1bb3c9ceaf07e1d0529c9806ce94be6adf9"
uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
version = "1.2.9"

[[ExtendableGrids]]
deps = ["AbstractTrees", "Dates", "DocStringExtensions", "ElasticArrays", "InteractiveUtils", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "1e8e50f054057f23e908fbd6935766dca6293cc2"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.7"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[GridVisualize]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GeometryBasics", "LinearAlgebra", "Observables", "OrderedCollections", "PkgVersion", "Printf", "Requires", "StaticArrays"]
git-tree-sha1 = "925ba2f11df005d894b113292d32fca9afe3f8c8"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.9"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "f19e978f81eca5fd7620650d7dbea58f825802ee"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[PlutoVista]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "GridVisualize", "UUIDs"]
git-tree-sha1 = "34fc7e41e6eefa58fef0786ab62a20262df88764"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.6"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SimplexGridFactory]]
deps = ["DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GridVisualize", "LinearAlgebra", "Printf", "Test"]
git-tree-sha1 = "af52ec74a4b6cfcc5b6d60d259099fa0596de2c1"
uuid = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
version = "0.5.9"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[Triangle_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bfdd9ef1004eb9d407af935a6f36a4e0af711369"
uuid = "5639c1d2-226c-5e70-8d55-b3095415a16a"
version = "1.6.1+0"

[[Triangulate]]
deps = ["DocStringExtensions", "Libdl", "Printf", "Test", "Triangle_jll"]
git-tree-sha1 = "2b4f716b192c0c615d96d541ee029e85666388cb"
uuid = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
version = "2.1.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╟─07194f25-5735-452f-9bed-cf791958d44d
# ╠═d0f3483a-2bf4-4e2d-80b0-6c869b45cda8
# ╠═fd8f6c48-0030-4f6f-9fbf-9f2340decbf2
# ╠═bbeeebc0-ae16-4022-975d-35588e62166a
# ╠═e6c7e47d-6159-4bb6-a71d-7c633e2fefbd
# ╠═9800c97e-c25b-4d08-a014-cccd746f7d71
# ╠═9a219111-0275-4cd9-b97e-648d3fcfcbb9
# ╠═85be1677-87ff-49dc-af9a-557e575bc55f
# ╠═33ecc78c-16ac-46ac-8d83-cbb26288a6fb
# ╟─9e640256-61a4-4fb9-a449-d5f948fb2d26
# ╟─bf0004f4-b3c8-4385-b87d-df2ef1420408
# ╠═1b70e44b-cbe8-4812-b267-ef9bce20a5b7
# ╠═10ebaeab-3ece-480f-9c4d-a676814e2f7f
# ╠═b711885e-cc61-4e08-9367-437731323863
# ╠═e4029d2c-18b7-4592-b47d-294e4af93497
# ╠═0c3f5b1c-7bcf-491e-8318-73f2c880e0b5
# ╠═f9130b5d-b4ea-4ebd-b73f-82a3c1bc1307
# ╠═73847699-40c0-4043-a16f-37883def6858
# ╠═bff0a3f1-c545-4fcc-8be8-ef460f8479bd
# ╠═66ce2939-2957-43e0-b7ea-cdd10d9fc17e
# ╠═a310027d-3f9e-4de5-bc9e-5457cb19eef5
# ╟─752bdb8f-de5a-4863-b2d1-53d69aff7dcb
# ╠═5aef30a0-8712-4c6c-9465-25da0624b408
# ╠═4debce38-96c5-4661-a965-1ca723fa8a36
# ╠═52f4c693-0218-4f56-ac95-abd6d99a0b67
# ╠═f43bd8ff-a2ed-4740-a61e-4f5d2c615c10
# ╠═b9ced73f-e7a0-4573-9413-c2a231d21c78
# ╠═2c285e4f-9f10-474b-9482-83a5c4cbfe09
# ╠═97d8db91-16e6-4070-81a0-7bcaff6cc9f6
# ╟─78dba5b2-52c9-40cd-bc1d-d5c343271f97
# ╠═b9cc0359-7286-4c02-ba10-35303da26a50
# ╠═605c914d-607b-4d2e-80c5-d14cb6918e32
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
##############################################################################
# GAMS specific functions
##############################################################################

function reshop_setup_gams()
    ctx = ccall((:ctx_alloc, libreshop), Ptr{context}, (Cuint,), 0)

    gamscntr_template_file = joinpath(solverdata_dir, "gamscntr.dat")

    if !isfile(gamscntr_template_file)
        reshop_init_gams_solverdata()
    end
    if !isfile(gamscntr_template_file)
        error("Could not create template GAMS control file! Make sure that GAMS is properly installed and available via the system path")
    end

    gamscntr_template = read(gamscntr_template_file, String)

    gams_dir = mktempdir(pwd())
    cur_dir = gams_dir

    open(joinpath(gams_dir, "gamscntr.dat"), "w") do gamscntr_file
        println(gamscntr_file, replace(gamscntr_template, r"@@SUB@@" => cur_dir))
    end

    # we need an empty Matrixfile
    touch(joinpath(cur_dir, "gamsmatr.dat"))

    res = ccall((:gams_setgamscntr, libreshop), Cint, (Ptr{context}, Cstring), ctx, joinpath(cur_dir, "gamscntr.dat"))
    res != 0 && error("return code $res from ReSHOP")

    # hm bad hack
    gamsdir = split(gamscntr_template, term_str)[29]

    CONFIG[:debug] && println("DEBUG: gamsdir is ``$gamsdir''")

    res = ccall((:gams_setgamsdir, libreshop), Cint, (Ptr{context}, Cstring), ctx, gamsdir)
    res != 0 && error("return code $res from ReSHOP")

    ENV["PATH"] *= ":" * gamsdir

    # This is needed to prevent the listing of the Process directory
    ENV["DEBUG_PGAMS"] = '0'

    return (ctx, gams_dir)
end

function reshop_init_gams_solverdata(force=false)
    curdir = pwd()

    if !isdir(solverdata_dir)
        error("No directory named $solverdata_dir. ")
    end

    gamscntr = joinpath(solverdata_dir, "gamscntr.dat")
    if isfile(gamscntr)
        !force  && return
    end

    gms_file = joinpath(solverdata_dir, "dummy.gms")

    println("Running the gams command to initial bootstrap file\n")

    mktempdir(solverdata_dir) do substr
        cd(solverdata_dir)
        # use run(`gams $gms_file scrdir=$substr lo=0 curdir=$substr`) ?
        run(`gams $gms_file scrdir=$substr`)

        open(gamscntr, "w") do out_gamscntr
            input = read(joinpath(substr, "gamscntr.dat"), String)
            # the order of subsitution matters here
            input = replace(input, substr => "@@SUB@@")
            input = replace(input, pwd() => "@@SUB@@")
            println(out_gamscntr, input)
        end
    end

    println("\nGAMS boostrap file successfully generated")
    cd(curdir)

end


text:
"""

```
shellradius(i,step)
```

Compute the point in which the radius comprises half of the
volume of the shell.
"""
function shellradius(i, step)
    rmin = (i - 1) * step
    return (0.5 * ((rmin + step)^3 + rmin^3))^(1 / 3)
end


text:
"""
    cor(x[, y], ::Type{<:Correlation})

Compute the correlation matrix of a given type.

The possible correlation types are:
  * [`Pearson`](@ref)
  * [`Spearman`](@ref)
  * [`Kendall`](@ref)

# Examples
```jldoctest
julia> x = [-1.62169     0.0158613   0.500375  -0.794381
             2.50689     3.31666    -1.3049     2.16058
             0.495674    0.348621   -0.614451  -0.193579
             2.32149     2.18847    -1.83165    2.08399
            -0.0573697   0.39908     0.270117   0.658458
             0.365239   -0.321493   -1.60223   -0.199998
            -0.55521    -0.898513    0.690267   0.857519
            -0.356979   -1.03724     0.714859  -0.719657
            -3.38438    -1.93058     1.77413   -1.23657
             1.57527     0.836351   -1.13275   -0.277048];

julia> cor(x, Pearson)
4×4 Matrix{Float64}:
  1.0        0.86985   -0.891312   0.767433
  0.86985    1.0       -0.767115   0.817407
 -0.891312  -0.767115   1.0       -0.596762
  0.767433   0.817407  -0.596762   1.0

julia> cor(x, Spearman)
4×4 Matrix{Float64}:
  1.0        0.866667  -0.854545   0.709091
  0.866667   1.0       -0.781818   0.684848
 -0.854545  -0.781818   1.0       -0.612121
  0.709091   0.684848  -0.612121   1.0

julia> cor(x, Kendall)
4×4 Matrix{Float64}:
  1.0        0.733333  -0.688889   0.555556
  0.733333   1.0       -0.688889   0.555556
 -0.688889  -0.688889   1.0       -0.422222
  0.555556   0.555556  -0.422222   1.0
```
"""
function cor end
cor(x,    ::Type{Pearson})  = cor(x)
cor(x, y, ::Type{Pearson})  = cor(x, y)
cor(x,    ::Type{Spearman}) = corspearman(x)
cor(x, y, ::Type{Spearman}) = corspearman(x, y)
cor(x,    ::Type{Kendall})  = corkendall(x)
cor(x, y, ::Type{Kendall})  = corkendall(x, y)

"""
    cor_fast(M::Matrix{S}, T::Type{<:Correlation}=Pearson) where {S<:Real}

Calculate the correlation matrix in parallel.
"""
function cor_fast(M::Matrix{S}, T::Type{<:Correlation}=Pearson) where {S<:Real}
    n, d = size(M)
    C = SharedMatrix{S}(d, d)
    @threads for i in collect(subsets(1:d, Val{2}()))
        C[i...] = cor(view(M, :, i[1]), view(M, :, i[2]), T)
    end
    C .= Symmetric(C, :U)
    C[diagind(C)] .= one(S)
    sdata(C)
end

"""
    cor_convert(X::Matrix{<:Real}, from::Type{<:Correlation}, to::Type{<:Correlation})

Convert from one type of correlation matrix to another.

The role of conversion in this package is typically used from either Spearman or
Kendall to Pearson where the Pearson correlation is used in the generation of
random multivariate normal samples. After converting, the correlation matrix
may not be positive semidefinite, so it is recommended to check using
`LinearAlgebra.isposdef`, and subsequently calling [`cor_nearPD`](@ref).

See also: [`cor_nearPD`](@ref), [`cor_fastPD`](@ref)

The possible correlation types are:
  * [`Pearson`](@ref)
  * [`Spearman`](@ref)
  * [`Kendall`](@ref)

# Examples
```jldoctest
julia> r = [ 1.0       -0.634114   0.551645   0.548993
            -0.634114   1.0       -0.332105  -0.772114
             0.551645  -0.332105   1.0        0.143949
             0.548993  -0.772114   0.143949   1.0];

julia> cor_convert(r, Pearson, Spearman)
4×4 Matrix{Float64}:
  1.0       -0.616168   0.533701   0.531067
 -0.616168   1.0       -0.318613  -0.756979
  0.533701  -0.318613   1.0        0.13758
  0.531067  -0.756979   0.13758    1.0

julia> cor_convert(r, Spearman, Kendall)
4×4 Matrix{Float64}:
  1.0       -0.452063   0.385867    0.383807
 -0.452063   1.0       -0.224941   -0.576435
  0.385867  -0.224941   1.0         0.0962413
  0.383807  -0.576435   0.0962413   1.0

julia> r == cor_convert(r, Pearson, Pearson)
true
```
"""
function cor_convert end
cor_convert(x::Real, from::Type{C},        to::Type{C}) where {C<:Correlation} = x
cor_convert(x::Real, from::Type{Pearson},  to::Type{Spearman}) = cor_constrain(asin(x / 2) * 6 / π)
cor_convert(x::Real, from::Type{Pearson},  to::Type{Kendall})  = cor_constrain(asin(x) * 2 / π)
cor_convert(x::Real, from::Type{Spearman}, to::Type{Pearson})  = cor_constrain(sin(x * π / 6) * 2)
cor_convert(x::Real, from::Type{Spearman}, to::Type{Kendall})  = cor_constrain(asin(sin(x * π / 6) * 2) * 2 / π)
cor_convert(x::Real, from::Type{Kendall},  to::Type{Pearson})  = cor_constrain(sin(x * π / 2))
cor_convert(x::Real, from::Type{Kendall},  to::Type{Spearman}) = cor_constrain(asin(sin(x * π / 2) / 2) * 6 / π)
function cor_convert(X::VecOrMat{<:Real}, from::Type{<:Correlation}, to::Type{<:Correlation})
    cor_constrain(cor_convert.(copy(X), from, to))
end

"""
    cor_constrain!(C::Matrix{<:Real}[, uplo=:U])

Same as [`cor_constrain`](@ref), except that the matrix `C` is updated in place
to save memory.

# Examples
```jldoctest
julia> a = [ 0.802271   0.149801  -1.1072     1.13451
             0.869788  -0.824395   0.38965    0.965936
            -1.45353   -1.29282    0.417233  -0.362526
             0.638291  -0.682503   1.12092   -1.27018];

julia> cor_constrain!(a)

julia> a
4×4 Matrix{Float64}:
  1.0       0.149801  -1.0        1.0
  0.149801  1.0        0.38965    0.965936
 -1.0       0.38965    1.0       -0.362526
  1.0       0.965936  -0.362526   1.0
```
"""
function cor_constrain!(C::Matrix{<:Real}, uplo=:U)
    C .= clampcor.(C)
    C .= Symmetric(C, uplo)
    C[diagind(C)] .= one(eltype(C))
    nothing
end

"""
    cor_constrain(C::Matrix{<:Real}[, uplo=:U])

Constrain a matrix so that its diagonal elements are 1, off-diagonal elements
are bounded between -1 and 1, and a symmetric view of the upper (if `uplo = :U`)
or lower (if `uplo = :L`) triangle is made.

See also: [`cor_constrain!`](@ref)

# Examples
```jldoctest
julia> a = [ 0.802271   0.149801  -1.1072     1.13451
             0.869788  -0.824395   0.38965    0.965936
            -1.45353   -1.29282    0.417233  -0.362526
             0.638291  -0.682503   1.12092   -1.27018];

julia> cor_constrain(a)
4×4 Matrix{Float64}:
  1.0       0.149801  -1.0        1.0
  0.149801  1.0        0.38965    0.965936
 -1.0       0.38965    1.0       -0.362526
  1.0       0.965936  -0.362526   1.0

julia> cor_constrain(a, :L)
4×4 Matrix{Float64}:
  1.0        0.869788  -1.0   0.638291
  0.869788   1.0       -1.0  -0.682503
 -1.0       -1.0        1.0   1.0
  0.638291  -0.682503   1.0   1.0
```
"""
function cor_constrain(C::Matrix{<:Real}, uplo=:U)
    R = copy(C)
    cor_constrain!(R, uplo)
    R
end
cor_constrain(x::Real) = clamp(x, -one(eltype(x)), one(eltype(x)))

"""
    cov2cor(C::Matrix{<:AbstractFloat})

Transform a covariance matrix into a correlation matrix.

# Details

If ``X \\in \\mathbb{R}^{n \\times n}`` is a covariance matrix, then

```math
\\tilde{X} = D^{-1/2} X  D^{-1/2}, \\quad D = \\mathrm{diag(X)}
```

is the associated correlation matrix.
"""
function cov2cor(C::Matrix{<:AbstractFloat})
    D = sqrt(inv(Diagonal(C)))
    cor_constrain(D * C * D)
end

"""
    cov2cor!(C::Matrix{<:AbstractFloat})

Same as [`cov2cor`](@ref), except that the matrix `C` is updated in place
to save memory.
"""
function cov2cor!(C::Matrix{<:AbstractFloat})
    D = sqrt(inv(Diagonal(C)))
    C .= D * C * D
    cor_constrain!(C)
    nothing
end

Base.clamp(R::Matrix{<:Real}, L::Matrix{<:Real}, U::Matrix{<:Real}) = clamp.(R, L, U)

function cor_clamp(R, L, U)
    L2 = copy(L)
    U2 = copy(U)
    L2[diagind(L2)] .= -Inf
    U2[diagind(U2)] .= Inf
    R2 = clamp(R, L2, U2)
    cor_constrain(R2)
end


text:
import MultivariatePolynomials
const MP = MultivariatePolynomials
import DynamicPolynomials: @polyvar
using StaticPolynomials
const SP = StaticPolynomials
using StaticArrays
using Test
using LinearAlgebra

@testset "StaticPolynomials" begin
    include("codegen_tests.jl")
    include("basic_tests.jl")
    include("system_evaluation_tests.jl")
end


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test.avr")
sname = joinpath(dirname(@__FILE__), "..", "data", "tmp", "same.avr")

a, b = read_avr(fname)

write_avr(sname, a, b, 8192)

a2, b2 = read_avr(sname)

@test a == a2
@test b == b2


text:
#
# NOTE: This is not a formulation of any particular problem
# It is only for testing and illustration purposes
#

using JuMP
PMs = PowerModels

@testset "test multinetwork" begin

    function build_mn_data(base_data)
        mp_data = PowerModels.parse_file(base_data)
        return PowerModels.replicate(mp_data, 2)
    end

    function build_mn_data(base_data_1, base_data_2)
        mp_data_1 = PowerModels.parse_file(base_data_1)
        mp_data_2 = PowerModels.parse_file(base_data_2)

        @assert mp_data_1["per_unit"] == mp_data_2["per_unit"]
        @assert mp_data_1["baseMVA"] == mp_data_2["baseMVA"]

        mn_data = Dict{String,Any}(
            "name" => "$(mp_data_1["name"]) + $(mp_data_2["name"])",
            "multinetwork" => true,
            "per_unit" => mp_data_1["per_unit"],
            "baseMVA" => mp_data_1["baseMVA"],
            "nw" => Dict{String,Any}()
        )

        delete!(mp_data_1, "multinetwork")
        delete!(mp_data_1, "per_unit")
        delete!(mp_data_1, "baseMVA")
        mn_data["nw"]["1"] = mp_data_1

        delete!(mp_data_2, "multinetwork")
        delete!(mp_data_2, "per_unit")
        delete!(mp_data_2, "baseMVA")
        mn_data["nw"]["2"] = mp_data_2

        return mn_data
    end

    @testset "idempotent unit transformation" begin
        @testset "5-bus replicate case" begin
            mn_data = build_mn_data("../test/data/case5_dc.m")
            PowerModels.make_mixed_units(mn_data)
            PowerModels.make_per_unit(mn_data)

            @test compare_dict(mn_data, build_mn_data("../test/data/case5_dc.m"))
        end
        @testset "14+24 hybrid case" begin
            mn_data = build_mn_data("../test/data/case14.m", "../test/data/case24.m")
            PowerModels.make_mixed_units(mn_data)
            PowerModels.make_per_unit(mn_data)

            @test compare_dict(mn_data, build_mn_data("../test/data/case14.m", "../test/data/case24.m"))
        end
    end

    function post_mpopf_test(pm::GenericPowerModel)
        for (n, network) in pm.ref[:nw]
            PMs.variable_voltage(pm, n)
            PMs.variable_generation(pm, n)
            PMs.variable_branch_flow(pm, n)
            PMs.variable_dcline_flow(pm, n)

            PMs.constraint_voltage(pm, n)

            for i in ids(pm, n, :ref_buses)
                PMs.constraint_theta_ref(pm, n, i)
            end

            for i in ids(pm, n, :bus)
                PMs.constraint_kcl_shunt(pm, n, i)
            end

            for i in ids(pm, n, :branch)
                PMs.constraint_ohms_yt_from(pm, n, i)
                PMs.constraint_ohms_yt_to(pm, n, i)

                PMs.constraint_voltage_angle_difference(pm, n, i)

                PMs.constraint_thermal_limit_from(pm, n, i)
                PMs.constraint_thermal_limit_to(pm, n, i)
            end

            for i in ids(pm, n, :dcline)
                PMs.constraint_dcline(pm, n, i)
            end
        end

        # cross network constraint, just for illustration purposes
        # designed to be feasible with two copies of case5_asym.m
        t1_pg = var(pm, 1, :pg)
        t2_pg = var(pm, 2, :pg)
        @constraint(pm.model, t1_pg[2] == t2_pg[4])

        PMs.objective_min_fuel_cost(pm, keys(pm.ref[:nw]))
    end

    @testset "2 period 5-bus asymmetric case" begin
        mn_data = build_mn_data("../test/data/case5_asym.m")

        @testset "test ac polar opf" begin
            result = run_generic_model(mn_data, ACPPowerModel, ipopt_solver, post_mpopf_test, multinetwork=true)

            @test result["status"] == :LocalOptimal
            @test isapprox(result["objective"], 35184.2; atol = 1e0)
            @test isapprox(
                result["solution"]["nw"]["1"]["gen"]["2"]["pg"],
                result["solution"]["nw"]["2"]["gen"]["4"]["pg"];
                atol = 1e-3
            )
        end
    end

    @testset "hybrid network case - polar" begin
        mn_data = build_mn_data("../test/data/case14.m", "../test/data/case24.m")

        @testset "test ac polar opf" begin
            result = run_generic_model(mn_data, ACPPowerModel, ipopt_solver, post_mpopf_test, multinetwork=true)

            @test result["status"] == :LocalOptimal
            @test isapprox(result["objective"], 88289.0; atol = 1e0)
            @test isapprox(
                result["solution"]["nw"]["1"]["gen"]["2"]["pg"],
                result["solution"]["nw"]["2"]["gen"]["4"]["pg"];
                atol = 1e-3
            )
        end
    end

    @testset "hybrid network case - rect" begin
        mn_data = build_mn_data("../test/data/case14.m", "../test/data/case24.m")

        @testset "test ac polar opf" begin
            result = run_generic_model(mn_data, ACRPowerModel, ipopt_solver, post_mpopf_test, multinetwork=true)

            @test result["status"] == :LocalOptimal
            @test isapprox(result["objective"], 88289.0; atol = 1e0)
            @test isapprox(
                result["solution"]["nw"]["1"]["gen"]["2"]["pg"],
                result["solution"]["nw"]["2"]["gen"]["4"]["pg"];
                atol = 1e-3
            )
        end
    end

    @testset "5-bus asymmetric case" begin
        # this works, but should throw a warning
        mn_data = build_mn_data("../test/data/case5_asym.m")
        result = run_ac_opf(mn_data, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 17551; atol = 1e0)
    end

    function post_mppf_test(pm::GenericPowerModel)
        for (n, network) in pm.ref[:nw]
            PMs.variable_voltage(pm, n, bounded = false)
            PMs.variable_generation(pm, n, bounded = false)
            PMs.variable_branch_flow(pm, n, bounded = false)
            PMs.variable_dcline_flow(pm, n, bounded = false)

            PMs.constraint_voltage(pm, n)

            for i in ids(pm, n, :ref_buses)
                PMs.constraint_theta_ref(pm, n, i)
                PMs.constraint_voltage_magnitude_setpoint(pm, n, i)
            end

            for (i,bus) in ref(pm, n, :bus)
                PMs.constraint_kcl_shunt(pm, n, i)

                # PV Bus Constraints
                if length(ref(pm, n, :bus_gens, i)) > 0 && !(i in ids(pm, n, :ref_buses))
                    @assert bus["bus_type"] == 2

                    PMs.constraint_voltage_magnitude_setpoint(pm, n, i)
                    for j in ref(pm, n, :bus_gens, i)
                        PMs.constraint_active_gen_setpoint(pm, n, j)
                    end
                end
            end

            for i in ids(pm, n, :branch)
                PMs.constraint_ohms_yt_from(pm, n, i)
                PMs.constraint_ohms_yt_to(pm, n, i)
            end

            for (i,dcline) in ref(pm, n, :dcline)
                PMs.constraint_active_dcline_setpoint(pm, n, i)

                f_bus = ref(pm, :bus)[dcline["f_bus"]]
                if f_bus["bus_type"] == 1
                    PMs.constraint_voltage_magnitude_setpoint(pm, n, f_bus["index"])
                end

                t_bus = ref(pm, :bus)[dcline["t_bus"]]
                if t_bus["bus_type"] == 1
                    PMs.constraint_voltage_magnitude_setpoint(pm, n, t_bus["index"])
                end
            end
        end
    end

    @testset "test solution feedback" begin
        mn_data = build_mn_data("../test/data/case5_asym.m")

        opf_result = run_generic_model(mn_data, ACPPowerModel, ipopt_solver, post_mpopf_test, multinetwork=true)
        @test opf_result["status"] == :LocalOptimal
        @test isapprox(opf_result["objective"], 35184.2; atol = 1e0)
        #@test isapprox(opf_result["objective"], 35533.8; atol = 1e0) # case5_dc (out of date)

        PowerModels.update_data(mn_data, opf_result["solution"])

        pf_result = run_generic_model(mn_data, ACPPowerModel, ipopt_solver, post_mppf_test, multinetwork=true)
        @test pf_result["status"] == :LocalOptimal
        @test isapprox(pf_result["objective"], 0.0; atol = 1e-3)

        for (n, nw_data) in mn_data["nw"]
            #println(n)
            for (i,bus) in nw_data["bus"]
                #println(opf_result["solution"]["nw"][n]["bus"][i]["va"])
                #println(pf_result["solution"]["nw"][n]["bus"][i]["va"])
                #println()

                @test isapprox(opf_result["solution"]["nw"][n]["bus"][i]["va"], pf_result["solution"]["nw"][n]["bus"][i]["va"]; atol = 1e-3)
                @test isapprox(opf_result["solution"]["nw"][n]["bus"][i]["vm"], pf_result["solution"]["nw"][n]["bus"][i]["vm"]; atol = 1e-3)
            end

            for (i,gen) in nw_data["gen"]
                @test isapprox(opf_result["solution"]["nw"][n]["gen"][i]["pg"], pf_result["solution"]["nw"][n]["gen"][i]["pg"]; atol = 1e-3)
                # cannot check this value solution does not appeat to be unique; verify this!
                #@test isapprox(opf_result["solution"]["gen"][i]["qg"], pf_result["solution"]["gen"][i]["qg"]; atol = 1e-3)
            end

            for (i,dcline) in nw_data["dcline"]
                @test isapprox(opf_result["solution"]["nw"][n]["dcline"][i]["pf"], pf_result["solution"]["nw"][n]["dcline"][i]["pf"]; atol = 1e-3)
                @test isapprox(opf_result["solution"]["nw"][n]["dcline"][i]["pt"], pf_result["solution"]["nw"][n]["dcline"][i]["pt"]; atol = 1e-3)
            end
        end

    end

end


text:
"""
    update!(fssrb::PotentialSimulationSetupRB{T, 3, 4, S}, RBT::DataType)::Nothing

Loop over `even` grid points. A point is `even` if the sum of its cartesian indicies (of the not extended grid) is even.
Even points get the red black index (rbi) = 2. ( -> rbpotential[ inds..., rbi ]).
"""
@fastmath function update!( fssrb::PotentialSimulationSetupRB{T, 3, 4, S}, use_nthreads::Int,
                            update_even_points::Val{even_points},
                            depletion_handling::Val{depletion_handling_enabled},
                            bulk_is_ptype::Val{_bulk_is_ptype},
                            is_weighting_potential::Val{_is_weighting_potential})::Nothing where {T, S, even_points, depletion_handling_enabled, _bulk_is_ptype, _is_weighting_potential}
    @inbounds begin
        rb_tar_idx::Int, rb_src_idx::Int = even_points ? (rb_even::Int, rb_odd::Int) : (rb_odd::Int,rb_even::Int)

        gw1::Array{T, 2} = fssrb.geom_weights[1].weights  # r or x
        gw2::Array{T, 2} = fssrb.geom_weights[2].weights  # φ or y
        gw3::Array{T, 2} = fssrb.geom_weights[3].weights  # z or z

        # for idx3 in 2:(size(fssrb.potential, 3) - 1)
        @onthreads 1:use_nthreads for idx3 in workpart(2:(size(fssrb.potential, 3) - 1), 1:use_nthreads, Base.Threads.threadid())
            innerloops!( idx3, rb_tar_idx, rb_src_idx, gw1, gw2, gw3, fssrb, update_even_points, depletion_handling, bulk_is_ptype, is_weighting_potential)
        end
    end
    nothing
end

"""
    innerloops!(  ir::Int, rb_tar_idx::Int, rb_src_idx::Int, gw_r::Array{T, 2}, gw_φ::Array{T, 2}, gw_z::Array{T, 2}, fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical},
                                update_even_points::Val{even_points},
                                depletion_handling::Val{depletion_handling_enabled},
                                bulk_is_ptype::Val{_bulk_is_ptype}  )::Nothing where {T, even_points, depletion_handling_enabled, _bulk_is_ptype}

(Vectorized) inner loop for Cylindrical coordinates. This function does all the work in the field calculation.
"""
@fastmath function innerloops!( ir::Int, rb_tar_idx::Int, rb_src_idx::Int, gw_r::Array{T, 2}, gw_φ::Array{T, 2}, gw_z::Array{T, 2}, fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical},
                                update_even_points::Val{even_points},
                                depletion_handling::Val{depletion_handling_enabled},
                                bulk_is_ptype::Val{_bulk_is_ptype},
                                is_weighting_potential::Val{_is_weighting_potential}  )::Nothing where {T, even_points, depletion_handling_enabled, _bulk_is_ptype, _is_weighting_potential}
    @inbounds begin
        inr::Int = ir - 1

        pwwrr::T               = gw_r[1, inr]
        pwwrl::T               = gw_r[2, inr]
        r_inv_pwΔmpr::T        = gw_r[3, inr]
        Δr_ext_inv_r_pwmprr::T = gw_r[4, inr]
        Δr_ext_inv_l_pwmprl::T = gw_r[5, inr]
        Δmpr_squared::T        = gw_r[6, inr]

        for iφ in 2:(size(fssrb.potential, 2) - 1)
            # φ loop
            inφ::Int = iφ - 1
            pwwφr::T        = gw_φ[1, inφ]
            pwwφl::T        = gw_φ[2, inφ]
            pwΔmpφ::T       = gw_φ[3, inφ]
            Δφ_ext_inv_r::T = gw_φ[4,  iφ]
            Δφ_ext_inv_l::T = gw_φ[4, inφ]

            if inr == 1
                pwwφr = T(0.5)
                pwwφl = T(0.5)
                pwΔmpφ = T(2π)
                Δφ_ext_inv_r = inv(pwΔmpφ)
                Δφ_ext_inv_l = Δφ_ext_inv_r
            end

            rφi_is_even::Bool = iseven(ir + iφ)
            rφi_is_even_t = rφi_is_even ? Val{true}() : Val{false}()

            pwwrr_pwwφr::T = pwwrr * pwwφr
            pwwrr_pwwφl::T = pwwrr * pwwφl
            pwwrl_pwwφr::T = pwwrl * pwwφr
            pwwrl_pwwφl::T = pwwrl * pwwφl

            Δr_ext_inv_r_pwmprr_pwΔmpφ::T = Δr_ext_inv_r_pwmprr * pwΔmpφ
            Δr_ext_inv_l_pwmprl_pwΔmpφ::T = Δr_ext_inv_l_pwmprl * pwΔmpφ
            pwΔmpφ_Δmpr_squared::T = pwΔmpφ * Δmpr_squared
            r_inv_pwΔmpr_Δφ_ext_inv_r::T = r_inv_pwΔmpr * Δφ_ext_inv_r
            r_inv_pwΔmpr_Δφ_ext_inv_l::T = r_inv_pwΔmpr * Δφ_ext_inv_l

            @fastmath @inbounds @simd ivdep for iz in 2:(size(fssrb.potential, 1) - 1)
                inz::Int = nidx(iz, update_even_points, rφi_is_even_t)::Int
                # izr::Int = get_rbidx_right_neighbour(iz, update_even_points, rφi_is_even)::Int # this is somehow slower than the two lines below
                izr::Int = ifelse( rφi_is_even, iz, even_points ? iz - 1 : iz + 1)
                izr += ifelse(even_points, 1, 0)

                pwwzr::T        = gw_z[1, inz]
                pwwzl::T        = gw_z[2, inz]
                pwΔmpz::T       = gw_z[3, inz]
                Δz_ext_inv_r::T = gw_z[4, inz + 1]
                Δz_ext_inv_l::T = gw_z[4, inz]

                ϵ_rrr::T = fssrb.ϵ[  ir,  iφ, inz + 1]
                ϵ_rlr::T = fssrb.ϵ[  ir, inφ, inz + 1]
                ϵ_rrl::T = fssrb.ϵ[  ir,  iφ, inz ]
                ϵ_rll::T = fssrb.ϵ[  ir, inφ, inz ]
                ϵ_lrr::T = fssrb.ϵ[ inr,  iφ, inz + 1]
                ϵ_llr::T = fssrb.ϵ[ inr, inφ, inz + 1]
                ϵ_lrl::T = fssrb.ϵ[ inr,  iφ, inz ]
                ϵ_lll::T = fssrb.ϵ[ inr, inφ, inz ]

                vrr::T = fssrb.potential[     iz,     iφ, ir + 1, rb_src_idx]
                vrl::T = fssrb.potential[     iz,     iφ,    inr, rb_src_idx]
                vφr::T = fssrb.potential[     iz, iφ + 1,     ir, rb_src_idx]
                vφl::T = fssrb.potential[     iz,    inφ,     ir, rb_src_idx]
                vzr::T = fssrb.potential[    izr,     iφ,     ir, rb_src_idx]
                vzl::T = fssrb.potential[izr - 1,     iφ,     ir, rb_src_idx]

                pwwφr_pwwzr::T = pwwφr * pwwzr
                pwwφl_pwwzr::T = pwwφl * pwwzr
                pwwφr_pwwzl::T = pwwφr * pwwzl
                pwwφl_pwwzl::T = pwwφl * pwwzl
                pwwrl_pwwzr::T = pwwrl * pwwzr
                pwwrr_pwwzr::T = pwwrr * pwwzr
                pwwrl_pwwzl::T = pwwrl * pwwzl
                pwwrr_pwwzl::T = pwwrr * pwwzl

                # right weight in r: wrr
                wrr::T = ϵ_rrr * pwwφr_pwwzr
                wrr    = muladd(ϵ_rlr, pwwφl_pwwzr, wrr)
                wrr    = muladd(ϵ_rrl, pwwφr_pwwzl, wrr)
                wrr    = muladd(ϵ_rll, pwwφl_pwwzl, wrr)
                # left weight in r: wrr
                wrl::T = ϵ_lrr * pwwφr_pwwzr
                wrl    = muladd(ϵ_llr, pwwφl_pwwzr, wrl)
                wrl    = muladd(ϵ_lrl, pwwφr_pwwzl, wrl)
                wrl    = muladd(ϵ_lll, pwwφl_pwwzl, wrl)
                # right weight in φ: wφr
                wφr::T = ϵ_lrr * pwwrl_pwwzr
                wφr    = muladd(ϵ_rrr, pwwrr_pwwzr, wφr)
                wφr    = muladd(ϵ_lrl, pwwrl_pwwzl, wφr)
                wφr    = muladd(ϵ_rrl, pwwrr_pwwzl, wφr)
                # left weight in φ: wφl
                wφl::T = ϵ_llr * pwwrl_pwwzr
                wφl    = muladd(ϵ_rlr, pwwrr_pwwzr, wφl)
                wφl    = muladd(ϵ_lll, pwwrl_pwwzl, wφl)
                wφl    = muladd(ϵ_rll, pwwrr_pwwzl, wφl)
                # right weight in z: wzr
                wzr::T = ϵ_rrr * pwwrr_pwwφr
                wzr    = muladd(ϵ_rlr, pwwrr_pwwφl, wzr)
                wzr    = muladd(ϵ_lrr, pwwrl_pwwφr, wzr)
                wzr    = muladd(ϵ_llr, pwwrl_pwwφl, wzr)
                # left weight in z: wzr
                wzl::T = ϵ_rrl * pwwrr_pwwφr
                wzl    = muladd(ϵ_rll, pwwrr_pwwφl, wzl)
                wzl    = muladd(ϵ_lrl, pwwrl_pwwφr, wzl)
                wzl    = muladd(ϵ_lll, pwwrl_pwwφl, wzl)

                wrr *= Δr_ext_inv_r_pwmprr_pwΔmpφ * pwΔmpz
                wrl *= Δr_ext_inv_l_pwmprl_pwΔmpφ * pwΔmpz
                wφr *= r_inv_pwΔmpr_Δφ_ext_inv_r * pwΔmpz
                wφl *= r_inv_pwΔmpr_Δφ_ext_inv_l * pwΔmpz
                wzr *= Δz_ext_inv_r * pwΔmpφ_Δmpr_squared
                wzl *= Δz_ext_inv_l * pwΔmpφ_Δmpr_squared

                new_potential::T = _is_weighting_potential ? 0 : fssrb.ρ[iz, iφ, ir, rb_tar_idx]
                new_potential = muladd( wrr, vrr, new_potential)
                new_potential = muladd( wrl, vrl, new_potential)
                new_potential = muladd( wφr, vφr, new_potential)
                new_potential = muladd( wφl, vφl, new_potential)
                new_potential = muladd( wzr, vzr, new_potential)
                new_potential = muladd( wzl, vzl, new_potential)

                new_potential *= fssrb.volume_weights[iz, iφ, ir, rb_tar_idx]

                old_potential::T = fssrb.potential[iz, iφ, ir, rb_tar_idx]

                new_potential -= old_potential
                new_potential = muladd(new_potential, fssrb.sor_const[inr], old_potential)

                if depletion_handling_enabled
                    if inr == 1 vrl = vrr end
                    if _bulk_is_ptype # p-type detectors
                        if new_potential < fssrb.minimum_applied_potential
                            # new_potential = fssrb.minimum_applied_potential
                            new_potential -= fssrb.ρ[iz, iφ, ir, rb_tar_idx] * fssrb.volume_weights[iz, iφ, ir, rb_tar_idx] * fssrb.sor_const[inr]
                            if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                        else
                            vmin::T = ifelse( vrr <  vrl, vrr,  vrl)
                            vmin    = ifelse( vφr < vmin, vφr, vmin)
                            vmin    = ifelse( vφl < vmin, vφl, vmin)
                            vmin    = ifelse( vzr < vmin, vzr, vmin)
                            vmin    = ifelse( vzl < vmin, vzl, vmin)
                            if new_potential <= vmin # bubble point
                                new_potential -= fssrb.ρ[iz, iφ, ir, rb_tar_idx] * fssrb.volume_weights[iz, iφ, ir, rb_tar_idx] * fssrb.sor_const[inr]
                                if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                            else # normal point
                                if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit > 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] -= undepleted_bit end # unmark this point
                            end
                        end
                    else # n-type detectors
                        if new_potential > fssrb.maximum_applied_potential
                            # new_potential = fssrb.maximum_applied_potential
                            new_potential -= fssrb.ρ[iz, iφ, ir, rb_tar_idx] * fssrb.volume_weights[iz, iφ, ir, rb_tar_idx] * fssrb.sor_const[inr]
                            if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                        else
                            vmax::T = ifelse( vrr >  vrl, vrr,  vrl)
                            vmax    = ifelse( vφr > vmax, vφr, vmax)
                            vmax    = ifelse( vφl > vmax, vφl, vmax)
                            vmax    = ifelse( vzr > vmax, vzr, vmax)
                            vmax    = ifelse( vzl > vmax, vzl, vmax)
                            if new_potential >= vmax # bubble point
                                new_potential -= fssrb.ρ[iz, iφ, ir, rb_tar_idx] * fssrb.volume_weights[iz, iφ, ir, rb_tar_idx] * fssrb.sor_const[inr]
                                if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                            else # normal point -> unmark
                                if (fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & undepleted_bit > 0) fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] -= undepleted_bit end # unmark this point
                            end
                        end
                    end
                end

                fssrb.potential[iz, iφ, ir, rb_tar_idx]::T = ifelse(fssrb.pointtypes[iz, iφ, ir, rb_tar_idx] & update_bit > 0, new_potential, old_potential)
            end # z loop
        end # φ loop
    end # inbounds
end

function update!(   fssrb::PotentialSimulationSetupRB{T}; use_nthreads::Int = Base.Threads.nthreads(),
                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(), only2d::Val{only_2d} = Val{false}(),
                    is_weighting_potential::Val{_is_weighting_potential} = Val{false}())::Nothing where {T, depletion_handling_enabled, only_2d, _is_weighting_potential}
    update!(fssrb, use_nthreads, Val{true}(), depletion_handling, Val{fssrb.bulk_is_ptype}(), is_weighting_potential)
    apply_boundary_conditions!(fssrb, Val{true}(), only2d)
    update!(fssrb, use_nthreads, Val{false}(), depletion_handling, Val{fssrb.bulk_is_ptype}(), is_weighting_potential)
    apply_boundary_conditions!(fssrb, Val{false}(), only2d)
    nothing
end


text:
using NetMSA
using Test

@testset "NetMSA.jl" begin
    # Write your tests here.
    S1 = "abcbcdem";
    S2 = "acbcfg";
    S3 = "abchimn";
    S4 = "abcbcjkm";

    L = [S1, S2, S3, S4];
    M = [
        ['a' 'a' 'a' 'a']
        ['b' 'c' 'b' 'b']
        ['c' 'b' 'c' 'c']
        ['b' 'c' 'h' 'b']
        ['c' 'f' 'i' 'c']
        ['d' 'g' 'm' 'j']
        ['e' missing 'n' 'k']
        ['m' missing missing 'm']
        ]

    @test isequal(NetMSA.createPeerMatrix(L), M);

    @test NetMSA.getposition('b', 2, M).indexes == [1, 3, 4]

    @test NetMSA.mostfrequent(M[1, :]) == (4, 'a');

    @test NetMSA.mostfrequent(M[2, :]) == (3, 'b');

    @test NetMSA.full(M[1, :]) == true;

    @test NetMSA.aligned(M[2, :]) == false;

    @test NetMSA.aligned(M[8, :]) == true;

    @test NetMSA.weight(M[1, :], w1=0.25, w2=0.5, w3=1.0) == 1.0;

    @test NetMSA.weight(M[2, :], w1=0.25, w2=0.5, w3=1.0) == 0.1875;

    @test NetMSA.weight(M[4, :], w1=0.25, w2=0.5, w3=1.0) == 0.125;

    @test NetMSA.weight(M[6, :], w1=0.25, w2=0.5, w3=1.0) == 0.0;

    @test NetMSA.weight(M[8, :], w1=0.25, w2=0.5, w3=1.0) == 0.25;

    @test sum(NetMSA.weight.(eachrow(M[2:end, :]))) == 0.875;

    @test NetMSA.objective(M, 2) == 2.625;

    @test_throws ArgumentError NetMSA.objective(M, 2, endindex=9)

    @test length(NetMSA.createswarm(2, M)) == 2;

    @test length(NetMSA.createswarm(8, M)) == 1;

    p = NetMSA.Particle('b', NetMSA.getposition('b', 2, M));

    @test NetMSA.stopcriteria(p, 3, M) == true;

    @test NetMSA.criteria3(p, 3, M) == true;

    @test NetMSA.criteria2(p) == false;

    p.updated = 6;

    @test NetMSA.criteria2(p) == true;

    @test NetMSA.stopcriteria(p, 2, M) == true;

    newrow = [missing missing missing missing];

    @test isequal(NetMSA.remove_missing_rows(vcat(M, newrow)), M);

    N = [
        ['a' 'a' 'a' 'a']
        ['-' 'c' '-' '-']
        ['b' 'b' 'b' 'b']
        ['c' 'c' 'c' 'c']
        ['b' 'f' 'h' 'b']
        ['c' 'g' 'i' 'c']
        ['d' missing 'm' 'j']
        ['e' missing 'n' 'k']
        ['m' missing missing 'm']
        ]

    @test isequal(NetMSA.flydown(p, M), N);

    @test NetMSA.rowalignment(2, M).bestvalue == 9.0;

    @test NetMSA.rowalignment(2, M).best.row == 3;

    aligned_M = [
        ['a' 'a' 'a' 'a']
        ['b' '-' 'b' 'b']
        ['c' 'c' 'c' 'c']
        ['b' 'b' '-' 'b']
        ['c' 'c' '-' 'c']
        ['d' 'f' 'h' 'j']
        ['e' 'g' 'i' 'k']
        ['m' '-' 'm' 'm']
        ['-' '-' 'n' '-']
        ]

    @test isequal(NetMSA.matrixalignment(M), aligned_M);

end


text:
#=
    Optimal Control (feed-forward) approach for a mechanical oscillator with two masses, springs and attenuators

    Author: Stephan Scholz
    Year: 2021
=#

using JuMP,Ipopt, LinearAlgebra

#model initialization
mod = Model( optimizer_with_attributes( Ipopt.Optimizer , "max_iter" => 1000))

const d₁ = 0.5; # Damping
const d₂ = 0.5;
const c₁ = 0.1; # Spring constant
const c₂ = 0.1;
const m₁ = 1.0; # Mass
const m₂ = 2.0;

# System matrix
A = [0 0 1 0; 0 0 0 1; (-1/m₁)*(c₁ + c₂) c₂/m₁ (-1/m₁)*(d₁ + d₂) d₂/m₁; c₂/m₂ -c₂/m₂ d₂/m₂ -d₂/m₂]

B = [0 0; 0 0; 1/m₁ 0; 0 1/m₂]  # Input matrix
C = [1 0 0 0; 0 1 0 0]          # Output matrix

Nx = size(A)[1]
Nu = size(B)[2]

# Input constraints
umin = -10.0;
umax = 10.0;

K = 200; # Number of time steps
x0 = zeros(4); # Initial conditions; also possible: x0 = 3*rand(4);

∆T = 0.1
Tf = K * ∆T; # Final time

Q = 100 * Diagonal(ones(Nx))    # Weighing matrix for states / error: e' Q e or x' Q x
R = Diagonal(ones(Nu))          # Weighing matrix for inputs: u' R u

ref = [2.0, 4.0, 0.0, 0.0]      # Reference

@variable(mod, umin <= u[1:2, 1:K-1] <= umax) # Control signals

@variable(mod, x[1:Nx, 1:K])            # States x
@constraint(mod, x[1:Nx, 1] .== x0 )    # Initial values

@variable(mod, e[1:Nx, 1:K] )           # Error between reference and states
@constraint(mod,  e[1:Nx, 1:K] .== ref[1:Nx] - x[1:Nx, 1:K])

for k = 1: K-1
    @constraint(mod, x[:, k+1] .== x[:, k] + ∆T *(A * x[:,k] + B * u[:,k] )) # Discrete-time system dynamics
end

J = @NLexpression(mod, 0.5 * ∆T * sum( Q[1,1]*e[1,k]^2 + Q[2,2]*e[2,k]^2 + Q[3,3]*e[3,k]^2 + Q[4,4]*e[4,k]^2 + R[1,1]*u[1,i]^2 + R[2,2]*u[2,i] for k=1:K, i=1:K-1 ))
@NLobjective( mod, Min, J)
optimize!(mod)

# Optimization results
x_sol = JuMP.value.(x)
u_sol = JuMP.value.(u)

using Plots
tspan = 0 : ∆T : (K-1)*∆T
plot(tspan, x_sol', label=["x1" "x2" "x3" "x4"], title="States")
plot(tspan[1:end-1], u_sol', label=["u1" "u2"], title="Control signals")

# Rebuild solution as ODE with feed-forward input: u = optimization result

# Definition of ODE
function mech_oscillator(dx,x,p,t)
    n = floor(Int64,(t/Tf)*K) + 1
    u = u_sol2[:,n] # [u1; u2]
    dx .= A*x + B*u # Right-hand side of ODE
end

using DifferentialEquations

u_sol2 = hcat(u_sol, zeros(2,1)) # Input signal found by optimizer
trange = (0, (K-1)*∆T)

prob = ODEProblem(mech_oscillator, x0, trange)
sol  = DifferentialEquations.solve(prob, Euler(), dt = ∆T)
plot(sol, labels=["x1" "x2" "x3" "x4"], title="Optimal Control of a mechanical oscillator")


text:
using Push
using Base.Test

# Load the base configuration used by all tests.
cfg_path = joinpath(dirname(@__FILE__), "configuration/boolean.cfg")
cfg = Push.load_configuration(cfg_path)

# EQUALS.
s = Push.run("(FALSE FALSE BOOLEAN.=)", cfg)
@test s.boolean == [true]
s = Push.run("(FALSE TRUE BOOLEAN.=)", cfg)
@test s.boolean == [false]

# NOT.
s = Push.run("(FALSE BOOLEAN.NOT)", cfg)
@test s.boolean == [true]
s = Push.run("(TRUE BOOLEAN.NOT)", cfg)
@test s.boolean == [false]

# AND.
s = Push.run("(FALSE FALSE BOOLEAN.AND)", cfg)
@test s.boolean == [false]
s = Push.run("(FALSE TRUE BOOLEAN.AND)", cfg)
@test s.boolean == [false]
s = Push.run("(TRUE FALSE BOOLEAN.AND)", cfg)
@test s.boolean == [false]
s = Push.run("(TRUE TRUE BOOLEAN.AND)", cfg)
@test s.boolean == [true]

# OR.
s = Push.run("(FALSE FALSE BOOLEAN.OR)", cfg)
@test s.boolean == [false]
s = Push.run("(FALSE TRUE BOOLEAN.OR)", cfg)
@test s.boolean == [true]
s = Push.run("(TRUE FALSE BOOLEAN.OR)", cfg)
@test s.boolean == [true]
s = Push.run("(TRUE TRUE BOOLEAN.OR)", cfg)
@test s.boolean == [true]

# DUP.
s = Push.run("(FALSE BOOLEAN.DUP)", cfg)
@test length(s.boolean) == 2

# FLUSH.
s = Push.run("(FALSE TRUE FALSE FALSE BOOLEAN.FLUSH)", cfg)
@test isempty(s.boolean)

# FROM FLOAT.
s = Push.run("(0.0 BOOLEAN.FROMFLOAT)", cfg)
@test isempty(s.float) && s.boolean == [false]
s = Push.run("(1.871 BOOLEAN.FROMFLOAT)", cfg)
@test isempty(s.float) && s.boolean == [true]

# FROM INT.
s = Push.run("(0 BOOLEAN.FROMINT)", cfg)
@test isempty(s.integer) && s.boolean == [false]
s = Push.run("(32 BOOLEAN.FROMINT)", cfg)
@test isempty(s.integer) && s.boolean == [true]

# POP.
s = Push.run("(TRUE FALSE BOOLEAN.POP)", cfg)
@test s.boolean == [true]

# ROT.
s = Push.run("(FALSE TRUE TRUE BOOLEAN.ROT)", cfg)
@test s.boolean == [true, false, true]

# SWAP.
s = Push.run("(FALSE TRUE BOOLEAN.SWAP)", cfg)
@test s.boolean == [true, false]

# STACKDEPTH.
s = Push.run("(TRUE FALSE FALSE FALSE BOOLEAN.STACKDEPTH)", cfg)
@test length(s.boolean) == 4 && s.integer == [4]

# SHOVE.
s = Push.run("(3 FALSE FALSE FALSE TRUE BOOLEAN.SHOVE)", cfg)
@test isempty(s.integer) && s.boolean == [true, false, false, false]

# YANK.
s = Push.run("(3 TRUE FALSE FALSE FALSE BOOLEAN.YANK)", cfg)
@test isempty(s.integer) && s.boolean == [false, false, false, true]

# YANKDUP.
s = Push.run("(3 TRUE FALSE FALSE FALSE BOOLEAN.YANKDUP)", cfg)
@test isempty(s.integer) && s.boolean == [true, false, false, false, true]

# RAND.
s = Push.run("(BOOLEAN.RAND)", cfg)
@test length(s.boolean) == 1

# DEFINE.
s = Push.run("(X BOOLEAN.DEFINE X)", cfg)
@test isempty(s.boolean) && s.name == [:X, :X]
s = Push.run("(X TRUE BOOLEAN.DEFINE X)", cfg)
@test s.boolean == [true]


text:
# This is an impelementations of the bi-lanczos algorithms, producing an
# left and right decomposition of the matrix that is orthogonal, and transform
# the matrix to a triangular one.
# Major references are from Professor GreenBaum's book and my own derivations.

mutable struct IterativeBiLanczos

    A::Function                 # the A matrix.
    AH::Function                # The A hermitian matrix
    T::Dict{Tuple, Number}      # The tridiagonal matrix.
    V::Dict{Int64, Union{Number, AbstractVector}}   # The right krylov subspace.
    W::Dict{Int64, Union{Number, AbstractVector}}   # The left krylov subspace.

    tensor_shape::Tuple
    Av::AbstractVector          # Computed results
    AHw::AbstractVector

    itr::Int64                  # Iterations Counts
    """
        The initializations needs:
            * A linear operator: A in the form of a function, acts on tensor.
            * 2 initialization vector r1, r2, for the krylov subspace of A and A^H.
    """
    function IterativeBiLanczos(
        A::Function,
        AH::Function,
        r1::AbstractArray,
        r2::AbstractArray
    )
        this = new()
        @assert dot(r1, r2) != 0 "The 2 initialization vector for bi-lanczos cannot be orthogobal to each other. "
        r1 /= norm(r1)
        r2 /= dot(r2, r1)
        this.A = A; this.AH = AH
        Av = this.A(r1); AHw = this.AH(r2)   # try left and right apply the vector to the linear operator.
        @assert size(Av) == size(AHw) "The shape of left and right krylov ortho tensor should be in the same shape. "
        this.Av = Av; this.AHw = AHw
        # Initialization of T, tridiag.
        this.T = Dict{Tuple, typeof(Av[1])}()
        this.T[0, 1] = 0             # Gamma zero
        this.T[1, 0] = 0             # beta zero
        this.T[1, 1] = dot(Av, r2)   # alpha one
        # Initializaiton of V, W ortho subspace.
        this.V = Dict{Int64, Union{AbstractVector, Number}}();
        this.W = Dict{Int64, Union{AbstractVector, Number}}()
        this.V[0] = 0; this.W[0] = 0
        this.V[1] = r1; this.W[1] = r2
        # Misc Parameters.
        this.tensor_shape = size(Av)
        this.itr = 0
        return this
    end

    function IterativeBiLanczos(A::AbstractMatrix{T}) where {T<:Number}
        AFunc = (x) -> A*x
        AHFunc = (x) -> A'*x
        r1 = rand(T, size(A, 1))
        r2 = rand(T, size(A', 1))
        return IterativeBiLanczos(AFunc, AHFunc, r1, r2)
    end
end

function RightApply(this::IterativeBiLanczos, v::AbstractVector)
    return reshape(this.A(reshape(v, this.tensor_shape)), :)
end

function LeftApply(this::IterativeBiLanczos, v::AbstractVector)
    return reshape(this.AH(reshape(v, this.tensor_shape)), :)
end

function GetTMatrix(this::IterativeBiLanczos, start::Int64 = 1, stop = nothing)
    stop = stop === nothing ? this.itr : stop
    T = this.T
    valuesX = Vector{Int64}()
    valuesY = Vector{Int64}()
    values = Vector{Number}()
    for Idx in start: stop
        β = T[Idx, Idx + 1]
        push!(valuesX, Idx); push!(valuesY, Idx + 1); push!(values, β)
        α = T[Idx, Idx]
        push!(valuesX, Idx); push!(valuesY, Idx); push!(values, α)
        γ = T[Idx + 1, Idx]
        push!(valuesX, Idx + 1); push!(valuesY, Idx); push!(values, γ)
    end
    # Only on element in the matrix I think
    if this.itr == 0
        return T[1, 1]
    end
    push!(valuesX, stop + 1)
    push!(valuesY, stop + 1)
    push!(values, T[stop + 1, stop + 1])
    return sparse(valuesX, valuesY, values)
end

function GetWMatrix(this::IterativeBiLanczos)
    return GetWMatrix(this, 1, this.itr)
end

function GetWMatrix(this::IterativeBiLanczos, start::Int64, stop::Int64)
    return hcat([this.W[Idx] for Idx in start: stop + 1]...)
end

function GetVMatrix(this::IterativeBiLanczos)
    return GetVMatrix(this, 1, this.itr)
end

function GetVMatrix(this::IterativeBiLanczos, start::Int64, stop::Int64)
    return hcat([this.V[Idx] for Idx in start: stop + 1]...)
end

"""
    Returns the matrix, V, T, and W.
    A = VTW^T.
"""
function Decomposition(this::IterativeBiLanczos)
    return nothing
end

"""
    Performs one step of the iterative bi-lanczos algorithm and
    return the v, w matrix.
"""
function (this::IterativeBiLanczos)()
    this.itr += 1
    j = this.itr

    Av = this.Av; AHw = this.AHw
    V = this.V; W = this.W
    T = this.T
    α = T[j, j]
    ṽ = similar(Av)
    w̃ = similar(AHw)
    @. ṽ =  Av - α*V[j] - T[j - 1, j]*V[j - 1]
    @. w̃ = AHw - conj(α)*W[j] - T[j, j - 1]*W[j - 1]
    # Invariant space might happen here, but due to numerical stability, we don't know

    γ = norm(ṽ)
    v = ṽ/γ  # v_{j + 1} set up here.
    β = dot(v, w̃)
    if β == 0
        error("BiLanczos serious break down, both v, w are non-zero but inner product is zero. ")
    end
    w = w̃/conj(β)

    this.T[j, j + 1] = β
    this.T[j + 1, j] = γ
    this.W[j + 1] = w
    this.V[j + 1] = v

    # update variables for the next run.
    this.Av = RightApply(this, V[j + 1])
    this.AHw = LeftApply(this, W[j + 1])
    this.T[j + 1, j + 1] = dot(RightApply(this, v), w)  # Update alpha.

    return v, w
end

"""
    Performs the iterations for j step, and returned the T, W, and V
    matrix generated during the running of the algorithms.
"""
function (this::IterativeBiLanczos)(j::Int64)
    start = this.itr + 1
    for _ in 1: j
        this()
    end
    return GetVMatrix(this, start, this.itr),
        GetWMatrix(this, start, this.itr),
        GetTMatrix(this, start)
end

# ==============================================================================
# Basic Testing


text:
# Autogenerated wrapper script for StarPU_jll for x86_64-unknown-freebsd-libgfortran5
export libstarpu

using Hwloc_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("StarPU")
JLLWrappers.@declare_library_product(libstarpu, "libstarpu-1.3.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Hwloc_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libstarpu,
        "lib/libstarpu-1.3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Solvers

using LinearAlgebra
using NLSolvers,Roots
using PositiveFactorizations
using DiffResults, ForwardDiff
using StaticArrays
import NaNMath
    function cholesky_linsolve(d,B,∇f)
        cholesky!(Positive, B)
        Bchol = Cholesky(B,'L',0)
        d .=  Bchol\∇f
    end

    function cholesky_linsolve(B,∇f)
        Bchol =cholesky(Positive, B)
        Bchol\∇f
    end

    #
    Base.summary(::NLSolvers.Newton{<:Direct, typeof(cholesky_linsolve)}) = "Newton's method with Cholesky linsolve"
    CholeskyNewton() = NLSolvers.Newton(linsolve=cholesky_linsolve)
    export CholeskyNewton

    """
        det_22(a,b,c,d)

    Calculates `a*b - c*d` without less rounding error than doing it naively
    """
    function det_22(a,b,c,d)
        t = c*d
        e = muladd(c,d,-t) #cd - cd
        f = muladd(a,b,-t) #ab - cd
        return f-e  #ab - cd + cd - cd
    end

    function solve_cubic_eq(poly::AbstractVector{T}) where {T<:Real}
        # copied from PolynomialRoots.jl, adapted to be AD friendly
        # Cubic equation solver for complex polynomial (degree=3)
        # http://en.wikipedia.org/wiki/Cubic_function   Lagrange's method
        # poly = (a,b,c,d) that represents ax3 + bx2 + cx2 + d

        _1 = one(T)
        third = _1/3
        a1  =  complex(one(T) / poly[4])
        E1  = -poly[3]*a1
        E2  =  poly[2]*a1
        E3  = -poly[1]*a1
        s0  =  E1
        E12 =  E1*E1
        A   =  det_22(2*E1,E12,9*E1,E2) + 27*E3
        #A   =  2*E1*E12 - 9*E1*E2 + 27*E3 # = s1^3 + s2^3
        B = det_22(E1,E1,3,E2)
        #B   =  E12 - 3*E2                 # = s1 s2
        # quadratic equation: z^2 - Az + B^3=0  where roots are equal to s1^3 and s2^3
        Δ2 = det_22(A,A,4*B*B,B)
        Δ = sqrt(Δ2)
        #Δ = (A*A - 4*B*B*B)^0.5
        if real(conj(A)*Δ)>=0 # scalar product to decide the sign yielding bigger magnitude
            s1 = exp(log(0.5 * (A + Δ)) * third)
        else
            s1 = exp(log(0.5 * (A - Δ)) * third)
        end
        if s1 == 0
            s2 = s1
        else
            s2 = B / s1
        end
        zeta1 = complex(-0.5, sqrt(T(3.0))*0.5)
        zeta2 = conj(zeta1)
        return (third*(s0 + s1 + s2), third*(s0 + s1*zeta2 + s2*zeta1), third*(s0 + s1*zeta1 + s2*zeta2))
    end

    """
        roots3(pol)
    solves a cubic equation of the form pol[1] + pol[2]*x + pol[3]*x^2 + pol[4]*x^3
    """
    function roots3(pol)
        return SVector(solve_cubic_eq(pol))
    end

    function roots3(a,b,c,d)
        x = SVector(a,b,c,d)
        return roots3(x)
    end

    """
    function x_sol(res::NLSolvers.ConvergenceInfo)

        Returns the scalar or vector x that solves the system of equations or is the minimizer of an optimization procedure.
    """

    function x_sol(res::NLSolvers.ConvergenceInfo{<:NLSolvers.LineSearch, <:Any, <:NLSolvers.NEqOptions})
        res.info.solution
    end

    function x_sol(res::NLSolvers.ConvergenceInfo{<:NLSolvers.TrustRegion, <:Any, <:NLSolvers.NEqOptions})
        return res.info.zero
    end

    function x_sol(res::NLSolvers.ConvergenceInfo{<:Any, <:Any, <:NLSolvers.OptimizationOptions})
        return res.info.minimizer
    end

    @inline log(x::Number) = Base.log(x)
    @inline log(x::Complex) = Base.log(x)
    @inline log(x::Real) = ifelse(x>=zero(x),Base.log(x),zero(x)/zero(x))
    @inline log(x::Float64) = NaNMath.log(x)
    @inline log(x::Float32) = NaNMath.log(x)
    @inline log(x::Int8) = NaNMath.log(float(x))
    @inline log(x::Int16) = NaNMath.log(float(x))
    @inline log(x::Int32) = NaNMath.log(float(x))
    @inline log(x::Int64) = NaNMath.log(float(x))

    include("ADNewton.jl")
    include("nlsolve.jl")
    include("fixpoint/fixpoint.jl")
    include("optimize.jl")

end # module


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"6.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using FinEtools
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionCircle, CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio, default_layout_3d, plot_from_json
using PlotlyJS
using JSON

L = 42
xyz = [0 0 0;
0 L/4 L*1/4;
L/4 L/4 L*2/4;
L/4 0 L*3/4;
0 0 L]
nL = 20

cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
fens, fes = frame_member(xyz, nL, cs)
plots = cat(plot_nodes(fens),
    plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4),
    dims = 1)
layout = default_layout_3d(autosize=true)
# layout = default_layout_3d(width=400, height=400)
layout[:showLegend] = true

config  = PlotConfig(plotlyServerURL="https://chart-studio.plotly.com", showLink=true)
pl = plot(plots, layout; config = config)

display(pl)
true


text:
using NumericalMethodsforEngineers, IterativeSolvers

a = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.]

c = a \ b

println("\nSolution Vector: $c")