text:
# --- extract values from Val() objects --- #

val(x) = x
val(::Val{x}) where x = x

# --- automatic differentiation --- #

# out-of-place
function autodiff_jacobian_func(f, ivar)
    (args...) -> begin
        x = args[ivar]
        fy = (x) -> f(Base.setindex(args, x, ivar)...)
        ForwardDiff.jacobian(fy, x)
    end
end

# in-place
function autodiff_jacobian_func(f, ny, ivar)
    (J, args...) -> begin
        x = args[ivar]
        y = similar(J, val(ny))
        fy = (y, x) -> f(y, Base.setindex(args, x, ivar)...)
        ForwardDiff.jacobian!(J, fy, y, x)
    end
end

# out-of-place
function autodiff_derivative_func(f, ivar)
    (args...) -> begin
        x = args[ivar]
        fy = (x) -> f(Base.setindex(args, x, ivar)...)
        ForwardDiff.derivative(fy, x)
    end
end

# in-place
function autodiff_derivative_func(f, ny, ivar)
    (dy, args...) -> begin
        x = args[ivar]
        y = similar(dy, val(ny))
        fy = (y, x) -> f(y, Base.setindex(args, x, ivar)...)
        ForwardDiff.derivative!(dy, fy, y, x)
    end
end

# --- block diagonal --- #

function block_diagonal(Jii)

    J = zeros(sum(size.(Jii, 1)), sum(size.(Jii, 2)))

    i = j = 1
    for iblock = 1:N
        nx, ny = size(Jii[iblock])
        ix, iy = i:i+nx-1, j:j+ny-1
        J[ix, iy] = Jii[iblock]
        i += nx; j += ny
    end

    return J
end

block_diagonal(Jii::NTuple{N,StaticMatrix}) where N = static_block_diagonal(Jii)

@generated function static_block_diagonal(Jii::NTuple{N,StaticMatrix}) where N

    # assemble first row of the jacobian matrix
    expr = quote
        J = Jii[1]
    end
    for j = 2:N
        if 1 == j
            expr = quote
                $expr
                Jij = Jii[1]
            end
        else
            expr = quote
                $expr
                nx = Size(Jii[1])[1]
                ny = Size(Jii[$j])[2]
                sz = Size(nx, ny)
                Jij = zero(similar_type(Jii[1], sz))
            end
        end
        expr = quote
            $expr
            J = hcat(J, Jij)
        end
    end

    # add remaining rows of the jacobian matrix
    for i = 2:N
        if i == 1
            expr = quote
                $expr
                Ji = Jii[$i]
            end
        else
            expr = quote
                $expr
                nx = Size(Jii[$i])[1]
                ny = Size(Jii[1])[2]
                sz = Size(nx, ny)
                Ji = zero(similar_type(Jii[$i], sz))
            end
        end
        for j = 2:N
            if i == j
                expr = quote
                    $expr
                    Jij = Jii[$i]
                end
            else
                expr = quote
                    $expr
                    nx = Size(Jii[$i])[1]
                    ny = Size(Jii[$j])[2]
                    sz = Size(nx, ny)
                    Jij = zero(similar_type(Jii[$i], sz))
                end
            end
            expr = quote
                $expr
                Ji = hcat(Ji, Jij)
            end
        end
        expr = quote
            $expr
            J = vcat(J, Ji)
        end
    end

    return expr
end

# --- eigenmode correlation --- #

function correlate_eigenmodes(U, M, V; tracked_modes=1:size(U,1), rtol=1e-3, atol=1e-3)

    # get number of (not necessarily unique) eigenvalues
    nev = size(U,1)

    # construct the correlation matrix
    C = U*M*V

    # get the row permutation that puts maximum values on the diagonals
    rval = Vector{real(eltype(C))}(undef, nev) # magnitude of values in chosen row
    cval = Vector{real(eltype(C))}(undef, nev) # magnitude of values in chosen column
    perm = zeros(Int, nev) # permutation
    corruption = Vector{real(eltype(C))}(undef, nev) # corruption index for each correlation

    # tracked modes
    for i in tracked_modes

        # get correlation matrix entry magnitudes for this row
        for j = 1:nev
            rval[j] = abs(C[i,j])
        end

        # rank modes correlations from best to worst for this row
        ranked_modes = sortperm(rval, rev=true)

        # choose the best unassigned mode correlation for this row
        j1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]

        # get correlation matrix entry magnitudes for the chosen column
        for k = 1:nev
            cval[k] = abs(C[k,j1])
        end

        # find approximate multiplicity of left and/or right eigenvector
        multiplicity = max(
            count(x -> isapprox(x, rval[j1]; rtol, atol), rval),
            count(x -> isapprox(x, cval[i]; rtol, atol), cval)
        )

        if multiplicity < nev
            # determine n+1 best unassigned fit, where n is the approximate multiplicity
            if j1 in ranked_modes[1:multiplicity]
                j2 = ranked_modes[multiplicity+1]
            else
                j2 = ranked_modes[1]
            end

            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = rval[j2]/rval[j1]
        else
            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = 0.0
        end
    end

    # remaining modes
    for i = 1:nev

        if !(i in tracked_modes)

            # get correlation matrix entry magnitudes for this row
            for j = 1:nev
                rval[j] = abs(C[i,j])
            end

            # rank modes correlations from best to worst for this row
            ranked_modes = sortperm(rval, rev=true)

            # choose the best unassigned mode correlation for this row
            j1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]

            # get correlation matrix entry magnitudes for the chosen column
            for k = 1:nev
                cval[k] = abs(C[k,j1])
            end

            # find approximate multiplicity of left and/or right eigenvector
            multiplicity = max(
                count(x -> isapprox(x, rval[j1]; rtol, atol), rval),
                count(x -> isapprox(x, cval[i]; rtol, atol), cval)
            )

            # determine n+1 best unassigned fit, where n is the approximate multiplicity
            if j1 in ranked_modes[1:multiplicity]
                j2 = ranked_modes[multiplicity+1]
            else
                j2 = ranked_modes[1]
            end

            # assign best eigenmode fit, create corruption index
            perm[i] = j1
            corruption[i] = rval[j2]/rval[j1]
        end

    end

    return perm, corruption
end


text:
function rkck(y::Vector{Float64},dy::Vector{Float64},n::Int64,t::Float64,
    h::Float64,yout::Vector{Float64},yerr::Vector{Float64},
    k::Float64,mparams::WK3.ModelParams)
    # function evaluation weights
    a2=0.2;a3=0.3;a4=0.6;a5=1.0;a6=0.875;b21=0.2;b31=3.0/40.0;b32=9.0/40.0;
    b41=0.3;b42=-0.9;b43=1.2;b51=-11.0/54.0;b52=2.5;b53=-70.0/27.0;b54=35.0/27.0;
    b61=1631/55296;b62=175/512;b63=575/13824;b64=44275/110592;b65=253/4096;c1=37/378;
    c3=250/621;c4=125/594;c6=512/1771;dc5=-277/14336;dc1=c1-2825/27648;dc3=c3-18575/48384;
    dc4=c4-13525/55296;dc6=c6-0.25;

    # allocators for intermediate derivatives and solution
    ak2=zeros(n);
    ak3=zeros(n);
    ak4=zeros(n);
    ak5=zeros(n);
    ak6=zeros(n);
    ytemp=zeros(n);

    # first step
    ytemp .= y .+ h.*b21.*dy;
    # second step
    wk3odes(t+a2*h,ytemp,ak2,mparams,k);
    ytemp .= y .+ h.*(b31.*dy.+b32.*ak2);
    # third step
    wk3odes(t+a3*h,ytemp,ak3,mparams,k);
    ytemp .= y .+ h.*(b41.*dy.+b42.*ak2.+b43.*ak3);
    # fourth step
    wk3odes(t+a4*h,ytemp,ak4,mparams,k);
    ytemp .= y .+ h.*(b51.*dy.+b52.*ak2.+b53.*ak3.+b54.*ak4);
    # fifth step
    wk3odes(t+a5*h,ytemp,ak5,mparams,k);
    ytemp .= y .+ h.*(b61.*dy.+b62.*ak2.+b63.*ak3.+b64.*ak4.+b65.*ak5);
    # sixth step
    wk3odes(t+a6*h,ytemp,ak6,mparams,k);
    # accumulate increments w/ proper weights
    yout .= y .+ h.*(c1.*dy.+c3.*ak3.+c4.*ak4.+c6.*ak6);
    # estimate error
    yerr .= h.*(dc1.*dy.+dc3.*ak3.+dc4.*ak4.+dc5.*ak5.+dc6.*ak6);
end

ext: jl


text:
# Displays exit code information
function reader(exit_code)
    if exit_code == 1
        return "CONVERGED: Provided vector b has zero norm. Returned x as vector of zeros."
    elseif exit_code == 2
        return "CONVERGED: Initial guess of vector x results in sufficient residual."
    elseif exit_code == 30
        return "CONVERGED: Sufficient residual achieved."
    elseif exit_code == 31
        return "CONVERGED: Sufficient residual achieved at first step."
    elseif exit_code == 32
        return "CONVERGED: Sufficient residual achieved at second step."
    elseif exit_code == -11
        return "DID NOT CONVERGE: Rho fell below tolerance."
    elseif exit_code == -12
        return "DID NOT CONVERGE: Omega fell below tolerance."
    elseif exit_code == -13
        return "DID NOT CONVERGE: Zero or infinity alpha."
    elseif exit_code == -2
        return "DID NOT CONVERGE: Maximum iterations reached before convergence."
    else
        return "Invalid error code."
    end
end

export reader


text:
function from_vk_call(x::Spec)
    prop = :(x.$(x.name))
    jtype = idiomatic_julia_type(x)
    @match x begin

        # array pointer
        GuardBy(is_arr) => @match jtype begin
            :(Vector{$_}) => :(unsafe_wrap($jtype, $prop, x.$(x.len); own = true))
        end

        GuardBy(is_length) => nothing
        _ => from_vk_call(prop, jtype, x.type)
    end
end

function from_vk_call(prop, jtype, t)
    @match t begin
        :Cstring => :(unsafe_string($prop))
        GuardBy(in(spec_flags.name)) || GuardBy(in(spec_enums.name)) => prop
        GuardBy(in(getproperty.(filter(!isnothing, spec_flags.bitmask), :name))) => :($jtype(UInt32($prop)))
        GuardBy(in(spec_handles.name)) => :($(remove_vk_prefix(t))($prop))
        GuardBy(is_ntuple) && if ntuple_type(t) ∈ filter(x -> x.is_returnedonly, spec_structs).name
        end => :(from_vk.($(remove_vk_prefix(ntuple_type(t))), $prop))
        if follow_constant(t) == jtype
        end => prop
        :(NTuple{$N,$T}) && if is_vulkan_type(T)
        end => broadcast_ex(from_vk_call(prop, idiomatic_julia_type(T), T))
        _ => :(from_vk($jtype, $prop))
    end
end

function vk_call(x::Spec)
    var = wrap_identifier(x.name)
    jtype = idiomatic_julia_type(x)
    @match x begin
        ::SpecStructMember && if x.type == :VkStructureType && parent(x) ∈ keys(structure_types)
        end => structure_types[parent(x)]
        ::SpecStructMember && if is_semantic_ptr(x.type)
        end => :(unsafe_convert($(x.type), $var))
        if is_fn_ptr(x.type)
        end => var
        GuardBy(is_size) && if x.requirement == POINTER_REQUIRED
        end => x.name # parameter converted to a Ref already
        GuardBy(is_opaque) => var
        GuardBy(is_length) => begin
            !x.autovalidity && @debug "Automatic validation was disabled for length parameter $x."
            @match x begin
                GuardBy(is_length_exception) || GuardBy(!is_inferable_length) => var
                _ => :(pointer_length($(wrap_identifier(first(arglen(x)))))) # Julia works with arrays, not pointers, so the length information can directly be retrieved from them
            end
        end
        GuardBy(is_pointer_start) => 0 # always set first* variables to 0, and the user should provide a (sub)array of the desired length
        if x.type ∈ spec_handles.name
        end => var # handled by unsafe_convert in ccall

        # constant pointer to a unique object
        if is_ptr(x.type) &&
           !is_arr(x) &&
           (x.is_constant || (func = func_by_name(x.func); func.type == FTYPE_QUERY && x ≠ last(children(func))))
        end => @match x begin
            if ptr_type(x.type) ∈ [spec_structs.name; spec_unions.name]
            end => var # handled by cconvert and unsafe_convert in ccall
            if x.requirement == OPTIONAL
            end => :($var == $(default(x)) ? $(default(x)) : Ref($var)) # allow optional pointers to be passed as C_NULL instead of a pointer to a 0-valued integer
            _ => :(Ref($var))
        end
        if x.type ∈ [spec_flags.name; spec_enums.name]
        end => var
        if x.type ∈ getproperty.(filter(!isnothing, spec_flags.bitmask), :name)
        end => :($(x.type)($var.val))
        if x.type ∈ extension_types
        end => var
        _ => @match jtype begin
            :String || :Bool || :(Vector{$et}) || if jtype == follow_constant(x.type)
            end => var # conversions are already defined
            if x.type in [spec_structs.name; spec_unions.name] && jtype == struct_name(x.type)
            end => :($var.vks)
            _ => :(to_vk($(x.type), $var)) # fall back to the to_vk function for conversion
        end
    end
end


text:
using Persa

struct RandModel <: Persa.Model
end

Persa.predict(model::RandModel, pref::Preference, user::Int, item::Int) = rand(pref.possibles)

pref = Preference(1:5)

data = [
   UserPreference(1, 1, Rating(4, pref)),
   UserPreference(1, 3, Rating(2, pref)),
   UserPreference(2, 1, Rating(4, pref)),
   UserPreference(2, 2, Rating(5, pref)),
   UserPreference(2, 3, Rating(3, pref)),
   UserPreference(3, 1, Rating(3, pref)),
   UserPreference(3, 2, Rating(4, pref)),
   UserPreference(4, 2, Rating(5, pref)),
   UserPreference(3, 3, Rating(4, pref))
]

dataset = Dataset(data, 4, 3, pref)

model = RandModel()

mach = machine(model, dataset)

mach[2, 3]


text:
using JuMP, Random, LinearAlgebra, JLD, DataFrames, Test, Suppressor, DelimitedFiles, CSV, StatsBase, MosekTools, Compat, Arpack, SCS
include("sdp_bounds.jl")


text:
using MultiScaleTreeGraph
using Documenter
# using CairoMakie
# using PlantGeom

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

makedocs(;
    modules = [MultiScaleTreeGraph],
    authors = "remi.vezy <VEZY@users.noreply.github.com> and contributors",
    repo = "https://github.com/VEZY/MultiScaleTreeGraph.jl/blob/{commit}{path}#{line}",
    sitename = "MultiScaleTreeGraph.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://VEZY.github.io/MultiScaleTreeGraph.jl",
        assets = String[]
    ),
    pages = [
        "Home" => "index.md",
        "Getting started" => "get_started.md",
        "The MTG format" => [
            "Concept" => "the_mtg/mtg_concept.md",
            "File format" => "the_mtg/mtg_format.md",
            "Our implementation" => "the_mtg/our_implementation.md"
        ],
        "Tutorials" => [
            "tutorials/0.read_write.md",
            "tutorials/1.manipulate_node.md",
            "tutorials/2.descendants_ancestors_filters.md",
            "tutorials/3.transform_mtg.md",
            "tutorials/4.convert_mtg.md",
            # "tutorials/5.plotting.md",
            "tutorials/6.add_remove_nodes.md",
        ],
        "API" => "api.md",
    ]
)

deploydocs(;
    repo = "github.com/VEZY/MultiScaleTreeGraph.jl"
)


text:
@testset "single" begin
    @test single([3]) === 3
    @test_throws Exception single([])
    @test_throws Exception single([3, 2])

    @test single((3,)) === 3
    @test_throws Exception single(())
    @test_throws Exception single((3, 2))

    @test single(Dict(1=>3)) === (1=>3)
    @test_throws Exception single(Dict{Int,Int}())
    @test_throws Exception single(Dict(1=>3, 2=>2))
end


text:
module ConicHulls

export Facet, nconic, nface, ftype, gtype
export ConicHull, create_hull, init_hull!, add!, verify

include("Common.jl")
include("Dets.jl")
include("Primitives.jl")
include("Verify.jl")
include("Hulls.jl")
include("Insertion.jl")
include("RefHull.jl")

using .Common
using .Primitives
using .Hulls
using .Insertion

Hulls.hulltype(NC) = Hulls.hulltype(AGen{Int})

Hulls.create_hull(NC::Int) = create_hull(NC, Hulls.hulltype(NC))
Hulls.create_hull(NC::Int, generators) = create_hull(NC, Hulls.hulltype(NC), generators)

end # module


text:
function __init__()
    _print_num_threads_warning()
    return nothing
end

function _print_num_threads_warning()
    sys_nc = VectorizationBase.num_cores()
    jl_nt = Threads.nthreads()
    return _print_num_threads_warning(sys_nc, jl_nt)
end

function _print_num_threads_warning(sys_nc::Integer, jl_nt::Integer)
    if jl_nt < sys_nc
        if !_is_suppress_warning()
            msg = string(
                "The system has $(_pluralize_cores(sys_nc)). ",
                "However, Julia was started with $(_pluralize_threads(jl_nt)). ",
                "We recommend starting Julia with at least $(_pluralize_threads(sys_nc)) ",
                "to take advantage of Gaius's multithreading algorithms. ",
                "To suppress this warning, set the environment variable ",
                "SUPPRESS_GAIUS_WARNING=true",
            )
            @warn(msg)
        end
    end
    return nothing
end

function _string_to_bool(s::AbstractString)
    b = tryparse(Bool, s)
    if b isa Nothing
        return false
    else
        return b::Bool
    end
end

function _is_suppress_warning()
    s = get(ENV, "SUPPRESS_GAIUS_WARNING", "")
    b = _string_to_bool(s)::Bool
    return b
end

function _pluralize(singular::S, plural::S, n::Integer) where {S <: AbstractString}
    if n == 1
        return singular
    else
        return plural
    end
end

function _pluralize_threads(nt::Integer)
    return "$(nt) $(_pluralize("thread", "threads", nt))"
end

function _pluralize_cores(nc::Integer)
    return "$(nc) $(_pluralize("core", "cores", nc))"
end


text:
"""
```
scenario_means_bands(m, scen::AbstractScenario,
    output_vars = [:forecastutobs, :forecastobs, :forecast4qobs,
                   :forecastutpseudo, :forecastpseudo, :forecast4qpseudo];
    verbose = :low, kwargs...)

scenario_means_bands(m, scen::AbstractScenario, output_var; kwargs...)

scenario_means_bands(m, scen::AbstractScenario, output_var, var_name,
    kwargs...)
```

Compute means and bands for model `m` and scenario `scen`. Other methods are for
one `output_var` and one `var_name` respectively.

### Keyword Arguments

- `minimize::Bool`: if `true`, choose shortest interval, otherwise just chop off
  lowest and highest (percent/2)
- `density_bands::Vector{Float64}`: a vector of percent values (between 0 and 1) for
  which to compute density bands
"""
function scenario_means_bands(m::AbstractDSGEModel, scen::AbstractScenario,
                              output_vars::Vector{Symbol} = [:forecastutobs, :forecastutpseudo,
                                                             :forecastobs, :forecastpseudo,
                                                             :forecast4qobs, :forecast4qpseudo];
                              verbose::Symbol = :low,
                              kwargs...)
    # Print
    println(verbose, :low, )
    info_print(verbose, :low, "Computing means and bands for scenario = " * string(scen.key) * "...")
    println(verbose, :low, "Start time: " * string(now()))
    println(verbose, :low, "Means and bands will be saved in " * workpath(m, "scenarios"))

    start_time = time_ns()
    # Revert model alt policy to historical rule
    m <= Setting(:alternative_policy, AltPolicy(:historical, solve, eqcond), false, "apol",
                 "Alternative policy")

    for output_var in output_vars
        print(verbose, :high, "Computing " * string(output_var) * "...")

        # Compute means and bands
        scenario_means_bands(m, scen, output_var; kwargs...)

        output_file = get_scenario_mb_output_file(m, scen, output_var)
        println(verbose, :high, "wrote " * basename(output_file))
    end

    # Print
    total_mb_time     = (time_ns() - start_time)/1e9
    total_mb_time_min = total_mb_time/60
    println(verbose, :low, "\nTotal time to compute scenario means and bands: " * string(total_mb_time_min) * " minutes")
    println(verbose, :low, "Computation of means and bands complete: " * string(now()))
end

function scenario_means_bands(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol; kwargs...)
    # Read metadata
    metadata = get_scenario_mb_metadata(m, scen, output_var)
    date_list = collect(keys(metadata[:date_inds]))
    variable_names = collect(keys(metadata[:indices]))

    # Get to work!
    mapfcn = use_parallel_workers(m) ? pmap : map
    mb_vec = mapfcn(var_name -> scenario_means_bands(m, scen, output_var, var_name; kwargs...),
                    variable_names)

    # Re-assemble pmap outputs
    means = DataFrame(date = date_list)
    bands = Dict{Symbol,DataFrame}()
    for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
        means[!, var_name] = var_means
        bands[var_name] = var_bands
        bands[var_name][!, :date] = date_list
    end
    mb = MeansBands(metadata, means, bands)

    # Write to file
    output_file = get_scenario_mb_output_file(m, scen, output_var)
    output_dir = dirname(output_file)
    isdir(output_dir) || mkpath(output_dir)
    JLD2.jldopen(output_file, "w") do file
        write(file, "mb", mb)
    end

    return mb
end

function scenario_means_bands(m::AbstractDSGEModel, scen::AbstractScenario, output_var::Symbol,
                              var_name::Symbol; minimize::Bool = false,
                              density_bands::Vector{Float64} = [0.5,0.6,0.7,0.8,0.9])
    # Determine class and product
    class = get_class(output_var)
    product = get_product(output_var)
    @assert product in [:forecast, :forecastut, :forecast4q] "Product can only be forecast(ut|4q)"

    # Read in scenario draws
    fcast_series, transform = read_scenario_output(m, scen, class, product, var_name)

    # Reverse transform
    transformed_series = scenario_mb_reverse_transform(fcast_series, transform, product)

    # Compute means and bands
    means = vec(mean(transformed_series, dims=1))
    bands = find_density_bands(transformed_series, density_bands, minimize = minimize)
    return means, bands
end

"""
```
scenario_mb_reverse_transform(fcast_series, transform, product)
```

### Inputs

- `fcast_series::Matrix{Float64}`: `ndraws` x `horizon` matrix of untransformed
  scenario forecasts
- `transform::Function`: reverse transform (possibly no transform at all or 4Q)
- `product::Symbol`: must be one of `:forecast`, `:forecastut`, or `:forecast4q`
"""
function scenario_mb_reverse_transform(fcast_series::Matrix{Float64}, transform::Function,
                                       product::Symbol)
    if product == :forecast4q
        transform4q_scen = get_scenario_transform4q(transform)

        y0s = if transform4q_scen == loggrowthtopct_4q_approx
            # Sum growth rates y_{t-3}, y_{t-2}, y_{t-1}, and y_t
            zeros(3)
        elseif transform4q_scen == logleveltopct_4q_approx
            # Divide log levels y_t by y_{t-4}
            zeros(4)
        else
            Float64[]
        end

        reverse_transform(fcast_series, transform4q_scen;
                          fourquarter = true, y0s = y0s)
    elseif product == :forecast
        transform_scen = get_scenario_transform(transform)
        reverse_transform(fcast_series, transform_scen; y0 = 0.0)

    else
        return fcast_series
    end
end


text:
using FinancialPortfolios, LinearAlgebra, Dictionaries
using Test

w = [0.5, 0.25, 0.25]
nm = ["a", "b", "c"]
plain_dict = Dict(Iterators.zip(nm,w))
plain_dict_sym = Dict(Iterators.zip(Symbol.(nm),w))
plain_dict_int = Dict(Iterators.zip(1:3,w))

r = [0.1, 0.1, -0.1]
const trueret1 = dot(w,r)
const new_weight = w .* (1 .+ r) ./ (1 + trueret1)
const trueret2 = dot(new_weight,r)
const ewret = sum(r) / 3

wgt_dict = Dict(Iterators.zip(nm,new_weight))
wgt_dict_sym = Dict(Iterators.zip(Symbol.(nm),new_weight))
wgt_dict_int = Dict(Iterators.zip(1:3,new_weight))

ret_dict = Dict(Iterators.zip(nm,r))
ret_dict_sym = Dict(Iterators.zip(Symbol.(nm),r))
ret_dict_int = Dict(Iterators.zip(1:3,r))

@testset "Vectors" begin
    FP = FinancialPortfolio(w)
    @test portfolioreturn(FP,r) ≈ trueret1
    update!(FP,r)
    @test portfolioreturn(FP,r) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,r) ≈ ewret
    @test sort(keys(FP)) == [1,2,3]
end

@testset "Dict-String" begin
    ret = ret_dict
    wgts = wgt_dict
    FP = FinancialPortfolio(plain_dict)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == nm
end

@testset "Dict-Symbol" begin
    ret = ret_dict_sym
    wgts = wgt_dict_sym
    FP = FinancialPortfolio(plain_dict_sym)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == Symbol.(nm)
end

@testset "Dict-Int" begin
    ret = ret_dict_int
    wgts = wgt_dict_int
    FP = FinancialPortfolio(plain_dict_int)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == [1,2,3]
end

@testset "HashDictionary-String" begin
    dict = dictionary(plain_dict)
    ret = dictionary(ret_dict)
    wgts = dictionary(wgt_dict)
    FP = FinancialPortfolio(dict)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == nm
end

@testset "HashDictionary-Symbol" begin
    dict = dictionary(plain_dict_sym)
    ret = dictionary(ret_dict_sym)
    wgts = dictionary(wgt_dict_sym)
    FP = FinancialPortfolio(dict)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == Symbol.(nm)
end

@testset "HashDictionary-Int" begin
    dict = dictionary(plain_dict_int)
    ret = dictionary(ret_dict_int)
    wgts = dictionary(wgt_dict_int)
    FP = FinancialPortfolio(dict)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
    @test length(FP) == length(w)
    @test FP == copy(FP)
    FP2 = similar(FP)
    @test portfolioreturn(FP2,ret) ≈ ewret
    @test sort(keys(FP)) == [1,2,3]
end

@testset "HashDictionary-Symbol-NamedTupleReturn" begin
    dict = dictionary(plain_dict_sym)
    ret = (;Iterators.zip(Symbol.(nm),r)...)
    wgts = dictionary(wgt_dict_sym)
    FP = FinancialPortfolio(dict)
    @test portfolioreturn(FP,ret) ≈ trueret1
    update!(FP,ret)
    @test portfolioreturn(FP,ret) ≈ trueret2
end


text:
using NumericalMethodsforEngineers

A = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.]
n = size(A, 1)
x = 0.0

# Convert to upper triang;e
for k in 1:n-1
  if abs(A[k, k]) > 1.0e-6
    for i in k+1:n
      x = A[i, k] / A[k, k]
      A[i, k] = 0.0
      for j in k+1:n
        A[i, j] -= A[k, j] * x
      end
      b[i] -= b[k] * x
    end
  else
    println("Zero pivot found in row $k")
  end
end

println("Modified Matrix")
A |> display
println()

for i in n:-1:1
  x = b[i]
  if i < n
    for j in i+1:n
      x -= A[i, j] * b[j]
    end
  end
  b[i] = x / A[i, i]
end

println("Solution vector:")
b |> display
d = copy(b)

# Restore A and b
A = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.]

c = A\b

println()
@test round(c, 14) == d
@test round(A * c, 14) == b


text:
@doc raw"""
    CircleGroup <: GroupManifold{Circle{ℂ},MultiplicationOperation}

The circle group is the complex circle ([`Circle(ℂ)`](@ref)) equipped with
the group operation of complex multiplication ([`MultiplicationOperation`](@ref)).
"""
const CircleGroup = GroupManifold{ℂ,Circle{ℂ},MultiplicationOperation}

CircleGroup() = GroupManifold(Circle{ℂ}(), MultiplicationOperation())

Base.show(io::IO, ::CircleGroup) = print(io, "CircleGroup()")

invariant_metric_dispatch(::CircleGroup, ::ActionDirection) = Val(true)

default_metric_dispatch(::MetricManifold{ℂ,CircleGroup,EuclideanMetric}) = Val(true)

function compose(G::CircleGroup, p::AbstractVector, q::AbstractVector)
    return map(compose, repeated(G), p, q)
end

compose!(G::CircleGroup, x, p, q) = copyto!(x, compose(G, p, q))

Base.identity(::CircleGroup, p::AbstractVector) = map(one, p)
Base.identity(G::GT, e::Identity{GT}) where {GT<:CircleGroup} = e

identity!(::CircleGroup, q::AbstractVector, p) = copyto!(q, 1)
identity!(::GT, q::AbstractVector, ::Identity{GT}) where {GT<:CircleGroup} = copyto!(q, 1)

Base.inv(G::CircleGroup, p::AbstractVector) = map(inv, repeated(G), p)
Base.inv(G::GT, e::Identity{GT}) where {GT<:CircleGroup} = e

function inverse_translate(
    ::CircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::LeftAction,
)
    return map(/, q, p)
end
function inverse_translate(
    ::CircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::RightAction,
)
    return map(/, q, p)
end

translate_diff(::GT, p, q, X, ::ActionDirection) where {GT<:CircleGroup} = map(*, p, X)
function translate_diff(
    ::GT,
    ::Identity{GT},
    q,
    X,
    ::ActionDirection,
) where {GT<:CircleGroup}
    return X
end

function translate_diff!(G::CircleGroup, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, translate_diff(G, p, q, X, conv))
end

function group_exp(G::CircleGroup, X)
    return map(X) do imθ
        θ = imag(imθ)
        sinθ, cosθ = sincos(θ)
        return Complex(cosθ, sinθ)
    end
end

group_exp!(G::CircleGroup, q, X) = (q .= group_exp(G, X))

function group_log(G::CircleGroup, q)
    return map(q) do z
        cosθ, sinθ = reim(z)
        θ = atan(sinθ, cosθ)
        return θ * im
    end
end

group_log!(G::CircleGroup, X::AbstractVector, q::AbstractVector) = (X .= group_log(G, q))


text:
using IteratorInterfaceExtensions
using TableTraitsUtils
import DataValues

IteratorInterfaceExtensions.isiterable(x::DataFrame) = true
TableTraits.isiterabletable(x::DataFrame) = true

function TableTraits.getiterator(df::DataFrame)
    col_names_raw = [i for i in Pandas.columns(df)]
    col_names = Symbol.(col_names_raw)

    column_data = [eltype(df[i])==String ? [df[i][j] for j=1:length(df)] : values(df[i]) for i in col_names_raw]

    return create_tableiterator(column_data, col_names)
end

TableTraits.supports_get_columns_copy_using_missing(df::DataFrame) = true

function TableTraits.get_columns_copy_using_missing(df::Pandas.DataFrame)
    # return a named tuple of columns here
    col_names_raw = [i for i in Pandas.columns(df)]
    col_names = Symbol.(col_names_raw)
    cols = (Array(eltype(df[i])==String ? [df[i][j] for j=1:length(df)] : df[i]) for i in col_names_raw)
    return NamedTuple{tuple(col_names...)}(tuple(cols...))
end

function _construct_pandas_from_iterabletable(source)
    y = create_columns_from_iterabletable(source, errorhandling=:returnvalue)
    y===nothing && return nothing
    columns, column_names = y[1], y[2]
    cols = Dict{Symbol,Any}(i[1]=>i[2] for i in zip(column_names, columns))

    for (k,v) in pairs(cols)
        if eltype(v)<:DataValues.DataValue
            T = eltype(eltype(v))
            if T<:AbstractFloat
                cols[k] = T[get(i, NaN) for i in v]
            elseif T<:Integer
                cols[k] = Float64[DataValues.isna(i) ? NaN : Float64(get(i)) for i in v]
            else
                throw(ArgumentError("Can't create a Pandas.DataFrame from a source that has missing data."))
            end
        end
    end

    return DataFrame(cols)
end


text:
using Test
using WaveProp.Geometry

@testset "Parametric entity tests" begin
    @testset "curves" begin
        f   = (x) -> SVector(cos(x[1]),sin(x[1]))
        d   = ReferenceLine()
        ent = ParametricEntity(f,d)
        s  =  SVector(rand())
        @test ent(s) == f(s)
        jac = jacobian(ent,s)
        @test jac[1] ≈ -sin(s[1]) && jac[2] ≈ cos(s[1])
        @test ent(s) ≈ normal(ent,s)
    end
    @testset "surfaces" begin
        f      = (x) -> SVector(x[1],x[2],sin(x[1]))
        d      = ReferenceSquare()
        ent    = ParametricEntity(f,d)
        s      =  SVector{2}(rand(2))
        @test ent(s) == f(s)
        @test jacobian(ent,s) ≈ [1 0; 0 1; cos(s[1]) 0]
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Revise
using OptiMode
using LinearAlgebra
using Statistics
using ArrayInterface
using RecursiveArrayTools
using StaticArrays
using HybridArrays
using SparseArrays
using FFTW
using LinearMaps
using GeometryPrimitives
using BenchmarkTools
using ChainRules
using Zygote
using ForwardDiff
using FiniteDifferences
using FiniteDiff
using UnicodePlots
using OhMyREPL
using Crayons.Box       # for color printing
using Zygote: @ignore, dropgrad
using Setfield: @set
using StaticArrays: Dynamic
using IterativeSolvers: bicgstabl
using Rotations: RotY, MRP
using RuntimeGeneratedFunctions
using Tullio
using GLMakie
RuntimeGeneratedFunctions.init(@__MODULE__)

gradRM(fn,in) 			= 	Zygote.gradient(fn,in)[1]
gradFM(fn,in) 			= 	ForwardDiff.gradient(fn,in)
gradFD(fn,in;n=3)		=	FiniteDifferences.grad(central_fdm(n,1),fn,in)[1]
gradFD2(fn,in;rs=1e-2)	=	FiniteDiff.finite_difference_gradient(fn,in;relstep=rs)

derivRM(fn,in) 			= 	Zygote.gradient(fn,in)[1]
derivFM(fn,in) 			= 	ForwardDiff.gradient(fn,in)
derivFD(fn,in;n=3)		=	FiniteDifferences.grad(central_fdm(n,1),fn,in)[1]
derivFD2(fn,in;rs=1e-2)	=	FiniteDiff.finite_difference_derivative(fn,in;relstep=rs)

AD_style = BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style = BOLD*RED_FG
MAN_style = BOLD*GREEN_FG
AD_style_N = NEGATIVE*BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style_N = NEGATIVE*BOLD*RED_FG
MAN_style_N = NEGATIVE*BOLD*GREEN_FG

Δx,Δy,Nx,Ny = 6.0, 4.0, 128, 128;
grid = Grid(Δx,Δy,Nx,Ny)
LNx = rotate(MgO_LiNbO₃,Matrix(MRP(RotY(π/2))),name=:LiNbO₃_X);
LNxN = NumMat(LNx;expr_module=@__MODULE__())
SiO₂N = NumMat(SiO₂;expr_module=@__MODULE__())
Si₃N₄N = NumMat(Si₃N₄;expr_module=@__MODULE__())
AlOxN = NumMat(αAl₂O₃;expr_module=@__MODULE__())
##
rwg(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,Si₃N₄N,SiO₂N,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).
# p = [
#        0.8,                #   top ridge width         `w_top`         [μm]
#        0.4,                #   ridge thickness         `t_core`        [μm]
#        0.75, #0.5,                #   ridge thickness         `t_core`        [μm]
#        π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
#                ];

p = [
       1.2,                #   top ridge width         `w_top`         [μm]
       0.5,                #   ridge thickness         `t_core`        [μm]
       0.85, #0.5,                #   ridge thickness         `t_core`        [μm]
       0.,           #   ridge sidewall angle    `θ`             [radian]
               ];

ω = inv(0.85)
# rwg1 = rwg(p)
ε⁻¹ = smooth(ω,p,:fεs,true,rwg,grid)
kguess =  k_guess(ω,ε⁻¹)
k⃗ = SVector(0.0,0.0,kguess)
M̂ = HelmholtzMap(k⃗, ε⁻¹, grid)
ms = ModeSolver(kguess, ε⁻¹, grid; nev=2, maxiter=300, tol=1e-6)
k1,Hv1 = solve_k(ω,p,rwg,grid;nev=1,eigind=1);
##
geom_fn = rwg
eigind = 1
nev = 2
solve(ω,p,rwg,grid;kguess=kguess,nev=2,eigind=1,maxiter=300,tol=1e-6,f_filter=nothing)
ε,ε⁻¹,nng,nng⁻¹ = deepcopy(smooth(ω,p,(:fεs,:fεs,:fnn̂gs,:fnn̂gs),[false,true,false,true],geom_fn,grid))

smooth(ω,p,:fnn̂gs,true,geom_fn,grid)

##
using KrylovKit
using IterativeSolvers
using DFTK: LOBPCG
evals,evecs,convinfo = eigsolve(x->ms.M̂*x,randn(ComplexF64,size(ms.H⃗,1)),size(ms.H⃗,2),:SR;maxiter=3000,tol=1e-6,krylovdim=40,verbosity=2)
# resIS = lobpcg(ms.eigs_itr,true;not_zeros=false,maxiter=100,tol=1e-6)
resIS = lobpcg(ms.M̂,true,4;not_zeros=false,maxiter=3000,tol=1e-6,P=ms.P̂)
resDF = LOBPCG(ms.M̂,randn(ComplexF64,size(ms.H⃗)),I,ms.P̂,1e-6,300; display_progress=true)

@show evals |> real
@show resIS.λ |> real
@show resDF.λ |> real

##
ω,p,fnames,invert_fn,f_geom,grid,smoothing_fn = ω,p,(:fnĝvds,:fnn̂gs),[false,false],geom_fn,grid,volfrac_smoothing
n_p = length(p)
n_fns=length(fnames)
om_p = vcat(ω,p)
xyz = Zygote.@ignore(x⃗(grid))			# (Nx × Ny × Nz) 3-Array of (x,y,z) vectors at pixel/voxel centers
xyzc = Zygote.@ignore(x⃗c(grid))
vxlmin,vxlmax = vxl_minmax(xyzc)

arr_flatB = Zygote.Buffer(om_p,9,size(grid)...,n_fns)
arr_flat = Zygote.forwarddiff(om_p) do om_p
        geom = f_geom(om_p[2:n_p+1])
        shapes = getfield(geom,:shapes)
        om_inv = inv(first(om_p))
        mat_vals = mapreduce(ss->[ map(f->(mat=SMatrix{3,3}(f(om_inv)); 0.5*(mat+mat')),getfield(geom,ss))... ], hcat, fnames)
        mat_vals_inv = inv.(mat_vals)
        # calcinv = repeat([invert_fn...]',size(mat_vals,1))
        sinds::Matrix{NTuple{4, Int64}} = Zygote.@ignore(proc_sinds(corner_sinds(shapes,xyzc)))
        smoothed_vals_nested = map(sinds,xyz,vxlmin,vxlmax) do sinds,xx,vn,vp
                Tuple(smoothing_fn(sinds,shapes,geom.material_inds,mat_vals,mat_vals_inv,invert_fn,xx,vn,vp))
        end
        smoothed_vals = hcat( [map(x->getindex(x,i),smoothed_vals_nested) for i=1:n_fns]...)
        smoothed_vals_rr = copy(reinterpret(eltype(first(smoothed_vals)),smoothed_vals))
        return smoothed_vals_rr  # new spatially smoothed ε tensor array
end
copyto!(arr_flatB,copy(arr_flat))
arr_flat_r = copy(arr_flatB)
Nx = size(grid,1)
Ny = size(grid,2)
# fn_arrs = [hybridize(view(reshape(arr_flat_r,3,3,size(grid)...,n_fns),1:3,1:3,axes(grid),n),grid) for n=1:n_fns]
views = [view(reshape(arr_flat_r,3,3,size(grid)...,n_fns),1:3,1:3,axes(grid),n) for n in range(n_fns)]
view(reshape(arr_flat_r,3,3,size(grid)...,n_fns),1:3,1:3,axes(grid),1)
arr_flat_r |> size
reshape(arr_flat_r,3,3,size(grid)...,n_fns)
view(reshape(arr_flat_r,3,3,size(grid)...,n_fns),1:3,1:3,axes(grid)...,1)
view(reshape(arr_flat_r,3,3,size(grid)...,n_fns),1:3,1:3,UnitRange.(axes(grid)),1)
##
ωₜ = 1.19
n_guess = 2.1
ω²,H⃗ = solve_ω²(ms,ω*n_guess; nev=1, eigind=1, maxiter=1000, tol=1e-6)
Δω² = ω²[1] - ωₜ^2
∂ω²∂k = 2 * HMₖH(H⃗[:,1],ms.M̂.ε⁻¹,ms.M̂.mag,ms.M̂.mn)
∂ω²∂k = 2 * HMₖH(H⃗[:,1],ms.M̂.ε⁻¹,ms.M̂.mag,view(ms.M̂.mn,:,1,:,:),view(ms.M̂.mn,:,2,:,:))

##
eigind = 1
_solve_Δω²(ms,kguess,ω,eigind=1,maxiter=300,tol=1e-6)
ω²,Hv = solve_ω²(ms,kguess,eigind=1,maxiter=300,tol=1e-6)

f1 = kk -> solve_ω²(ms,kk,eigind=1,maxiter=300,tol=1e-6)[1]
@show omsq = f1(kguess)
@show ∂ω²∂k_fd = derivFD2(f1,kguess;rs=1e-2)

@show ∂ω²∂k = 2 * HMₖH(ms.H⃗[:,eigind],ms.M̂.ε⁻¹,ms.M̂.mag,ms.M̂.mn) / 128^2
m = view(ms.M̂.mn,:,1,:,:)
n = view(ms.M̂.mn,:,2,:,:)
mag = ms.M̂.mag
ε⁻¹ = ms.M̂.ε⁻¹
mn2 = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3])))
permutedims(mn2,(2,1,3,4)) ≈ ms.M̂.mn

HMₖH(H1,ε⁻¹,mag,m,n)
##

E1 = E⃗(ms,1)
H1 = H⃗(ms)[1]
E1s = reinterpret(reshape, SVector{3,Complex{Float64}},  E1)
H1s = reinterpret(reshape, SVector{3,Complex{Float64}},  H1)
S1s = real.( cross.( conj.(E1s), H1s) )
S1 = reshape( reinterpret( Float64,  S1s), (3,size(grid)...))

##

ms.M̂.mn |> size

Hv = copy(ms.H⃗[:,1])
H = reshape(Hv,(2,size(grid)...))

∂ω²∂k / dot(H,H) / 2

dot(H,H) / size(H)[1]
dot(Hv,ms.M̂,Hv) / size(Hv)[1]
dot(H, -kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,ms.M̂.mn,ms.M̂.mag), (2:3) ), real(ε⁻¹)), (2:3)),ms.M̂.mn,ms.M̂.mag) ) / sum(abs2.(H))

kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,ms.M̂.mn,ms.M̂.mag), (2:3) ), real(ε⁻¹)), (2:3)),ms.M̂.mn,ms.M̂.mag)
ifft( ε⁻¹_dot( fft( kx_tc(H,ms.M̂.mn,ms.M̂.mag), (2:3) ), real(ε⁻¹)), (2:3))
size(H1)

dot(H, -kx_ct2( ifft( ε⁻¹_dot( fft( kx_tc2(H,m,n,mag), (2:3) ), real(ε⁻¹)), (2:3)),m,n,mag) ) / sum(abs2.(H))
dot(H, -kx_ct2( ifft( ε⁻¹_dot( fft( kx_tc2(H,n,m,mag), (2:3) ), real(ε⁻¹)), (2:3)),n,m,mag) ) / sum(abs2.(H))

##
dc = copy(ms.M̂.d)
Hc = copy(H)
kx_tc!(dc,Hc,ms.M̂.mn,ms.M̂.mag)
d2 = kx_tc(H,ms.M̂.mn,ms.M̂.mag)
dc ≈ d2
compare_fields(d2,dc)
##
d1 = copy(ms.M̂.d)
e1 = copy(ms.M̂.e)
𝓕! = ms.M̂.𝓕!
𝓕⁻¹! = ms.M̂.𝓕⁻¹!
H1 = copy(H)
kx_tc!(d1,H1,ms.M̂.mn,ms.M̂.mag)
mul!(d1,𝓕!,d1);
eid!(e1,ε⁻¹,d1);
# mul!(e1,𝓕⁻¹!,e1);

e2 =  ε⁻¹_dot( fft( kx_tc(H,ms.M̂.mn,ms.M̂.mag), (2:3) ), real(ε⁻¹))
e1 ≈ e2
compare_fields(e1,e2)

##
d1 = copy(ms.M̂.d)
e1 = copy(ms.M̂.e)
𝓕! = ms.M̂.𝓕!
𝓕⁻¹! = ms.M̂.𝓕⁻¹!
H1 = copy(H)
kx_tc!(d1,H1,ms.M̂.mn,ms.M̂.mag)
mul!(d1,𝓕!,d1);
eid!(e1,ε⁻¹,d1);
mul!(e1,𝓕⁻¹!,e1);

e2 = bfft( ε⁻¹_dot( fft( kx_tc(H,ms.M̂.mn,ms.M̂.mag), (2:3) ), real(ε⁻¹)), (2:3))
e1 ≈ e2
compare_fields(e1,e2)

##
function compare_fields(f1,f2)
        print("field comparison\n")
        @show maximum(abs2.(f1))
        @show maximum(abs2.(f2))
        @show maximum(real(f1))
        @show minimum(real(f1))
        @show maximum(real(f2))
        @show minimum(real(f2))
        @show maximum(imag(f1))
        @show minimum(imag(f1))
        @show maximum(imag(f2))
        @show minimum(imag(f2))
        @show maximum(abs.(real(f1)))
        @show maximum(abs.(real(f2)))
        @show maximum(abs.(imag(f1)))
        @show maximum(abs.(imag(f2)))
        print("\n")
        return
end
##

## Makie example

using CairoMakie
using FileIO

noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
noto_sans_bold = assetpath("fonts", "NotoSans-Bold.ttf")

f = Figure(backgroundcolor = RGBf(0.98, 0.98, 0.98),
    resolution = (1000, 700), font = noto_sans)

ga = f[1, 1] = GridLayout()
gb = f[2, 1] = GridLayout()
gcd = f[1:2, 2] = GridLayout()
gc = gcd[1, 1] = GridLayout()
gd = gcd[2, 1] = GridLayout()

# panel a
axtop = Axis(ga[1, 1])
axmain = Axis(ga[2, 1], xlabel = "before", ylabel = "after")
axright = Axis(ga[2, 2])

labels = ["treatment", "placebo", "control"]
data = randn(3, 100, 2) .+ [1, 3, 5]

for (label, col) in zip(labels, eachslice(data, dims = 1))
    scatter!(axmain, col, label = label)
    density!(axtop, col[:, 1])
    density!(axright, col[:, 2], direction = :y)
end

ylims!(axtop, low = 0)
xlims!(axright, low = 0)

leg = Legend(ga[1, 2], axmain)

hidedecorations!(axtop, grid = false)
hidedecorations!(axright, grid = false)
leg.tellheight = true

colgap!(ga, 10)
rowgap!(ga, 10)

Label(ga[1, 1:2, Top()], "Stimulus ratings", valign = :bottom,
    padding = (0, 0, 5, 0))

# panel b
xs = LinRange(0.5, 6, 50)
ys = LinRange(0.5, 6, 50)
data1 = [sin(x^1.5) * cos(y^0.5) for x in xs, y in ys] .+ 0.1 .* randn.()
data2 = [sin(x^0.8) * cos(y^1.5) for x in xs, y in ys] .+ 0.1 .* randn.()

ax1, hm = contourf(gb[1, 1], xs, ys, data1,
levels = 6)
ax1.title = "Histological analysis"
contour!(xs, ys, data1, levels = 5, color = :black)
hidexdecorations!(ax1)

_, hm2 = contourf(gb[2, 1], xs, ys, data2,
levels = 6)
contour!(xs, ys, data2, levels = 5, color = :black)

cb = Colorbar(gb[1:2, 2], hm, label = "cell group")
low, high = extrema(data1)
edges = range(low, high, length = 7)
centers = (edges[1:6] .+ edges[2:7]) .* 0.5
cb.ticks = (centers, string.(1:6))

cb.alignmode = Mixed(right = 0)
colgap!(gb, 10)
rowgap!(gb, 10)

# panel c
brain = load(assetpath("brain.stl"))

Axis3(gc[1, 1], title = "Brain activation")
m = mesh!(
    brain,
    color = [tri[1][2] for tri in brain for i in 1:3],
    colormap = Reverse(:magma),
)
Colorbar(gc[1, 2], m, label = "BOLD level")

# panel D

axs = [Axis(gd[row, col]) for row in 1:3, col in 1:2]
hidedecorations!.(axs, grid = false, label = false)

for row in 1:3, col in 1:2
    xrange = col == 1 ? (0:0.1:6pi) : (0:0.1:10pi)

    eeg = [sum(sin(pi * rand() + k * x) / k for k in 1:10)
        for x in xrange] .+ 0.1 .* randn.()

    lines!(axs[row, col], eeg, color = (:black, 0.5))
end

axs[3, 1].xlabel = "Day 1"
axs[3, 2].xlabel = "Day 2"
Label(gd[1, :, Top()], "EEG traces", valign = :bottom,
padding = (0, 0, 5, 0))
rowgap!(gd, 10)
colgap!(gd, 10)

for (i, label) in enumerate(["sleep", "awake", "test"])
        Box(gd[i, 3], color = :gray90)
        Label(gd[i, 3], label, rotation = pi/2, tellheight = false)
end
colgap!(gd, 2, 0)
n_day_1 = length(0:0.1:6pi)
n_day_2 = length(0:0.1:10pi)

colsize!(gd, 1, Auto(n_day_1))
colsize!(gd, 2, Auto(n_day_2))
for (label, layout) in zip(["A", "B", "C", "D"], [ga, gb, gc, gd])
Label(layout[1, 1, TopLeft()], label,
        textsize = 26,
        font = noto_sans_bold,
        padding = (0, 5, 5, 0),
        halign = :right)
end
colsize!(f.layout, 1, Auto(0.5))
rowsize!(gcd, 1, Auto(1.5))

f


text:
include("check_clock.jl")
include("update_scores.jl")
include("set_default_values.jl")
include("send_logs.jl")
include("get_rotate_offsets.jl")
include("calc_block_coords.jl")
include("fill_bag.jl")
include("summon_piece.jl")
include("make_buttons.jl")
include("make_scores.jl")
include("make_table.jl")
include("glue_piece.jl")
include("pin_piece.jl")
include("clear_rows.jl")
include("end_round.jl")
include("reset_shadow.jl")
include("can_move.jl")
include("score.jl")
include("step.jl")


text:
"""
    Gu, Gv = arrowscalegrd(xloc, yloc, uscale, vscale)

create u and v grd files for an arrow scale
"""
function arrowscalegrd(xloc, yloc, uscale, vscale)
    X = collect(-2:2) .+ xloc
    Y = collect(-2:2) .+ yloc
    U = zeros(Float64, 5,5)
    V = zeros(Float64, 5,5)
    U[3,3] = uscale
    V[3,3] = vscale

    Gu = GMT.surface([repeat(X,inner=(5,1)) repeat(Y,outer=(5,1)) vec(U)], I=1, R=[X[1],X[end],Y[1],Y[end]])
    Gv = GMT.surface([repeat(X,inner=(5,1)) repeat(Y,outer=(5,1)) vec(V)], I=1, R=[X[1],X[end],Y[1],Y[end]])

    return Gu, Gv
end

"""
    Gu, Gv = arrowgrd(tiles::Vector{VisClaw.AMRGrid}; cutoff=0.1, kwargs...)

create u and v grd files from AMR grid tiles
"""
function arrowgrd(tiles::Vector{VisClaw.AMRGrid}; cutoff=0.1, kwargs...)
    ## region
    region = VisClaw.getlims(tiles)

    ## number of tile
    ntile = length(tiles)

    ## preallocate
    xg = Vector{AbstractVector{Float64}}(undef,ntile)
    yg = Vector{AbstractVector{Float64}}(undef,ntile)
    ug = Vector{AbstractVector{Float64}}(undef,ntile)
    vg = Vector{AbstractVector{Float64}}(undef,ntile)

    ## make vectors of X, Y, U, V for all tiles
    for k = 1:ntile
        X, Y = VisClaw.meshtile(tiles[k])
        xg[k] = vec(X)
        yg[k] = vec(Y)
        ug[k] = vec(tiles[k].u)
        vg[k] = vec(tiles[k].v)
    end
    ## cat all tiles
    xg  = vcat(xg...)
    yg  = vcat(yg...)
    ug  = vcat(ug...)
    vg  = vcat(vg...)

    ## replace NaNs into 0.0
    V = @. sqrt(ug^2 + vg^2)
    ind = @. isnan(V) | (V<cutoff)
    xg[ind] .= 0.0
    yg[ind] .= 0.0
    ug[ind] .= 0.0
    vg[ind] .= 0.0

    ## makegrd
    Gu = GMT.xyz2grd([xg yg ug], I=tiles[1].dx, R=region, kwargs...)
    Gv = GMT.xyz2grd([xg yg vg], I=tiles[1].dy, R=region, kwargs...)

    ## return
    return Gu, Gv
end

arrowgrd(amrs::VisClaw.AMR, istep::Integer; cutoff=0.1, kwargs...) =
arrowgrd(amrs.amr[istep]; cutoff=cutoff, kwargs...)


text:
include("_dcd_off.jl")
include("_load_lp_1.jl")

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using HTTP: HTTP, URI, request
using JSON3: JSON3
using LibPQ: LibPQ, Connection, execute, load!
using DataFrames: DataFrames, DataFrame
using CSV: CSV

conn = Connection("")
# Voting Districts
response = request("GET",
                   URI(scheme = "https",
                       host = "tigerweb.geo.census.gov",
                       path = "/arcgis/rest/services/Census2020/Legislative/MapServer/3/query",
                       query = ["where" => "state=72",
                                "outFields" => "VTD",
                                "f" => "geojson"]))
@assert response.status == 200
json = JSON3.read(response.body)
parse_voting_district(node) = (vtd = node.properties.VTD, geometry = JSON3.write(node.geometry))
voting_districts = parse_voting_district.(json.features);
sort!(getproperty.(voting_districts, :vtd))
filter!(elem -> elem.vtd ≠ "ZZZZZZ", voting_districts);
# execute(conn, "TRUNCATE TABLE precintos_electorales_pr_2020.vtd;")
load!(voting_districts,
      conn,
      string("INSERT INTO precintos_electorales_pr_2020.voting_districts VALUES(",
             "\$1, ST_SetSRID(ST_GeomFromGeoJSON(\$2), 4326)) ",
             "ON CONFLICT DO NOTHING;"))
# Precintos
precintos = CSV.read(joinpath("data", "precintos.tsv"), DataFrame)
# execute(conn, "TRUNCATE TABLE precintos_electorales_pr_2020.precintos;")
load!(precintos,
      conn,
      string("INSERT INTO precintos_electorales_pr_2020.precintos VALUES(",
             "\$1, \$2, \$3, \$4)",
             "ON CONFLICT DO NOTHING;"))


text:
"""RMP制御器"""

# """
# RMPいろいろ
# """
# module RMP

using LinearAlgebra
using ForwardDiff  # 自動微分パッケージ

using Parameters

# export pullbacked_rmp
# export  natural

# export OriginalRMPAttractor
# export OriginalRMPCollisionAvoidance
# export OriginalJointLimitAvoidance
# export RMPfromGDSAttractor
# export RMPfromGDSCollisionAvoidance

"""pullback演算"""
function pullbacked_rmp(f, M, J, dJ, dx)
    pulled_f = J' * (f .- M * dJ * dx)
    pulled_M = J' * M * J
    return pulled_f, pulled_M
end

"""pullback演算"""
function pullbacked_rmp(f, M, J)
    pulled_f = J' * f
    pulled_M = J' * M * J
    return pulled_f, pulled_M
end

"""リーフノード"""
@with_kw mutable struct Leaf{T}
    parent::Int64
    ϕ::Vector{T}
    J::Matrix{T}  # parent -> childへのタスク写像のヤコビアン
    J_dot::Union{Matrix{T}, Nothing}
    x::Vector{T}
    x_dot::Vector{T}
    rmp_param
    f::Vector{T}
    M::Matrix{T}
end

"""ノード"""
@with_kw mutable struct Node{T}
    parent
    children::Vector{Leaf{T}}
    ϕ::Vector{T}
    J::Matrix{T}  # parent -> childへのタスク写像のヤコビアン
    J_dot::Union{Matrix{T}, Nothing}
    x::Vector{T}
    x_dot::Vector{T}
    f::Vector{T}
    M::Matrix{T}
end

"""ルートノード"""
@with_kw mutable struct Root{T}
    children
    q::Vector{T}
    q_dot::Vector{T}
    f::Vector{T}
    M::Matrix{T}
end

"""push
root -> leaf
"""
function pushforward(leaf::Leaf{T}, q::Vector{T}, q_dot::Vector{T}) where T
    leaf.x = leaf.ϕ(q)
    leaf.x_dot = leaf.J(q) * q_dot
end

"""push
node -> leaf
"""
function pushforward(node::Node{T}, q::Vector{T}, q_dot::Vector{T}) where T
    node.x = node.ϕ(q)
    node.x_dot = node.J(q) * q_dot
    node.f = zero(node.f)
    node.M = zero(node.M)
    for child in node.children
        pushforward(child, node.x, node.x_dot)
    end
end

"""push
root -> node
"""
function pushforward(root::Root{T}, q_ddot::Vector{T}, Δt::T) where {T}
    root.q_dot = root.q_dot .+ Δt .* q_ddot
    root.q = root.q .+ .+ Δt .* q_dot
    root.f = zero(root.f)
    root.M = zero(root.M)
    for child in root.children
        pushforward(child, root.q, root.q_dot)
    end
end

function pullback(leaf::Leaf{T}) where {T}
    leaf.f, leaf.M = natural(leaf.rmp_param, leaf.x, leaf.x_dot, leaf.x)

end

### 以下RMP

"""ソフトマックス関数"""
function soft_max(s::T, α::T) where T
    s + 1/α * log(1 + exp(-2 * α * s))
end

"""fromGDSのアトラクターパラメータ"""
@with_kw struct RMPfromGDSAttractor{T}
    max_speed::T
    gain::T
    f_alpha::T
    sigma_alpha::T
    sigma_gamma::T
    wu::T
    wl::T
    alpha::T
    epsilon::T
end

"""目標吸引ポテンシャル2"""
potential_2 = soft_max

"""ポテンシャルの勾配"""
function ∇potential_2(x, η)
    x_norm = norm(x)
    return (1-exp(-2*η*x_norm)) / (1+exp(-2*η*x_norm)) .* x ./ x_norm
end

"""?"""
α_or_γ(x, σ) = exp(-(norm(x))^2 / (2*σ^2))

"""重み行列（fromGDSのアトラクターで使用）"""
w(x, sigma_gamma, wu, wl) = (wu-wl) * α_or_γ(x, sigma_gamma) + wl

"""fromGDSのアトラクター慣性行列"""
function inertia_matrix(x, p::RMPfromGDSAttractor{T}, x₀) where T
    dim = length(x)
    z = x .- x₀
    ∇pot = ∇potential_2(z, p.alpha)
    α = α_or_γ(z, p.sigma_alpha)
    return w(z, p.sigma_gamma, p.wu, p.wl) .* ((1-α) .* ∇pot * ∇pot' .+ (α + p.epsilon).*Matrix{T}(I, dim, dim))
end

"""力用"""
function xMx(x, p::RMPfromGDSAttractor{T}, dx, x₀) where T
    dx' * inertia_matrix(x, p, x₀) * dx
end

"""曲率項"""
function ξ(p::RMPfromGDSAttractor{T}, x, dx, x₀) where T
    # 第一項を計算
    dim = length(x)  # タスク空間の次元
    #A = Matrix{T}(undef, dim, dim)
    A = zeros(T, dim, dim)
    for i in 1:dim
        _M(x) = inertia_matrix(x, p, x₀)[:, i]
        _jacobian_M = ForwardDiff.jacobian(_M, x)
        #println(_jacobian_M)
        A[:, i] = _jacobian_M * dx
    end
    A *= dx

    # 第二項を計算
    _xMx(x) = xMx(x, p, dx, x₀)
    B = 1/2 .* ForwardDiff.gradient(_xMx, x)  # 便利

    return A .- B
end

"""fromGDSのアトラクター力"""
function f(p::RMPfromGDSAttractor{T}, x, dx, x₀, M) where T
    z = x .- x₀
    damp = p.gain / p.max_speed
    #return M * (-p.gain .* soft_normal(z, p.f_alpha) .- damp .* dx) .- ξ(p, x, dx, x₀)
    return M * (-p.gain .* ∇potential_2(z, p.alpha) .- damp .* dx) .- ξ(p, x, dx, x₀)
end

"""fromGDSのアトラクタの自然形式RMPを取得"""
function natural(p::RMPfromGDSAttractor{T}, x, dx, x₀) where T
    M = inertia_matrix(x, p, x₀)
    return f(p, x, dx, x₀, M), M
end

### 障害物回避 ###
"""fromGDSの障害物回避rmpのパラメータ"""
@with_kw struct RMPfromGDSCollisionAvoidance{T}
    rw::T
    sigma::T
    alpha::T
end

#
# w(s) = s^(-4)

# """重み関数の微分"""
# dwds(s) = -4 * s^(-5)

"""距離に関する重み関数"""
function w2(s, rw=1.0)
    if rw - s > 0.0
        return (rw - s)^2 / s
    else
        return zero(s)
    end

end

"""重み関数のs微分"""
function dw2(s, rw)
    if rw - s > 0.0
        return (-2*(rw - s) * s + (rw - s)) / s^2
    else
        return zero(s)
    end
end

"""速度に関する重み関数"""
function u2(ds, σ)
    if ds < 0.0
        return 1 - exp(-ds^2 / (2*σ^2))
    else
        return zero(ds)
    end
end

"""重み関数のds微分"""
function du2(ds, σ)
    if ds < 0.0
        return -exp(-ds^2 / (2*σ^2)) * (-ds / σ^2)
    else
        return zero(ds)
    end
end

function δ(s, ds, σ)
    u2(ds, σ) + 1/2 * ds * du2(ds, σ)
end

"""曲率項"""
function ξ(s, ds, σ, rw)
    1/2 * u2(ds, σ) * dw2(s, rw) * ds^2
end

"""障害物回避ポテンシャル"""
function Φ1(s, α, rw)
    1/2 * α * w2(s, rw)^2
end

"""障害物回避ポテンシャルの勾配"""
function ∇Φ1(s, α, rw)
    α * w2(s, rw) * dw2(s, rw)
end

"""fromGDSの障害物回避力"""
function f(p::RMPfromGDSCollisionAvoidance{T}, s, ds) where T
    return -∇Φ1(s, p.alpha, p.rw) - ξ(s, ds, p.sigma, p.rw)
end

"""fromGDSの障害物回避慣性行列"""
function inertia_matrix(p::RMPfromGDSCollisionAvoidance{T}, s, ds) where T
    return w2(s, p.rw) * δ(s, ds, p.sigma)
end

"""fromGDSの障害物回避rmpの自然形式を取得"""
function natural(p::RMPfromGDSCollisionAvoidance{T}, x, dx, x₀, dx₀=zero(x₀)) where T

    s_vec = x .- x₀
    ds_vec = dx .- dx₀
    s = norm(s_vec)
    ds = (1/s .* dot(s_vec, ds_vec))[1]

    m = inertia_matrix(p, s, ds)
    _f = f(p, s, ds)

    J = -s_vec' ./ s
    J̇ = -s^(-2) .* (ds_vec' .- s_vec' .* ds)

    _f, M = pullbacked_rmp(_f, m, J, J̇, dx)

    return _f, M
end

"""RMPfromGDSのジョイント制限回避のパラメータ"""
@with_kw struct RMPfromGDSJointLimitAvoidance{T}
    gamma_p::T
    gamma_d::T
    lambda::T
    sigma::T
end

function α_upper(dq::T, sigma::T) where T
    1.0 - exp(-max(dq, 0)^2 / (2 * sigma^2))
end

function α_lower(dq::T, sigma::T) where T
    1.0 - exp(-min(dq, 0)^2 / (2 * sigma^2))
end

function _s(q::T, q_u::T, q_l::T) where T
    (q - q_l) / (q_u - q_l)
end

function _dsdq(q::T, q_u::T, q_l::T) where T
    1 / (q_u - q_l)
end

function _d(s::T) where T
    4 * s * (1-s)
end

function _dddq(s::T, dsdq::T) where T
    (4.0 - 8 * s) * dsdq
end

function _b(q::T, dq::T, q_u::T, q_l::T, sigma::T) where T
    s = _s(q, q_u, q_l)
    d = _d(s)
    α_u = α_upper(dq, sigma)
    α_l = α_lower(dq, sigma)
    return s*(α_u * d + (1-α_u)) + (1-s)*(α_l * d + (1-α_l))
end

function _dbdq(q::T, dq::T, q_u::T, q_l::T, sigma::T) where T
    s = _s(q, q_u, q_l)
    d = _d(s)
    α_u = α_upper(dq, sigma)
    α_l = α_lower(dq, sigma)
    dsdq = _dsdq(q, q_u, q_l)
    dddq = _dddq(s, dsdq)
    return (dsdq*(α_u * d + (1-α_u)) + s * dddq) + -dsdq*(α_l * d + (1-α_l)) + (1-s) * dddq
end

"""慣性行列の対角要素"""
function _a(q::T, dq::T, q_u::T, q_l::T, sigma::T) where T
    b = _b(q, dq, q_u, q_l, sigma)
    return b^(-2)
end

function _dadq(q::T, dq::T, q_u::T, q_l::T, sigma::T) where T
    -2 * (_b(q, dq, q_u, q_l, sigma))^(-3) * _dbdq(q, dq, q_u, q_l, sigma)
end

"""fromGDSのジョイント制限回避慣性行列"""
function inertia_matrix(
    p::RMPfromGDSJointLimitAvoidance{T}, q::Vector{T}, dq::Vector{T},
    q_max::Vector{T}, q_min::Vector{T}
    ) where T
    dim = length(q)
    A = zeros(T, dim, dim)

    for i in 1:dim
        A[i, i] = _a(q[i], dq[i], q_max[i], q_min[i], p.sigma)
    end
    #println(A)
    return p.lambda * A
end

"""ジョイント制限回避の曲率項"""
function ξ(q::Vector{T}, dq::Vector{T}, q_max::Vector{T}, q_min::Vector{T}, sigma::T) where T

    #z = Vector{T}(undef, 7)
    dim = length(q)
    z = zero(q)

    for i in 1:dim
        z[i] = 1/2 * _dadq(q[i], dq[i], q_max[i], q_min[i], sigma) * dq[i]^2
    end
    return z
end

"""fromGDSのジョイント制限回避慣性力"""
function f(
    p::RMPfromGDSJointLimitAvoidance{T}, q::Vector{T}, dq::Vector{T},
    q_neutral::Vector{T}, A::Matrix{T},
    ) where T
    ξ_A = ξ(q, dq, q_max, q_min, p.sigma)
    #ξ_A = zeros(T, 7)
    return A * (p.gamma_p * (q_neutral .- q) - p.gamma_d * dq) - ξ_A
end

"""fromGDSのジョイント制限回避RMP"""
function natural(
    p::RMPfromGDSJointLimitAvoidance{T}, q::Vector{T}, dq::Vector{T},
    q_max::Vector{T}, q_min::Vector{T}, q_neutral::Vector{T}
    ) where T
    A = inertia_matrix(p, q, dq, q_max, q_min)
    _f = f(p, q, dq, q_neutral, A)
    return _f, A
end

#end


text:
__precompile__()

module LoopThrottle

export
    @throttle

function loop_throttle_params(; max_rate = 1., min_sleep_time = 0.01)
    min_sleep_time >= 0.001 || error("min_sleep_time must be at least 0.001")
    Float64(max_rate), Float64(min_sleep_time)
end

"""
    @throttle

Throttle a `for` loop or `while` loop by calling `sleep` at the beginning of each
loop iteration (if necessary), so that a designated variable increases
at a rate of at most `max_rate` (compared to wall time).

The first argument that `@throttle` takes is the variable to be rate-limited. This
variable must be in scope inside the loop body. The second argument is a for loop
or while loop. Optionally, the following keyword arguments may be passed in after
the loop argument:
* `max_rate`: specifies the maximum rate at which the designated variable should
increase. Defaults to `1.`.
* `min_sleep_time`: specifies the minimum time to sleep (in seconds). Defaults to
`0.01`, and must be at least 0.001 (due to the accuracy of the `sleep` function).

# Examples
```julia
x = 0
@throttle t for t = 1 : 0.01 : 2
    x += 1
end max_rate = 2.
```
will finish in approximately 0.5 seconds.

```julia
x = 0.
@throttle x for i = 0 : 10
    x += 0.1
end
```
will use the default `max_rate` value of `1.` and thus finish in approximately 1 second.
```julia
i = 0
@throttle i while i <= 10
    println(i)
    i += 1
end min_sleep_time=1.5 max_rate=1
```
will print the numbers from 0 to 10 at an average rate of one per second, while never
sleeping for less than 1.5 second.
"""
macro throttle(t::Symbol, loopexpr::Expr, params::Expr...)
    foreach(params) do expr
        @assert expr.head == :(=)
        @assert expr.args[1] isa Symbol
        expr.head = :kw
        expr.args[2] = esc(expr.args[2])
    end

    setup = quote
        max_rate, min_sleep_time = loop_throttle_params($(params...))
        firstloop = true
        local t0
        local walltime0
    end

    @assert loopexpr.head ∈ (:while, :for)
    loopcondition = loopexpr.args[1]
    loopbody = loopexpr.args[2]
    newloopbody = quote
        if firstloop
            t0 = $(esc(t))
            walltime0 = time()
            firstloop = false
        else
            if !isinf(max_rate)
                Δwalltime = time() - walltime0
                Δt = $(esc(t)) - t0
                sleeptime = Δt / max_rate - Δwalltime
                if sleeptime > min_sleep_time
                    sleep(sleeptime)
                end
            end
        end

        $(esc(loopbody))
    end
    loop = Expr(loopexpr.head, :($(esc(loopcondition))), newloopbody)

    quote
        $setup
        $loop
    end
end

end # module


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0103 Predefined types

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 01.01
using Random
Random.seed!(123)
using InteractiveUtils # loaded automatically when working... interactively

# ## Primitive types

# Primitive types have a fixed number of bits associated to them. Examples of them are `Int64`, `Float64`, `Char`, `UInt64`, `UFloat64`, `Int32`, `Float32`,...
# Even primitive types can be custom defined. See the "custom types" segment !

# ### Char and Strings - `Char`, `String`

a = "Hello World "
b = a[2]
typeof(b)
b = a[3:end]
typeof(b)
#a[2] = 'E' # error !

## info...
ismutable(a) # long story.... https://github.com/JuliaLang/julia/issues/30210
length(a)
findfirst(isequal('o'), a)
findnext(isequal('o'), a, 6)
occursin("world",a)
occursin(lowercase("world"),lowercase(a))
using Unicode
occursin(Unicode.normalize("world", casefold=true),Unicode.normalize(a, casefold=true))
endswith(a,"ld ")
startswith(a,"Mooo")
occursin(r"H*.d ", a)

## modifications..
lowercase(a)
lowercasefirst(a)
split(a, " ")
replace(a,"World" => "Universe")
strip(a)

## concatenation..
a = "Hello"; b= "World"
c = join([a,b]," ")
c = a*" "*b
c = string(a," ",b)
c = "$a $b"              # interpolation

## Conversion..
a = parse(Int64,"2012")
b = string(2019)

# !!! warning
#     Attention not to confuse the `string` function with the `String` type and the `String()` constructor!

## to know more...
methodswith(String,supertypes=true); # where any argument is a String or any parent type (such e.g. AbstractString)
## See also https://docs.julialang.org/en/v1/manual/strings/

# ## Arrays - `Array{T,N}`

# `Array{T,NDims}` A _parameteric_ type where the type of the content and the number of dimensions define the specific type

# !!! tip
#     `Vector{T}` is an alias for `Array{T,1}` and `Matrix{T}` is an alias for `Array{T,2}`, but there isn't anything "special" for 1 or 2 dimensions compared to more dimensions

# ### Vectors - `Array{T,1}`

# One-dimensions arrays in julia are treated as column vector, and, depending on the inner type, they can be stored efficiently contiguously in memory. However they are NOT the same of a single column of a two dimensions array.
# A row vector is necessarily instead a single row of a 2 dimensions array.

a = [1,2,3]
b = [1 ; 2 ; 3 ;;] # This syntax requires Julia >= 1.7
a == b

# #### Initialisation
a = [1,2,3]; #= or =# a = [1;2;3]
a = [1; 6:-2:2; 10] # notes: (a) the semicolon, (b) the range includes BOTH the extremes
a = [[1,2],[3,4]]   # nested vectors. Each elements can have a different lenght, but rules of linear algebra doesn't apply
# !!! danger
#     Don't confuse nested vectors with multi-dimensional arrays!

# Empty (zero-elements) arrays:
a = []
a = Int64[]
a = [1,2,3]
b = [1,"pizza","beer"]
a = Array{Int64,1}()
# !!! warning
#     Watch out for the difference between `a = Array{Int64,1}()` and `a = Array{Int64,1}`
a = Vector{Int64}()
# n-elements initialisation:
n = 3
T = Int64
zeros(n)            # zeros (Float64)
zeros(T,n)          # zeros (casted as type T)
ones(n)             # ones  (Float64)
ones(T,n)           # ones  (casted of type T)
Array{T,1}(undef,n) # garbage
fill(2,3)

# #### Accessing Vectors
a = [10,20,30,40,50,60,70,80,90]
a[1]
a[end]
a[[1; 6:-2:2; 8]] # indicised by a vector of positions

# #### Collecting iterators into vectors
aRange = 3:2:7
a = collect(aRange)
a = [3:2:7;] # alternative
typeof(aRange)
typeof(aRange) <: AbstractArray # Everywhere an AbstractArray is expected, you can provide a range instead

# #### Common operations with vectors
a = [1,2,3]
reverse(a), a[end:-1:1] # Other way to revert an array
vcat([1,2,3],[4,5],6)

push!(a,4)  # add as individual elements
append!(a,5) # add as many new elements

# !!! tip "Functions with exclamation marks"
#     By **convention** functions that modify one of their parameters (and usually this is the first one) are named with an exclamation mark at the end. Remember (most) unicode characters are valid in functions or variable names.

push!([[1,2],[3,4,5]],[6,7])
append!([1,2,3,4,5],[6,7])

pop!(a)
a
popfirst!(a)
deleteat!(a,2)
pushfirst!([2,3],1)
a = [2,1,3,1]
sort(a)   # also `sort!(a)``
unique(a) # also `unique!(a)`
in(1,a)   # also available as operator: `if 2 in a [...] end`
length(a) # number of elements contained in all the dimensions
size(a),size(a,1) # number of elements by dimension
minimum(a)
min(a...)
# !!! tip "..."
#     "..." is called the *splat operator* and it is used to convert the elements in a vector into a tuple of separate elements in a function call, like the example above
min(4,7,3)
minimum([4,7,9])
argmin([4,2,5,2])
sum([1,2,3])
cumsum([1,2,3])
empty!(a) # only for Vectors
using Random
shuffle([1,2,3]) # also shuffle!([1,2,3])
isempty(a)
findall(x -> x == 1, [2,1,3,1]) # anonymous function returning an array of bools, findall then return the indexes
findfirst(x -> x == 1, [2,1,3,1])
myComparitionWith1(i) = i==1
findall(x -> myComparitionWith1(x), [2,1,3,1])
filter(i -> i > 2, [1,2,3,4])

## delete [7,2,5] from an 1:10 array:
data     = [1:10;]
toDelete = [7,5,2]
deleteat!(data, findall(x -> x in toDelete, data))

for (i,value) in enumerate([10,20,30]) # iterator that returns an index/element tuple
    println("$i - $value")
end

names = ["Marc", "Anne"]
sex   = ['M','F']
age   = [25,20]

for zippedElements in zip(names,sex,age) # iterator that returns tuples made with one element per argument of zip
    println(zippedElements)
end

# ### Multidimensional arrays - `Array{T,N}`

# #### Initialisation
a = [[1,2,3] [4,5,6]] # By column, i.e. elements of the first column, elements of the second column, ...
a = [1 4; 2 5; 3 6]   # By row, i.e. elements of the first row, elements of the second row, ...

## Empty (zero-elements) arrays:
a = Array{Float64}(undef, 0, 0, 0) # using explicitly the constructor and explicitly giving zero for each wanted dimension
## n-elements initialisation:
(T,n,m,g,j) = (Int64,1,2,3,'a')
a = zeros(n,m,g)            # n,m,g-elements zeros array
a = ones(n,m,g)             # n,m,g-elements ones array
a = Array{T,3}(undef,n,m,g) # n,m,g-elements array whose content is garbage
a = fill(j,n,m,g)           # n,m,g-elements array of identical j elements
a = rand(n,m,g)             # n,m,g-elements array of of random numbers
a = [3x + 2y + z for x in 1:2, y in 2:3, z in 1:2] # from using list comprehension

# #### Accessing n-dimensional arrays
# Access by indicating position
a = [1 2 3 4; 5 6 7 8; 9 10 11 12]
a[2,1]   # comma to separate the dimensions
# !!! warning
#     Don't confuse `a[i,j]` for selecting an element of a Matrix with `a[i][j]` to select the inner component of a nested array
a[1,2:3]   # with a range on the second dimension
a[[1,3],1] # with a vector of positions in the first dimension
a[2,:]     # with a full range (all values) in the second dimension, i.e. all columns value for row 2
# !!! warning
#     Note that when the data as only one element on a given dimension, julia reduces the dimensions automatically: the result of `a[2,:]` is NOT a row vector (that is a one-row matrix) but a one dimensional array
# Access by a mask (boolean selection)
b = [true false true false; true true true false; true false true false]
a[b] # always flatted array returned (need eventually reshaping, see later)

# #### Funcionality related to dimensions

size(a)              # returns a tuple (i.e. an immutable list) with the sizes of the n dimensions
ndims(a)             # return the number of dimensions of the array (e.g. `2` for a matrix)
reshape(a, 2,3,2)
2*3*2 == length(a)
b = rand(2,1,3)
dropdims(b,dims=(2)) # remove the specified dimensions, provided that the specified dimension have only a single element
permutedims(a)  # "swap" the dimensions
reshape(a,4,3)  # keep the column mayor order
for slice in eachslice(a,dims=1)
    println(slice)
end
a = reshape(1:24, 3,4,2)
for slice in eachslice(a,dims=1)
    println(slice)
end
a = [1 2;3 4; 5 6]
selectdim(a,1,3) # Select an hyperplane on dimension 1 (rows) at position 3. Returns a view

# #### Flat to vector..
a = [1 2; 3 4]
vec(a)        # shadow copy (different view of the underlying data)
reshape(a,4)  # shadow copy
a[:]          # allocate, as all slice operations do

# #### Other functionality related to Arrays
vcat([1 2; 3 4], [5 6; 7 8]) # works also for DataFrames
hcat([1,2,3],[4,5,6])        # works also for DataFrames
a = [1 2; 3 4]
b = similar(a) # garbage inside
cat(a,a,a,dims=3)

# Sort by column (field)
a = [[3,2,1] [20,30,20] [1000,3000,2000] [300,100,200]]
idx = sortperm(a[:,3], rev=true) # return the positions that sort the 3rd column
sortedMatrix = a[idx,:] # selected by using the sorted positions array on the row dimension
sortslices(a, dims=2)   # by cols, using the first row to sort
sortslices(a, dims=1)   # by rows, using the first column to sort
sortslices(a, dims=1, by = x -> (x[2],x[4])) # by rows, using second and fourth columns

# #### Basic linear algebra
using LinearAlgebra
a = [-1,2,3]
b = [4,5,6]
transpose(a)
a'
norm(a)   # l-2 by default
norm(a,1)

# Vector products:
dot(a,b)    # dot, aka "inner" product
a' * b  == dot(a,b)
cross(a,b)  # cross product
a .* b      # element-wise product
a * b'

A = [1 2 3; 6 4 5; 7 8 9]

A^(-1) # inverse
A^2
det(A) # determinant
transpose(A)
A'
# !!! warning
#     Be aware that `transpose` works only for numerical types. When the matrix contains other types (e.g. strings), use `permutedims`

diag(A)
I # operator that automatically scale to the context without actually building the matrix
A*I
B = [1 2; 3 4]; B*I
(evalues, evectors) = eigen(A)

# ## Tuples - `Tuple{T1,T2,...}`

# A "collection" similar to Array but:
# - Immutable
# - Can efficiently host heterogeneous types, as type information is stored for each individual element
# - Linear algebra doesn't apply (use StaticArray.jl package for that)

# Can be tought as anonymous (immutable) structures
# Used to unpack multiple values, e.g. to store on inddividual variables the output of functions with multiple return value

# ### Initialisation

t = (1,2.5,"a",3)
t = 1,2.5,"a",3
typeof(t)

# ### Indexing
t[1]

# ### Conversion
v = [1,2,3]
t = (v...,) # note the comma
v2 = [t...]
v3 = [i[1] for i in t]
v4 = collect(t)
v == v2 == v3 == v4

# ## Named tuples - `NamedTuple{T1,T2,...}`

# As the name suggests, named tuples are collection similar to ordinary tuples, but whose indexing can accept also a name:

nt = (a=1, b=2.5)
#nt = ("a"=1, "b"=2.5)    # Error !
typeof(nt)
nt[1]
nt.a
keys(nt)
values(nt)

for (k,v) in pairs(nt)
    println("$k - $v")
end

# !!! warning
#     The keys of NamedTuples are _symbols_, not strings. We'll see symbols in the metaprogramming segment.

# ### Conversion
k = [:a,:b,:c]
v = [1,2,3]
nt = NamedTuple(Dict(:a=>1,:b=>2,:c=>3))             # Order not guaranteed! We are "lucky" here
nt = NamedTuple(Dict([k=>v for (k,v) in zip(k,v)]))  # Same...
v2 = [nt...]
v3 = [i[1] for i in nt]
v4 = collect(nt)
v == v2 == v3 == v4

# ## Dictionaries - `Dict{Tkey,TValue}`

# Dictionary are mutable, key-referenced containers:
#
# |               | Mutable      | Immutable    |
# | ------------- | ------------ | ------------ |
# | Use position  | Arrays       | Tuples       |
# | Use keys      | Dictionaries | Named tuples |

# !!! warning
#     Note that order is not preserved. For insertion-order preservation see `OrderedDict` and for sorted dictionaries see `SortedDict`, both from the [DataStructures.jl](https://github.com/JuliaCollections/DataStructures.jl) package.

# ### Initialisation

mydict = Dict(); #= or better =#  mydict = Dict{String,Int64}()
mydict = Dict('a'=>1, 'b'=>2, 'c'=>3)

# ### Indexing
mydict['a']
#mydict['d']     # error!
get(mydict,'d',0) # specific a default if key not found

# ### Adding/deleting/checking
mydict['d'] = 4

delete!(mydict,'d')
haskey(mydict, 'a')
in(('a' => 1), mydict)
typeof('a' => 1)

# ### Conversion
# Array - > Dictionary
map((i,j) -> mydict[i]=j, ['e','f','g'], [4,5,6])
mydict
k = [:a,:b,:c]
v = [1,2,3]
mydict = Dict([k=>v for (k,v) in zip(k,v)])
# Dictionary -> Arrays
collect(keys(mydict)) # keys or values alore return an iterator
collect(values(mydict))

# ### Iteration
for (k,v) in mydict
   println("$k is $v")
end

# ## Sets - `Set{T}`
s = Set(); #= or better =# Set{Int64}()
s = Set([1,2,3,4]) # Only a single `2` will be stored
s
push!(s,5)
delete!(s,1)
s2 = Set([4,5,6,7])
intersect(s,s2)
union(s,s2)
setdiff(s,s2)

# ## Date and time - `Date`, `DateTime`

using Dates # a standard library module dealing with dates and times, including periods and calendars

# While a `DateTime` is a more informative object it is also a much more complex one, as it has to deal with problems
# as the time zones and the daylight saving

# ### Creation of a date or time object ("input")
# From current (local) date/time...
todayDate = today()
nowTime = now()
typeof(todayDate)

typeof(nowTime)
Date     <: Dates.AbstractTime
DateTime <: Dates.AbstractTime
nowTimeUnix = time()  # The so-called "Unix time, a 64bit integer counting the number of seconds since the beginning of the year 1970
nowTime = Dates.unix2datetime(nowTimeUnix) # attention this is not local but UTC (Coordinated Universal Time - the Greenwitch time )!
nowTime = Dates.now(Dates.UTC) # an other way to have UTC time

# !!! tip
#     For Time Zone functionalities and conversion, use the external package [TimeZone.jl](https://github.com/JuliaTime/TimeZones.jl/)

# From a String...
christmasDay      = Date("25 Dec 2030", "d u yyyy")
newYearDay        = Date("2031/01/01", "yyyy/m/d")
christmasLunch    = DateTime("2030-12-25T12:30:00", ISODateTimeFormat)   # well known string datetime ISO8601 Format
newYearEvenDinner = DateTime("Sat, 30 Dec 2030 21:30:00", RFC1123Format) # an othe well known format

# Date and time formatters:
# - `y`  Year digit (ef yyyy => 2030, yy => 30)
# - `m`  Month digit (eg m => 3, mm => 03)
# - `u`  Month name (eg "Jan")
# - `U`  Month name long (eg "January")
# - `e`  Day of week (eg "Tue")
# - `E`  Day of week long (eg "Tuesday")
# - `d`  Day of month (eg d => 3, dd => 03)
# - `H`  Hour digit (eg H => 8, HH => 08)
# - `M`  Minute digit (eg M => 0, MM => 00)
# - `S`  Second digit (eg S => 0, SS => 00)
# - `s`  Millisecond digit (eg .000, fixed 3 digits)
#
# Note that the doubling for the digits matters only for using the formatters in the output (see later)

# From a tuple of integers: y, m, d, H, M, S, s ...
d  = Date(2030, 12)  # no need to give it all
dt = DateTime(2030, 12, 31, 9, 30, 0, 0)

# ### Date/Time extraction of information ("output")...

# To String represerntation...
Dates.format(newYearDay, "dd/m/yy")
Dates.format(christmasLunch, "dd/mm/yy H:M:SS")

# Other...
## Date and DateTime...
year(christmasDay)
isleapyear(christmasDay)
month(christmasLunch)
monthname(christmasDay)
day(christmasDay)
dayofweek(christmasDay)
dayname(christmasDay)
daysofweekinmonth(christmasDay) # there are 4 Wednesdays in December 2030
dayofweekofmonth(christmasDay)  # and the 25th is the 4th of them

## Only datetime..
hour(christmasLunch)
minute(christmasLunch)
second(christmasLunch)

# ### Periods and datetime arithmetics
hollidayPeriod = newYearDay - christmasDay  # between dates is in days
longPeriod = Date(2035,6,1) - christmasDay
mealPeriod = DateTime(2030,12,31,23,30) - newYearEvenDinner # between datetime is in milliseconds
#newYearDay - newYearEvenDinner # error! no mixed
convert(DateTime,newYearDay)
convert(Date,newYearEvenDinner) # possible information loss
mealPeriod = convert(DateTime,newYearDay) - newYearEvenDinner
typeof(hollidayPeriod)
typeof(mealPeriod)

# Period hierarchy:
# - `Period`
#   - `DatePeriod`
#     - `Year`
#     - `Month`
#     - `Week`
#     - `Day`
#   - `TimePeriod`
#     - `Hour`
#     - `Minute`
#     - `Second`
#     - `Millisecond`
#     - `Microsecond`
#     - `Nanosecond`

#convert(Dates.Year,longPeriod)      # going up: error or inexacterror
convert(Dates.Millisecond,longPeriod) # going down:  fine
convert(Dates.Millisecond,mealPeriod)

canLongPeriod = Dates.canonicalize(longPeriod)
typeof(canLongPeriod)

# That the best we can get. We can't "easily" decompose a "period" in  years or months... how many days in a month ?
# 31 or 30 ? And in an year ? A `Period` doesn't store information on when it starts.
# However we can make math with periods based on a specific date/time:

nextChristmas                = christmasDay + Year(1) # We can use the constructors of the various periods
christmasPresentsOpeningTime = christmasLunch + Hour(3)
thisWeekdayNextCentury       = dayname(today()+Year(100))

# Ranges
semesters = Dates.Date(2020,1,1):Dates.Month(6):Dates.Date(2022,1,1)
collect(semesters)

# ### Adjustments

# Iterate the past/future days of a date untill some condition is true
sundayBefChristmas = toprev(d -> Dates.dayname(d) == "Sunday", christmasDay)
lastDayOfThisMonth = tonext(d -> Dates.day(d+Day(1)) == 1, today())

# Find first or last weekday of {month,year} of a given date:
lastTuesdayOfThisMonth = tolast(today(), 2, of=Month) # "2" stands for Tuesday
firstSundayOfThisYear  = tofirst(today(), 7, of=Year) # "7" stands for Sunday

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using FrictionModels
using Test

@testset "CoulombModel" begin
    v = -0.5
    mdl = CoulombModel(0.25)
    N = 100.0
    F = friction(mdl, N, v)
    @test F.f == sign(v) * N * mdl.coefficient
end


text:
#specify packages to be used for code
using SimplexGridFactory, ExtendableGrids
using LinearAlgebra
using TetGen
using BenchmarkTools
using Polyester
using GridVisualize, GLMakie
using GeometryBasics
using Colors

#create geometry for domain of interest
function tetrahedralization_of_cube(; xlims=(-1,1), ylims=(-1,1), zlims=(-1,1), maxvolume=.01)
    #set default parameters for a unit cube w/ max element volume size of 0.01
    #these values can be varied when calling the function

    #use TetGen to build geometry
    builder=SimplexGridBuilder(Generator=TetGen)

    #create pts
    p1=point!(builder, xlims[1], ylims[1], zlims[1])
    p2=point!(builder, xlims[2], ylims[1], zlims[1])
    p3=point!(builder, xlims[2], ylims[2], zlims[1])
    p4=point!(builder, xlims[1], ylims[2], zlims[1])
    p5=point!(builder, xlims[1], ylims[1], zlims[2])
    p6=point!(builder, xlims[2], ylims[1], zlims[2])
    p7=point!(builder, xlims[2], ylims[2], zlims[2])
    p8=point!(builder, xlims[1], ylims[2], zlims[2])

    #create faces
    facetregion!(builder,1)
    facet!(builder, p1 ,p2, p3 ,p4)
    facetregion!(builder,2)
    facet!(builder, p5 ,p6, p7 ,p8)
    facetregion!(builder,3)
    facet!(builder, p1, p2, p6 ,p5)
    facetregion!(builder,4)
    facet!(builder, p2, p3 ,p7 ,p6)
    facetregion!(builder,5)
    facet!(builder, p3, p4 ,p8 ,p7)
    facetregion!(builder,6)
    facet!(builder, p4, p1 ,p5 ,p8)

    #mesh over domain
    meshy = simplexgrid(builder,maxvolume=maxvolume)

    #return mesh for whole domain
    return meshy
end

#calculate function values at each of the chebyshev node points
function ufunc(xx,yy,zz)

    #pre-allocate matrix to hold chebyshev node point function evaluations
    uGrid = zeros(length(xx),length(yy),length(zz))

    #for our specific test function, we need to find the midpt of our conecentric spheres
    mid = mean(xx)

    #find function values at each of the chebyshev pts
    #plot concentric spheres a test function for the domain section of interest
    for i = 1:length(xx)
        for j = 1:length(yy)
            for k = 1:length(zz)
                uGridx = xx[i]-mid
                uGridy = yy[j]-mid
                uGridz = zz[k]-mid
                uGrid[i,j,k] = sqrt(uGridx^2+uGridy^2+uGridz^2)
            end
        end
    end

    #return matrix holding function evaluations of all chebyshev pts on a grid
    return uGrid
end

#find wj weight values in x-, y-, z-directions for barycentric lagrange polynomial implementation
#reference paper: https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
function bary(xx,yy,zz)

    #pre-allocate vectors to hold weights in x-, y-, z-directions
    wx = zeros(1,length(xx))
    wy = zeros(1,length(yy))
    wz = zeros(1,length(zz))

    #calculate weights wj in x-direction
    for j = 1:length(xx)
        wjx = 1
        for k = 1:length(xx)
            #exclude the case where we divide by zero
            if !(xx[j] ≈ xx[k])
                wjx = wjx * 1 / (xx[j]-xx[k])
            end
        end
        wx[j] = wjx
    end

    #calculate weights wj in y-direction
    for j = 1:length(yy)
        wjy = 1
        for k = 1:length(yy)
            #exclude the case where we divide by zero
            if !(yy[j] ≈ yy[k])
                wjy = wjy * 1 / (yy[j]-yy[k])
            end
        end
        wy[j] = wjy
    end

    #calculate weights wj in z-direction
    for j = 1:length(zz)
        wjz = 1
        for k = 1:length(zz)
            #exclude the case where we divide by zero
            if !(zz[j] ≈ zz[k])
                wjz = wjz * 1 / (zz[j] - zz[k])
            end
        end
        wz[j] = wjz
    end

    #return vectors holding wj wight values in x-, y-, z-directions
    return wx, wy, wz
end

#calculate function values at the nodal pts found in the original TetGen 3D mesh generation
function unitPlot(xVec,yVec,zVec,xx,yy,zz,uGrid)

    #find the wj weight values in the x-, y-, z-directions
    wx, wy, wz = bary(xx, yy, zz)

    #pre-allocate a vector to hold function values at each of the TetGen nodal pts
    plotty = zeros(length(xVec))

    #find function values for TetGen nodal pts
    @batch for g = 1:length(xVec)

        #specify the x,y,z coordinates of the TetGen nodal pts for this loop
        x=xVec[g]
        y=yVec[g]
        z=zVec[g]

        #set various sum variables to zero for each loop
            #this ensures the summations in equation 4.2 from
            #https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
            #can be calculated
        sum = 0
        sumid = 0
        sumjd = 0
        sumkd = 0

        #sum over wj weights in the z-direction while dividing by z-grid pts
            #this becomes part of the denominator of eqn 4.2
        for k =1:length(zz)
            #avoid dividing by zero
            if !(z ≈ zz[k])
                sumkd = sumkd + wz[k] / (z - zz[k])
            end
        end

        #sum over wj weights in the y-direction while dividing by y-grid pts
            #this becomes part of the denominator of eqn 4.2
        for j =1:length(yy)
            #avoid dividing by zero
            if !(y ≈ yy[j])
                sumjd = sumjd + wy[j] / (y - yy[j])
            end
        end

        #sum over wj weights in the x-direction while dividing by x-grid pts
            #this becomes part of the denominator of eqn 4.2
        for i =1:length(xx)
            #avoid dividing by zero
            if !(x ≈ xx[i])
                sumid = sumid + wx[i] / (x - xx[i])
            end
        end

        #put numerator and denominator of eqn 4.2 together
            #multiply summations together as needed
        for i = 1:length(xx)
            for j = 1:length(yy)
                for k =1:length(zz)
                    #avoid dividing by zero
                    if !(z ≈ zz[k]) && !(y ≈ yy[j]) && !(x ≈ xx[i])
                        sum = sum + uGrid[i,j,k] * (wz[k]/(z - zz[k]))/sumkd * (wy[j]/(y - yy[j]))/sumjd * (wx[i]/(x - xx[i]))/sumid
                    end
                end
            end
        end

        #store the approximation for each TetGen nodal point in the plotty vector
        plotty[g] = sum
    end

    #return a vector containing function values for the TetGen nodal pts
    return plotty
end

#set limits for each of the unit cubes
xyzlims = ([-3*pi,-pi], [-pi,pi], [pi,3*pi])

#specify isosurfaces of interest for each unit block
fmat = [1 1.15 1.25 1.5; 1 1.15 1.25 1.5; 1 1.15 1.25 1.5]
list_of_flevels = [fmat[i,:] for i in 1:size(fmat,1)]

#no planar surfaces are of interest, other than the isosurfaces specified
    #this is important to note, as a non-empty planes vector will alter the
    #output of the marching_tetrahedra algorithm
planes = []

#create empty vectors to hold isosurface mesh data and the colors corresponding to each mesh
list_of_meshes=[]
list_of_colors=[]

#find isosurfaces for each unit cube of interest
for j = 1:3

    #mesh over 3D domain
    meshy = tetrahedralization_of_cube(xlims=xyzlims[j], ylims=xyzlims[j], zlims=xyzlims[j], maxvolume=0.25)
    #extract coordinates and connectivity matrix from mesh
    coord = meshy.components[Coordinates]
    cellnodes = meshy.components[CellNodes] #connectivity matrix
    #find point-wise coordinates on mesh split into x-, y-, z-components
    xVec,yVec,zVec = (meshy.components[Coordinates][i,:] for i = 1:3)

    #create chebyshev nodes
    n = 7
    θ = LinRange(0, pi, n+1)
    xn = sort(cos.(θ))

    #set chebyshev nodes over interval of interest
    xdiff = 0.5 * (xyzlims[j][2] - xyzlims[j][1])
    xx = yy = zz = @. xn * xdiff + (xdiff + xyzlims[j][1])

    #calculate function values at chebyshev node points
    uGrid = ufunc(xx,yy,zz)
    #find function values at node points from original TetGen mesh
        #note that this is done using the chebyshev node grid and the function
        #values at these nodes
    func = unitPlot(xVec,yVec,zVec,xx,yy,zz,uGrid)

    #find each isosurface for a given unit block
    for i=1:length(list_of_flevels)

        #set isosurface value of interest
        level = list_of_flevels[i]

        #extract details about location isosurface pts using marching tetrahedra algorithm
        tet_deets = GridVisualize.marching_tetrahedra(meshy,func,planes,level)
        pts, trngls, fvals = tet_deets

        #output mesh that encompasses isosurface
        makie_triangles = Makie.to_triangles(hcat(getindex.(trngls,1), getindex.(trngls,2), getindex.(trngls,3)))
        makie_pts = Makie.to_vertices(pts)
        iso_mesh = GeometryBasics.normal_mesh(makie_pts, makie_triangles)

        #add newly found mesh to list of meshes
        push!(list_of_meshes, iso_mesh)
        #add a color to correspond w/ mesh of interest
        push!(list_of_colors, RGBA(0.1f0, 0.0f0, 0.0f0, 0.1f0))
    end
end

#plot isosurface meshes for each unit cube
fig, ax, plt = Makie.mesh(list_of_meshes[1],color = list_of_colors[1])
for i = 2:length(list_of_meshes)
    Makie.mesh!(list_of_meshes[i],color = list_of_colors[i])
end

#display final figure
fig

#store data in struct
#create struct
struct outputInfo
    mesh_list
    iso_list
    domain_lims
end

#specify data inputs to struct
isoGraphInfo = outputInfo(list_of_meshes,list_of_flevels,xyzlims)

max_stars_repo_name: panalluri/VisualizationTrixi

max_issues_repo_name: panalluri/VisualizationTrixi

max_forks_repo_name: panalluri/VisualizationTrixi


text:
using Documenter, MIToS

include("literate.jl")

makedocs(
    doctest = true,
    format = :html,
    sitename = "MIToS",
    authors = "Diego Javier Zea",
    modules = [MIToS],
    pages = [
        "Home" => "index.md",
        "Installation.md",
        "Example.md",
        "Modules" => [  "MSA.md",
                        "Information.md",
                        "SIFTS.md",
                        "PDB.md",
                        "Pfam.md"    ],
        "Scripts.md",
        "API" => [  "MSA_API.md",
                    "Information_API.md",
                    "SIFTS_API.md",
                    "PDB_API.md",
                    "Pfam_API.md",
                    "Utils_API.md"          ],
        "Cookbook" => [ "01_Change_B_factors.md",
                        "02_Linking_structural_and_evolutionary_information.md"]
    ]
)

deploydocs(
    repo   = "github.com/diegozea/MIToS.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)