text:
# TODO: docstrings
@doc raw"""

"""
function device_duration_look_ahead_outage!(
    optimization_container::PSI.OptimizationContainer,
    constraint_info::Vector{DeviceDurationConstraintInfo},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
)
    time_steps = PSI.model_time_steps(optimization_container)
    varon = PSI.get_variable(optimization_container, var_names[1])
    varstart = PSI.get_variable(optimization_container, var_names[2])
    varstop = PSI.get_variable(optimization_container, var_names[3])

    name_up = PSI.middle_rename(cons_name, PSI.PSI_NAME_DELIMITER, "up")
    name_down = PSI.middle_rename(cons_name, PSI.PSI_NAME_DELIMITER, "dn")

    set_names = [PSI.get_component_name(ic) for ic in constraint_info]
    con_up = PSI.add_cons_container!(optimization_container, name_up, set_names, time_steps)
    con_down =
        PSI.add_cons_container!(optimization_container, name_down, set_names, time_steps)

    for t in time_steps
        for cont in constraint_info
            name = PSI.get_component_name(cont)
            # Minimum Up-time Constraint
            expr_up =
                JuMP.GenericAffExpr{Float64, PSI._variable_type(optimization_container)}(0)
            for i in (t - cont.duration_data.up + 1):t
                if i in time_steps
                    JuMP.add_to_expression!(expr_up, varon[name, i])
                end
            end
            if t <= cont.duration_data.up
                expr_up += cont.initial_duration[1].value
            end
            con_up[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                varstop[name, t] * cont.duration_data.up - expr_up <=
                cont.timeseries[t] * cont.duration_data.up
            )

            # Minimum Down-time Constraint
            expr_dn =
                JuMP.GenericAffExpr{Float64, PSI._variable_type(optimization_container)}(0)
            for i in (t - cont.duration_data.down + 1):t
                if i in time_steps
                    JuMP.add_to_expression!(expr_dn, (1 - varon[name, i]))
                end
            end
            if t <= cont.duration_data.down
                expr_dn += cont.initial_duration[2].value
            end
            con_down[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                varstart[name, t] * cont.duration_data.down - expr_dn <=
                cont.timeseries[t] * cont.duration_data.down
            )
        end
    end
    return
end

# TODO: docstrings
@doc raw"""

"""
function device_duration_parameters_outage!(
    optimization_container::PSI.OptimizationContainer,
    constraint_info::Vector{DeviceDurationConstraintInfo},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
    param_reference::PSI.UpdateRef,
)
    time_steps = PSI.model_time_steps(optimization_container)
    varon = PSI.get_variable(optimization_container, var_names[1])
    varstart = PSI.get_variable(optimization_container, var_names[2])
    varstop = PSI.get_variable(optimization_container, var_names[3])

    name_up = PSI.middle_rename(cons_name, PSI.PSI_NAME_DELIMITER, "up")
    name_down = PSI.middle_rename(cons_name, PSI.PSI_NAME_DELIMITER, "dn")

    set_names = [PSI.get_component_name(ic) for ic in constraint_info]
    con_up = PSI.add_cons_container!(optimization_container, name_up, set_names, time_steps)
    con_down =
        PSI.add_cons_container!(optimization_container, name_down, set_names, time_steps)

    container_outage = PSI.add_param_container!(
        optimization_container,
        param_reference,
        set_names,
        time_steps,
    )
    param = PSI.get_parameter_array(container_outage)
    multiplier = PSI.get_multiplier_array(container_outage)

    for t in time_steps
        for cont in constraint_info
            name = PSI.get_component_name(cont)

            param[name, t] =
                PJ.add_parameter(optimization_container.JuMPmodel, cont.timeseries[t])
            multiplier[name, t] = cont.multiplier

            # Minimum Up-time Constraint
            expr_up =
                JuMP.GenericAffExpr{Float64, PSI._variable_type(optimization_container)}(0)
            for i in (t - cont.duration_data.up + 1):t
                if i in time_steps
                    JuMP.add_to_expression!(expr_up, varon[name, i])
                end
            end
            if t <= cont.duration_data.up
                expr_up += cont.initial_duration[1].value
            end

            con_up[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                varstop[name, t] * cont.duration_data.up - expr_up <=
                param[name, t] * cont.duration_data.up
            )

            # Minimum Down-time Constraint
            expr_dn =
                JuMP.GenericAffExpr{Float64, PSI._variable_type(optimization_container)}(0)
            for i in (t - cont.duration_data.down + 1):t
                if i in time_steps
                    JuMP.add_to_expression!(expr_dn, (1 - varon[name, i]))
                end
            end
            if t <= cont.duration_data.down
                expr_dn += cont.initial_duration[2].value
            end
            con_down[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                varstart[name, t] * cont.duration_data.down - expr_dn <=
                param[name, t] * cont.duration_data.down
            )
        end
    end
    return
end


text:
using Plasm
using PyCall
p = PyCall.pyimport("pyplasm")

function array2list(cells)
	return PyObject([Any[cell[h] for h=1:length(cell)] for cell in cells])
end

function doublefirst(cells)
	return p["AL"]([cells[1],cells])
end

ext: jl


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DifferentialEquations
using DynamicalSystems
using DelayEmbeddings
using DelimitedFiles
using Random

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we try to fool the reconstruction algorithms by feeding in stochastic
# time series.

## Generate data
Random.seed!(1234)
N = 5000 # time series length
tr = zeros(N,3) # multivariate dataset we will use
alpha, p = .9, 1.
u0 = .2
tr[:,1] = ar_process(u0, alpha, p, N)
tr[:,2] = randn(N)
tr[:,3] = rand(N)

tr = regularize(Dataset(tr))

#method = "mi_min"
method = "ac_zero"

w1 = estimate_delay(tr[:,1], method)
w2 = estimate_delay(tr[:,2], method)
w3 = estimate_delay(tr[:,3], method)

w = maximum(hcat(w1,w2,w3))
taus = 0:100

#Standard TDE
@time Y_tde, τ, _ = optimal_traditional_de(tr[:,1], "afnn", "ac_zero"; τs = taus)
d_tde = size(Y_tde,2)
τ_tde = [(i-1)*τ for i = 1:d_tde]
L_tde = compute_delta_L(tr[:,1], τ_tde, taus[end]; w = w)

#MDOP embedding
println("Computation time MDOP:")
@time begin
    tw1 = mdop_maximum_delay(tr[:,1])
    tw2 = mdop_maximum_delay(tr[:,2])
    tw3 = mdop_maximum_delay(tr[:,3])
    lm1 = DelayEmbeddings.findlocalminima(tw1[2])
    if lm1[1] ≤ 2
        lmm1 = try lm1[2]
        catch
            lm1[1]
        end
    else
        lmm1 = lm1[1]
    end
    lm2 = DelayEmbeddings.findlocalminima(tw2[2])
    if lm2[1] ≤ 2
        lmm2 = try lm2[2]
        catch
            lm2[1]
        end
    else
        lmm2 = lm2[1]
    end
    lm3 = DelayEmbeddings.findlocalminima(tw3[2])
    if lm3[1] ≤ 2
        lmm3 = try lm3[2]
        catch
            lm3[1]
        end
    else
        lmm3 = lm3[1]
    end
    tau_max = maximum(hcat(lmm1, lmm2, lmm3))
    taus_mdop = 0:tau_max

    Y_mdop, τ_vals_mdop, ts_vals_mdop, FNNs_mdop , βs_mdop = mdop_embedding(tr;
                                                        τs = taus_mdop , w = w)
    L_mdop = compute_delta_L(tr, τ_vals_mdop, ts_vals_mdop, taus_mdop[end]; w = w)
end

#Garcia & Almeida
println("Computation time Garcia & Almeida method:")
@time Y_GA, τ_vals_GA, ts_vals_GA, FNNs_GA , ns_GA = garcia_almeida_embedding(tr;
                                                    τs = taus , w = w, T = w)
L_GA = compute_delta_L(tr, τ_vals_GA, ts_vals_GA, taus[end]; w = w)

#Pecuzal
println("Computation time pecuzal method:")
@time Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(tr;
                                                            τs = taus , w = w)
L_pec = sum(Ls_pec)

## compute other evaluation statistics

@time mfnn1, mfnn2, mfnn3, mfnn4, f1, f2, f3, f4, RQA_ref, RQA1, RQA2,
                                            RQA3, RQA4, R_ref, R1, R2, R3, R4 =
        perform_recurrence_analysis(tr, Dataset(Y_tde), Dataset(Y_mdop),
                    Dataset(Y_GA), Dataset(Y_pec); ε = 0.08, w = w, kNN = 10)

writedlm("./scripts/Fooling systems/noise results/Y_GA.csv",Y_GA)
writedlm("./scripts/Fooling systems/noise results/Y_mdop.csv",Y_mdop)
writedlm("./scripts/Fooling systems/noise results/Y_pec.csv",Y_pec)
writedlm("./scripts/Fooling systems/noise results/taus_GA.csv", τ_vals_GA)
writedlm("./scripts/Fooling systems/noise results/taus_mdop.csv", τ_vals_mdop)
writedlm("./scripts/Fooling systems/noise results/taus_pec.csv", τ_vals_pec)
writedlm("./scripts/Fooling systems/noise results/ts_GA.csv", ts_vals_GA)
writedlm("./scripts/Fooling systems/noise results/ts_mdop.csv", ts_vals_mdop)
writedlm("./scripts/Fooling systems/noise results/ts_pec.csv", ts_vals_pec)


text:
using Photometry
using Test
using Random
using Statistics

Random.seed!(8462852)

@testset "Aperture Photometry" begin
    include("aperture/overlap.jl")
    include("aperture/circular.jl")
    include("aperture/photometry.jl")
    include("aperture/elliptical.jl")
    include("aperture/rectangular.jl")
    include("aperture/plots.jl")
end

@testset "Background Estimation" begin
    include("background/background.jl")
    include("background/estimators.jl")
    include("background/interpolators.jl")
end

@testset "Source Detection" begin
    include("detection/detection.jl")
end


text:
using NLsolve

function f!(F, x)
    F[1]=3*x[1]^2+2*x[1]*x[2]+2*x[2]^2+x[3]+3*x[4]-6
    F[2]=2*x[1]^2+x[1]+x[2]^2+3*x[3]+2*x[4]-2
    F[3]=3*x[1]^2+x[1]*x[2]+2*x[2]^2+2*x[3]+3*x[4]-1
    F[4]=x[1]^2+3*x[2]^2+2*x[3]+3*x[4]-3
end

r = mcpsolve(f!, [0., 0., 0., 0.], [Inf, Inf, Inf, Inf],[1.25, 0., 0., 0.5], reformulation = :smooth, autodiff = :forward)

## fef

using SteadyStateDiffEq, DiffEqBase, NLsolve

# `t` is ignored and set to zero
function f(t,u,du)
  du[1] = 2 - 2u[1]
  du[2] = u[1] - 4u[2]
end
u0 = zeros(2)

# Build it directly
prob = SteadyStateProblem(f,u0)
sol = solve(prob,SSRootfind())

# Or just use an existing ODEProblem
prob = ODEProblem(f,u0,(0.0,1.0))
prob = SteadyStateProblem(prob)

# You can pass in an NLsolve function. Needs to be generalized a bit.
sol = solve(prob,SSRootfind(nlsolve = (f,u0) -> NLsolve.nlsolve(f,u0,autodiff=true,method=:newton,iterations=Int(1e6))))


text:
# Functions for calculating time response of a system

# XXX : `step` is a function in Base, with a different meaning than it has
# here. This shouldn't be an issue, but it might be.
@doc """`[y, t, x] = step(sys[, Tf])` or `[y, t, x] = step(sys[, t])`

Calculate the step response of system `sys`. If the final time `Tf` or time
vector `t` is not provided, one is calculated based on the system pole
locations.

`y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`""" ->
function Base.step(sys::StateSpace, t::AbstractVector)
    lt = length(t)
    ny, nu = size(sys)
    nx = sys.nx
    u = ones(size(t))
    x0 = zeros(nx, 1)
    if nu == 1
        y, t, x = lsim(sys, u, t, x0, :zoh)
    else
        x = Array(Float64, lt, nx, nu)
        y = Array(Float64, lt, ny, nu)
        for i=1:nu
            y[:,:,i], t, x[:,:,i] = lsim(sys[:,i], u, t, x0, :zoh)
        end
    end
    return y, t, x
end
function Base.step(sys::TransferFunction{SisoGeneralized}, t::AbstractVector)
    lsim(sys, ones(length(t), sys.nu), t)
end

Base.step(sys::LTISystem, Tf::Real) = step(sys, _default_time_vector(sys, Tf))
Base.step(sys::LTISystem) = step(sys, _default_time_vector(sys))
Base.step(sys::TransferFunction, t::AbstractVector) = step(ss(sys), t::AbstractVector)

@doc """`[y, t, x] = impulse(sys[, Tf])` or `[y, t, x] = impulse(sys[, t])`

Calculate the impulse response of system `sys`. If the final time `Tf` or time
vector `t` is not provided, one is calculated based on the system pole
locations.

`y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`""" ->
function impulse(sys::StateSpace, t::AbstractVector)
    lt = length(t)
    ny, nu = size(sys)
    nx = sys.nx
    u = zeros(size(t))
    if iscontinuous(sys)
        # impulse response equivalent to unforced response of
        # ss(A, 0, C, 0) with x0 = B.
        imp_sys = ss(sys.A, zeros(nx, 1), sys.C, zeros(ny, 1))
        x0s = sys.B
    else
        imp_sys = sys
        x0s = zeros(nx, nu)
        u[1] = 1/sys.Ts
    end
    if nu == 1
        y, t, x = lsim(sys, u, t, x0s, :zoh)
    else
        x = Array(Float64, lt, nx, nu)
        y = Array(Float64, lt, ny, nu)
        for i=1:nu
            y[:,:,i], t, x[:,:,i] = lsim(sys[:,i], u, t, x0s[:,i], :zoh)
        end
    end
    return y, t, x
end
function impulse(sys::TransferFunction{SisoGeneralized}, t::AbstractVector)
    u = zeros(length(t), sys.nu)
    u[1,:] = 1/(t[2]-t[1])
    lsim(sys::TransferFunction{SisoGeneralized}, u, t)
end

impulse(sys::LTISystem, Tf::Real) = impulse(sys, _default_time_vector(sys, Tf))
impulse(sys::LTISystem) = impulse(sys, _default_time_vector(sys))
impulse(sys::TransferFunction, t::AbstractVector) = impulse(ss(sys), t)

@doc """`y, t, x = lsim(sys, u, t[, x0, method])`

`y, t, x, uout = lsim(sys, u::Function, t[, x0, method])`

Calculate the time response of system `sys` to input `u`. If `x0` is ommitted,
a zero vector is used.

`y`, `x`, `uout` has time in the first dimension.

Continuous time systems are discretized before simulation. By default, the
method is chosen based on the smoothness of the input signal. Optionally, the
`method` parameter can be specified as either `:zoh` or `:foh`.

`u` can be a function or a matrix/vector of precalculated control signals.
If `u` is a function, then `u(i,x)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u(t,x) = -L*x` calculated by `lqr`.
To simulate a unit step, use `(i,x)-> 1`, for a ramp, use `(i,x)-> i*h`, for a step at `t=5`, use (i,x)-> (i*h >= 5) etc.

Usage example:
```julia
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = eye(2)
R = eye(1)
L = lqr(sys,Q,R)

u(t,x) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0)
plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")
```
""" ->
function lsim(sys::StateSpace, u::AbstractVecOrMat, t::AbstractVector,
        x0::VecOrMat=zeros(sys.nx, 1), method::Symbol=_issmooth(u) ? :foh : :zoh)
    ny, nu = size(sys)
    nx = sys.nx

    if length(x0) != nx
        error("size(x0) must match the number of states of sys")
    elseif !any(size(u) .== [(length(t), nu) (length(t),)])
        error("u must be of size (length(t), nu)")
    end

    dt = Float64(t[2] - t[1])
    if !iscontinuous(sys) || method == :zoh
        if iscontinuous(sys)
            dsys = c2d(sys, dt, :zoh)[1]
        else
            if sys.Ts != dt
                error("Time vector must match sample time for discrete system")
            end
            dsys = sys
        end
    else
        dsys, x0map = c2d(sys, dt, :foh)
        x0 = x0map*[x0; u[1,:].']
    end
    x = ltitr(dsys.A, dsys.B, map(Float64,u), map(Float64,x0))
    y = (sys.C*(x.') + sys.D*(u.')).'
    return y, t, x
end

function lsim(sys::StateSpace, u::Function, t::AbstractVector,
        x0::VecOrMat=zeros(sys.nx, 1), method::Symbol=:zoh)
    ny, nu = size(sys)
    nx = sys.nx
    if length(x0) != nx
        error("size(x0) must match the number of states of sys")
    elseif size(u(1,x0)) != (nu,) && size(u(1,x0)) != (nu,1)
        error("return value of u must be of size nu")
    end

    dt = Float64(t[2] - t[1])
    if !iscontinuous(sys) || method == :zoh
        if iscontinuous(sys)
            dsys = c2d(sys, dt, :zoh)[1]
        else
            if sys.Ts != dt
                error("Time vector must match sample time for discrete system")
            end
            dsys = sys
        end
    else
        dsys, x0map = c2d(sys, dt, :foh)
        x0 = x0map*[x0; u(1,x0)]
    end
    x,uout = ltitr(dsys.A, dsys.B, u, length(t), map(Float64,x0))
    y = (sys.C*(x.') + sys.D*(uout.')).'
    return y, t, x, uout
end

lsim(sys::TransferFunction, u, t, args...) = lsim(ss(sys), u, t, args...)

function lsim(sys::TransferFunction{SisoGeneralized}, u, t)
    ny, nu = size(sys)
    if !any(size(u) .== [(length(t), nu) (length(t),)])
        error("u must be of size (length(t), nu)")
    end
    ny, nu = size(sys)
    y = Array{Float64}(length(t),ny)
    for i = 1:nu
        for o = 1:ny
            dt = Float64(t[2]-t[1])
            y[:,o] += lsimabstract(sys.matrix[o,i], u[:,i], dt, t[end])
        end
    end
    return y, t
end

@doc """`ltitr(A, B, u[,x0])`

`ltitr(A, B, u::Function, iters[,x0])`

Simulate the discrete time system `x[k + 1] = A x[k] + B u[k]`, returning `x`.
If `x0` is not provided, a zero-vector is used.

If `u` is a function, then `u(i,x)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u=-Lx` calculated by `lqr`. In this case, an integrer `iters` must be provided that indicates the number of iterations.
""" ->
function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::AbstractVecOrMat{T},
        x0::VecOrMat{T}=zeros(T, size(A, 1), 1))
    n = size(u, 1)
    x = Array(T, size(A, 1), n)
    for i=1:n
        x[:,i] = x0
        x0 = A * x0 + B * u[i,:].'
    end
    return x.'
end

function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::Function, iters::Int,
    x0::VecOrMat{T}=zeros(T, size(A, 1), 1))
    x = Array(T, size(A, 1), iters)
    uout = Array(T, size(B, 2), iters)

    for i=1:iters
        x[:,i] = x0
        uout[:,i] = u(i,x0)
        x0 = A * x0 + B * uout[:,i]
    end
    return x.', uout.'
end

# HELPERS:

# TODO: This is a poor heuristic to estimate a "good" time vector to use for
# simulation, in cases when one isn't provided.
function _default_time_vector(sys::LTISystem, Tf::Real=-1)
    Ts = _default_Ts(sys)
    if Tf == -1
        Tf = 100*Ts
    end
    return 0:Ts:Tf
end

function _default_Ts(sys::LTISystem)
    if !iscontinuous(sys)
        Ts = sys.Ts
    elseif !isstable(sys)
        Ts = 0.05
    else
        ps = pole(sys)
        r = minimum([abs(real(ps));0])
        if r == 0.0
            r = 1.0
        end
        Ts = 0.07/r
    end
    return Ts
end

_default_Ts(sys::TransferFunction{SisoGeneralized}) = 0.07

#TODO a reasonable check
_issmooth(u::Function) = false

# Determine if a signal is "smooth"
function _issmooth(u, thresh::AbstractFloat=0.75)
    u = [zeros(1, size(u, 2)); u]       # Start from 0 signal always
    dist = maximum(u) - minimum(u)
    du = abs(diff(u))
    return !isempty(du) && all(maximum(du) <= thresh*dist)
end


text:
using Test
using Bijectors
using Bijectors: Exp, Log, Logit, AbstractNamedBijector, NamedBijector, NamedInverse, NamedCoupling, NamedComposition, Shift

@testset "NamedBijector" begin
    b = NamedBijector((a = Exp(), b = Log()))
    @test b((a = 0.0, b = exp(1.0))) == (a = 1.0, b = 1.0)
end

@testset "NamedComposition" begin
    b = NamedBijector((a = Exp(), ))
    x = (a = 0., b = 1.)

    nc1 = NamedComposition((b, b))
    @test nc1(x) == b(b(x))
    @test logabsdetjac(nc1, x) ≈ logabsdetjac(b, x) + logabsdetjac(b, b(x))

    nc2 = b ∘ b
    @test nc1 == nc2

    inc2 = inv(nc2)
    @test (inc2 ∘ nc2)(x) == x
    @test logabsdetjac((inc2 ∘ nc2), x) ≈ 0.0
end

@testset "NamedCoupling" begin
    nc = NamedCoupling(Val(:b), Val((:a, )), a -> Logit(zero(a), a))
    @inferred NamedCoupling(Val(:b), Val((:a, )), Shift)

    nc = NamedCoupling(:b, (:a, ), a -> Logit(0., a)) # <= not type-inferrable but eh

    @test Bijectors.target(nc) == :b
    @test Bijectors.deps(nc) == (:a, )

    @inferred Bijectors.target(nc)
    @inferred Bijectors.deps(nc)

    x = (a = 1.0, b = 0.5, c = 99999.)
    @test Bijectors.coupling(nc)(x.a) isa Logit
    @test inv(nc)(nc(x)) == x

    @test logabsdetjac(nc, x) == logabsdetjac(Logit(0., 1.), x.b)
    @test logabsdetjac(inv(nc), nc(x)) == -logabsdetjac(nc, x)

    x = (a = 0.0, b = 2.0, c = 1.0)
    nc = NamedCoupling(:c, (:a, :b), (a, b) -> Logit(a, b))
    @test nc(x).c == 0.0
    @test inv(nc)(nc(x)) == x

    x = (a = 0.0, b = 2.0, c = 1.0)
    nc = NamedCoupling(:c, (:b, ), b -> Shift(b))
    @test nc(x).c == 3.0
    @test inv(nc)(nc(x)) == x
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
@show VERSION
@time using DifferentialEquations
@time using Plots
@time using ProgressMeter

# %%
# Descriptions of discretized wave equations with various boundary conditions

function f2_dirichlet!(dv, v, u, p, t)
    f2!(dv, v, u, p, t)
    dirichlet_boundary_condition!(dv, v, u, p, t)
end

function f2_neumann!(dv, v, u, p, t)
    f2!(dv, v, u, p, t)
    neumann_boundary_condition!(dv, v, u, p, t)
end

function f2_absorbing!(dv, v, u, p, t)
    f2!(dv, v, u, p, t)
    absorbing_boundary_condition!(dv, v, u, p, t)
end

function f2!(dv, v, u, p, t)
    h = p[1]; h⁻² = 1/h^2
    @. @views dv[2:end-1,2:end-1] = h⁻²*(
        u[1:end-2,2:end-1] + u[3:end,2:end-1] + u[2:end-1,1:end-2] + u[2:end-1,3:end] - 4u[2:end-1,2:end-1]) # interior
    @. @views dv[1      ,2:end-1] = h⁻²*(u[2      ,2:end-1] + u[1      ,1:end-2] + u[1    ,3:end] - 3u[1      ,2:end-1]) # left edge
    @. @views dv[end    ,2:end-1] = h⁻²*(u[end-1  ,2:end-1] + u[end    ,1:end-2] + u[end  ,3:end] - 3u[end    ,2:end-1]) # right edge
    @. @views dv[2:end-1,1      ] = h⁻²*(u[2:end-1,2      ] + u[1:end-2,1      ] + u[3:end,1    ] - 3u[2:end-1,1      ]) # bottom edge
    @. @views dv[2:end-1,end    ] = h⁻²*(u[2:end-1,end-1  ] + u[1:end-2,end    ] + u[3:end,end  ] - 3u[2:end-1,end    ]) # top edge
    dv[1  ,1  ] = h⁻²*(u[1  ,2    ] + u[2    ,1  ] - 2u[1  ,1  ]) # bottom-left corner
    dv[end,1  ] = h⁻²*(u[end,2    ] + u[end-1,1  ] - 2u[end,1  ]) # bottom-right corner
    dv[1  ,end] = h⁻²*(u[1  ,end-1] + u[2    ,end] - 2u[1  ,end]) # top-left corner
    dv[end,end] = h⁻²*(u[end,end-1] + u[end-1,end] - 2u[end,end]) # top-right corner
end

function dirichlet_boundary_condition!(dv, v, u, p, t)
    h = p[1]; h⁻² = 1/h^2
    @. @views dv[1      ,2:end-1] += h⁻²*(-u[1      ,2:end-1]) # left edge
    @. @views dv[end    ,2:end-1] += h⁻²*(-u[end    ,2:end-1]) # right edge
    @. @views dv[2:end-1,1      ] += h⁻²*(-u[2:end-1,1      ]) # bottom edge
    @. @views dv[2:end-1,end    ] += h⁻²*(-u[2:end-1,end    ]) # top edge
    dv[1  ,1  ] += h⁻²*(-2u[1  ,1  ]) # bottom-left corner
    dv[end,1  ] += h⁻²*(-2u[end,1  ]) # bottom-right corner
    dv[1  ,end] += h⁻²*(-2u[1  ,end]) # top-left corner
    dv[end,end] += h⁻²*(-2u[end,end]) # top-right corner
end

function neumann_boundary_condition!(dv, v, u, p, t)
    h = p[1]; h⁻² = 1/h^2
    @. @views dv[1      ,2:end-1] += h⁻²*(u[2      ,2:end-1] - u[1      ,2:end-1]) # left edge
    @. @views dv[end    ,2:end-1] += h⁻²*(u[end-1  ,2:end-1] - u[end    ,2:end-1]) # right edge
    @. @views dv[2:end-1,1      ] += h⁻²*(u[2:end-1,2      ] - u[2:end-1,1      ]) # bottom edge
    @. @views dv[2:end-1,end    ] += h⁻²*(u[2:end-1,end-1  ] - u[2:end-1,end    ]) # top edge
    dv[1  ,1  ] += h⁻²*(u[1  ,2    ] + u[2    ,1  ] - 2u[1  ,1  ]) # bottom-left corner
    dv[end,1  ] += h⁻²*(u[end,2    ] + u[end-1,1  ] - 2u[end,1  ]) # bottom-right corner
    dv[1  ,end] += h⁻²*(u[1  ,end-1] + u[2    ,end] - 2u[1  ,end]) # top-left corner
    dv[end,end] += h⁻²*(u[end,end-1] + u[end-1,end] - 2u[end,end]) # top-right corner
end

function absorbing_boundary_condition!(dv, v, u, p, t)
    h = p[1]; h⁻¹ = 1/h
    neumann_boundary_condition!(dv, v, u, p, t)
    @. @views dv[1      ,2:end-1] -= h⁻¹*2v[1      ,2:end-1] # left edge
    @. @views dv[end    ,2:end-1] -= h⁻¹*2v[end    ,2:end-1] # right edge
    @. @views dv[2:end-1,1      ] -= h⁻¹*2v[2:end-1,1      ] # bottom edge
    @. @views dv[2:end-1,end    ] -= h⁻¹*2v[2:end-1,end    ] # top edge
    dv[1  ,1  ] -= h⁻¹*4v[1  ,1  ] # bottom-left corner
    dv[end,1  ] -= h⁻¹*4v[end,1  ] # bottom-right corner
    dv[1  ,end] -= h⁻¹*4v[1  ,end] # top-left corner
    dv[end,end] -= h⁻¹*4v[end,end] # top-right corner
end

# %%
function gifanim(x, y, sol;
        umin = minimum(minimum.(u.x[2] for u in sol.u)),
        umax = maximum(maximum.(u.x[2] for u in sol.u)),
        clim = (umin, umax),
        color = :bwr,
        bg=:lightgrey,
        size = (300, 316)
    )
    @show umin, umax
    @show clim
    sleep(0.1)

    ts = range(sol.prob.tspan...; length=length(sol.t)) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
    prog = Progress(length(ts), 0)
    @gif for t in ts
        heatmap(x, y, sol(t).x[2]'; colorbar=false, clim, size, color, bg)
        title!("t = $(round(t, digits=1))"; titlefontsize=12, aspectratio=1)
        plot!(; xlim=extrema(x), ylim=extrema(y))
        next!(prog)
    end
end

function gifanim3d(x, y, sol;
        umin = minimum(minimum.(u.x[2] for u in sol.u)),
        umax = maximum(maximum.(u.x[2] for u in sol.u)),
        clim = (umin, umax),
        color = :bwr,
        bg=:lightgrey,
        camera = (30, 80),
        size = (600, 500)
    )
    @show umin, umax
    @show clim
    sleep(0.1)

    ts = range(sol.prob.tspan...; length=length(sol.t)) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
    prog = Progress(length(ts), 0)
    @gif for t in ts
        surface(x, y, sol(t).x[2]'; colorbar=false, clim, size, color, bg, camera)
        title!("t = $(round(t, digits=1))"; titlefontsize=12)
        plot!(; xlim=extrema(x), ylim=extrema(y), zlim=clim)
        next!(prog)
    end
end

# %%
x = y = 0:0.1:π
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_dirichlet!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%
x = y = 0:0.1:π
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_neumann!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%
x = y = 0:0.1:π
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_absorbing!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%
x = y = range(0, π; length=201)
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_dirichlet!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%
x = y = range(0, π; length=201)
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_neumann!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%
x = y = range(0, π; length=201)
h = step(x)
p = (; h)

f(x, y) = exp(-10(x - 2)^2 - 10(y - 2)^2)
u0 = f.(x', y)
v0 = zero(u0)
tspan = (0.0, 6.0)

prob = SecondOrderODEProblem(f2_absorbing!, v0, u0, tspan, p)
@time sol = solve(prob; saveat=0.05)

gifanim(x, y, sol; clim=(-0.3, 0.3))

# %%
gifanim3d(x, y, sol; clim=(-0.3, 0.3))

# %%

max_stars_repo_path: 0006/Absorbing boundary condition - 2D case.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using DistributedFactorGraphs
using IncrementalInference
using Test
using TimeZones

##

@testset "FileDFG Tests" begin
    for filename in ["/tmp/fileDFG", "/tmp/FileDFGExtension.tar.gz"]
        global dfg
        dfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)

        if typeof(dfg) <: CloudGraphsDFG
            @warn "TEST: Removing all data for user '$(dfg.userId)', robot '$(dfg.robotId)'!"
            clearRobot!!(dfg)
            # Need to recreate this, otherwise there is going to be an issue when creating the nodes.
            createDfgSessionIfNotExist(dfg)
        end

        # Same graph as iifInterfaceTests.jl
        numNodes = 5

        # Arranging some random data
        verts = map(n -> addVariable!(dfg, Symbol("x$n"), ContinuousScalar, tags = [:POSE]), 1:numNodes)
        map(v -> setSolvable!(v, Int(round(rand()))), verts)
        map(v -> getSolverData(verts[4]).solveInProgress=Int(round(rand())), verts)
        map(v -> setSolvedCount!(v, Int(round(10*rand()))), verts)

        # Add some data entries
        map(v -> addDataEntry!(v, BlobStoreEntry(:testing, uuid4(), :store, "", "", "", "", now(localzone()))), verts)
        map(v -> addDataEntry!(v, BlobStoreEntry(:testing2, uuid4(), :store, "", "", "", "", ZonedDateTime(2014, 5, 30, 21, tz"UTC-4"))), verts)

        # Add some PPEs
        ppe1 = MeanMaxPPE(:default, [1.0], [0.2], [3.456])
        ppe2 = MeanMaxPPE(:other, [1.0], [0.2], [3.456], ZonedDateTime(2014, 5, 30, 21, tz"UTC-4"))
        map(v -> addPPE!(dfg, getLabel(v), deepcopy(ppe1)), verts)
        map(v -> addPPE!(dfg, getLabel(v), deepcopy(ppe2)), verts)

        #call update to set it on cloud
        updateVariable!.(dfg, verts)

        facts = map(n -> addFactor!(dfg, [verts[n], verts[n+1]], LinearRelative(Normal(50.0,2.0))), 1:(numNodes-1))
        map(f -> setSolvable!(f, Int(round(rand()))), facts)
        map(f -> f.solverData.eliminated = rand() > 0.5, facts)
        map(f -> f.solverData.potentialused = rand() > 0.5, facts)
        updateFactor!.(dfg, facts)

        #test multihypo
        addFactor!(dfg, [:x1, :x2, :x3], LinearRelative(Normal(50.0,2.0)), multihypo = [1, 0.3, 0.7])

        # Save and load the graph to test.
        saveDFG(filename, dfg)

        @info "Going to load $filename"
        copyDfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)
        @test_throws AssertionError loadDFG!(DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg),"badfilename")

        retDFG = loadDFG!(copyDfg, filename)

        @test issetequal(ls(dfg), ls(retDFG))
        @test issetequal(lsf(dfg), lsf(retDFG))
        for var in ls(dfg)
            @test getVariable(dfg, var) == getVariable(retDFG, var)
        end
        for fact in lsf(dfg)
            @test compareFactor(getFactor(dfg, fact),
                getFactor(retDFG, fact),
                skip=[:timezone, :zone]) # Timezones
                    # :hypotheses, :certainhypo, :multihypo, # Multihypo
                    # :eliminated,
                    # :potentialused])
        end

        # Check data entries
        for v in ls(dfg)
            @test getDataEntries(getVariable(dfg, v)) == getDataEntries(getVariable(retDFG, v))
            @test issetequal(listPPEs(dfg, v), listPPEs(retDFG, v))
            for ppe in listPPEs(dfg, v)
                @test getPPE(dfg, v, ppe) == getPPE(retDFG, v, ppe)
            end
        end

        # test the duplicate order #581
        saveDFG(dfg, filename)
    end
end

##

@testset "FileDFG Regression Tests" begin
    @info "If any of these tests fail, we have breaking changes"
    for file in filter(f -> endswith(f, ".tar.gz"), readdir(joinpath(@__DIR__, "data")))
        loadFile = joinpath(@__DIR__, "data", file)
        global dfg
        dfgCopy = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)
        retDFG = loadDFG!(dfgCopy, loadFile)
        # It should have failed if there were any issues.
        # Trivial check as well
        @test issetequal(ls(retDFG), [:x1, :x2, :x3, :x4, :x5])
        @test issetequal(lsf(retDFG), [:x3x4f1, :x4x5f1, :x1x2f1, :x2x3f1])
    end
end

##


text:
struct EquationSolveError
    eq
    var
    rhs
end

function Base.showerror(io::IO, ese::EquationSolveError)
    print(io, "EquationSolveError: While solving\n\n\t")
    print(io, ese.eq)
    print(io, "\nfor ")
    printstyled(io, var, bold=true)
    print(io, ", obtained RHS\n\n\tt")
    println(io, rhs)
end

function masked_cumsum!(A::Vector)
    acc = zero(eltype(A))
    for i in eachindex(A)
        iszero(A[i]) && continue
        A[i] = (acc += A[i])
    end
end

function contract_variables(graph::BipartiteGraph, var_eq_matching::Matching, eliminated_variables)
    var_rename = ones(Int64, ndsts(graph))
    eq_rename = ones(Int64, nsrcs(graph))
    for v in eliminated_variables
        eq_rename[var_eq_matching[v]] = 0
        var_rename[v] = 0
    end
    masked_cumsum!(var_rename)
    masked_cumsum!(eq_rename)

    dig = DiCMOBiGraph{true}(graph, var_eq_matching)

    # Update bipartite graph
    var_deps = map(1:ndsts(graph)) do v
        [var_rename[v′] for v′ in neighborhood(dig, v, Inf; dir=:in) if var_rename[v′] != 0]
    end

    new_fadjlist = Vector{Int}[
        let new_list = Vector{Int}()
            for v in graph.fadjlist[i]
                if var_rename[v] != 0
                    push!(new_list, var_rename[v])
                else
                    append!(new_list, var_deps[v])
                end
            end
            new_list
        end for i = 1:nsrcs(graph) if eq_rename[i] != 0]

    return BipartiteGraph(new_fadjlist, ndsts(graph) - length(eliminated_variables))
end

"""
    algebraic_variables_scc(sys)

Find strongly connected components of algebraic variables in a system.
"""
function algebraic_variables_scc(sys)
    s = structure(sys)
    if !(s isa SystemStructure)
        sys = initialize_system_structure(sys)
        s = structure(sys)
    end

    # skip over differential equations
    algvars = isalgvar.(Ref(s), 1:ndsts(s.graph))

    var_eq_matching = complete(maximal_matching(s, e->s.algeqs[e], v->algvars[v]))
    var_sccs = find_var_sccs(complete(s.graph), var_eq_matching)

    return var_eq_matching, var_sccs
end


text:
include("util.jl")

lines = getlines("inputs/3.txt")

function line_parser(line)
    (id::String, location::String) = split(line, " @ ")
    id = SubString(id, 2)
    (offset::String, size::String) = split(location, ": ")
    (offset_x::String, offset_y::String) = split(offset, ",")
    (size_x::String, size_y::String) = split(size, "x")
    return (parse(Int64, id),
    parse(Int64, offset_x),
    parse(Int64, offset_y),
    parse(Int64, size_x),
    parse(Int64, size_y))
end

function full_size(line::NTuple{5, Int64})
    width = 1 + line[2] + line[4]
    height = 1 + line[3] + line[5]
    return (width, height)
end

function full_width(line::NTuple{5, Int64})
    return full_size(line)[1]
end

function full_height(line::NTuple{5, Int64})
    return full_size(line)[2]
end

rows = line_parser.(lines)

max_width = maximum(full_width.(rows))
max_height = maximum(full_height.(rows))

cloth = zeros(max_height, max_width)

function add_claim_to_cloth(row::NTuple{5,Int64})
    rows = (1 + row[3]):(row[5]+row[3])
    cols = (1 + row[2]):(row[2]+row[4])
    cloth[rows,cols] .+= 1.;
end

add_claim_to_cloth.(rows)

function num_overlapped_sq_inches()
    overlapped_inches = 0
    for inch in cloth
        if inch > 1.
            overlapped_inches += 1
        end
    end
    return overlapped_inches
end

num_overlapped_sq_inches()

function intact_claim(row::NTuple{5, Int64})::Bool
    rows = (1 + row[3]):(row[5]+row[3])
    cols = (1 + row[2]):(row[2]+row[4])
    subsquare = cloth[rows,cols]
    if maximum(subsquare) <= 1.
        return true
    else
        return false
    end
end

intact_claims = filter(intact_claim, rows)


text:
using Test
using ArnoldiMethod: sylv
using StaticArrays
using LinearAlgebra

@testset "Tiny sylvester equation $T, $s" for T in (Float64, ComplexF64), s = ((2,2),(2,1),(1,2))
    p, q = s
    A = @SMatrix rand(T, p, p)
    B = @SMatrix rand(T, q, q)
    C = @SMatrix rand(T, p, q)

    X, singular = sylv(A, B, C)

    @test A * X - X * B ≈ C
    @test !singular
end

@testset "Singular sylvester $T" for T in (Float64, ComplexF64)
    let
        A = @SMatrix T[1 2; 0 1]
        B = @SMatrix T[1 3; 0 1]
        C = @SMatrix rand(T, 2, 2)
        X, singular = sylv(A, B, C)
        @test singular
    end
    let
        A = @SMatrix T[1]
        B = @SMatrix T[1 3; 0 1]
        C = @SMatrix rand(T, 1, 2)
        X, singular = sylv(A, B, C)
        @test singular
    end
    let
        A = @SMatrix T[1 2; 0 1]
        B = @SMatrix T[1]
        C = @SMatrix rand(T, 2, 1)
        X, singular = sylv(A, B, C)
        @test singular
    end
end


text:
using Wavelets

wf = wavelet(WT.db2, WT.Filter)
wl = wavelet(WT.db2, WT.Lifting)
N = 1024*1024;
x0 = rand(N);
L = maxtransformlevels(N)  # int(log2(N)-2)
tn = 100


text:
using Winston

println("-----------------------------------")
println("  Welcome to Potential Plot.jl")
println("-----------------------------------")

mt_int = ARGS[1]
mt_float = ARGS[2]
xi = ARGS[3]

Data = readcsv("Data/Cosmo_$(mt_int)_$(mt_float)_$(xi).csv")

phi = Data[:,1];
V = Data[:,2];

# Set Figure Axis
function KeyGenerator(phi)
    for (i, elem) in enumerate(phi)
        if abs(elem-5) < 1e-03
            key_index = i
            return key_index
        end
    end
end
key_index = KeyGenerator(phi)
y_key = V[key_index];

# Potential Plot
p = FramedPlot(
    title="Potential Plots",
    xlabel="\\phi",
    ylabel="V",
    xrange=(0,5),
    yrange=(0,y_key+0.1y_key));
C = Curve(phi, V)
setattr(C, "label", "V")
lgnd = Legend(.9, .9, [C])
add(p, C, lgnd)
savefig(p, "Fig/Potential_$(mt_int)_$(mt_float)_$(xi).svg", (1000, 600))
run(`inkscape -z Fig/Potential_$(mt_int)_$(mt_float)_$(xi).svg -e Fig/Potential_$(mt_int)_$(mt_float)_$(xi).png -d 300 --export-background=WHITE`)
run(`rm Fig/Potential_$(mt_int)_$(mt_float)_$(xi).svg`)

max_stars_repo_name: Axect/RGE

max_issues_repo_name: Axect/RGE

max_forks_repo_name: Axect/RGE


text:
module MyModel

    using Mimi
    include("region_parameters.jl")
    include("gross_economy.jl")
    include("emissions.jl")

    export construct_MyModel
    function construct_MyModel()

        m = Model()

        set_dimension!(m, :time, collect(2015:5:2110))
        set_dimension!(m, :regions, [:Region1, :Region2, :Region3])	 # Note that the regions of your model must be specified here

        add_comp!(m, grosseconomy)
        add_comp!(m, emissions)

        update_param!(m, :grosseconomy, :l, l)
        update_param!(m, :grosseconomy, :tfp, tfp)
        update_param!(m, :grosseconomy, :s, s)
        update_param!(m, :grosseconomy, :depk,depk)
        update_param!(m, :grosseconomy, :k0, k0)
        update_param!(m, :grosseconomy, :share, 0.3)

        update_param!(m, :emissions, :sigma, sigma)
        connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

        return m
    end

end #module

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
# -------------------------------------------
# Document: common.jl
# Author: Sasha Petrenko <sap625@mst.edu>
# Description:
#   Types and functions that are used throughout AdaptiveResonance.jl.
# -------------------------------------------

# -------------------------------------------
# Aliases
# -------------------------------------------
#   **Taken from StatsBase.jl**
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}

# Real-numbered aliases
const RealArray{T<:Real, N} = AbstractArray{T, N}
const RealVector{T<:Real} = AbstractArray{T, 1}
const RealMatrix{T<:Real} = AbstractArray{T, 2}

# Integered aliases
const IntegerArray{T<:Integer, N} = AbstractArray{T, N}
const IntegerVector{T<:Integer} = AbstractArray{T, 1}
const IntegerMatrix{T<:Integer} = AbstractArray{T, 2}

# Specifically floating-point aliases
const RealFP = Union{Float32, Float64}

# System's largest native floating point variable
const Float = (Sys.WORD_SIZE == 64 ? Float64 : Float32)

# Acceptable iterators for ART module training and inference
const ARTIterator = Union{UnitRange, ProgressBar}

"""
    DataConfig

Container to standardize training/testing data configuration.

# Parameters
- `setup::Bool`: flag if data has been setup yet or not.
- `mins::RealVector`: list of minimum values for each feature.
- `maxs::RealVector`: list of maximum values for each feature.
- `dim::Int`: dimensionality of the feature vectors (i.e., number of features).
- `dim_comp::Int` complement coded feature dimensionality, twice the size of `dim`.
"""
mutable struct DataConfig
    setup::Bool
    mins::RealVector
    maxs::RealVector
    dim::Int
    dim_comp::Int
end # DataConfig

"""
    DataConfig()

Default constructor for a data configuration, not set up.
"""
function DataConfig()
    DataConfig(
        false,                      # setup
        Array{Float}(undef, 0),     # min
        Array{Float}(undef, 0),     # max
        0,                          # dim
        0                           # dim_comp
    )
end # DataConfig()

"""
    DataConfig(mins::RealVector, maxs::RealVector)

Convenience constructor for DataConfig, requiring only mins and maxs of the features.

This constructor is used when the mins and maxs differ across features. The dimension is inferred by the length of the mins and maxs.
"""
function DataConfig(mins::RealVector, maxs::RealVector)
    # Verify that the mins and maxs are the same length
    length(mins) != length(maxs) && error("Mins and maxs must be the same length.")
    # Get the dimension from one of the arrays
    dim = length(mins)
    # Initialize a Dataconfig with the explicit config
    DataConfig(
        true,   # setup
        mins,   # min
        maxs,   # max
        dim,    # dim
        dim*2   # dim_comp
    )
end # DataConfig(mins::RealVector, maxs::RealVector)

"""
    DataConfig(min::Real, max::Real, dim::Int)

Convenience constructor for DataConfig, requiring only a global min, max, and dim.

This constructor is used in the case that the feature mins and maxs are all the same respectively.
"""
function DataConfig(min::Real, max::Real, dim::Int)
    DataConfig(
        true,               # setup
        repeat([min], dim), # min
        repeat([max], dim), # max
        dim,                # dim
        dim*2               # dim_comp
    )
end # DataConfig(min::Real, max::Real, dim::Int)

"""
    DataConfig(data::RealMatrix)

Convenience constructor for DataConfig, requiring only the data matrix.
"""
function DataConfig(data::RealMatrix)
    # Create an empty dataconfig
    config = DataConfig()

    # Runthe setup upon the config using the data matrix for reference
    data_setup!(config, data)

    # Return the constructed DataConfig
    return config
end # DataConfig(min::Real, max::Real, dim::Int)

"""
    element_min(x::RealVector, W::RealVector)

Returns the element-wise minimum between sample x and weight W.
"""
function element_min(x::RealVector, W::RealVector)
    # Compute the element-wise minimum of two vectors
    return minimum([x W], dims = 2)
end # element_min(x::RealVector, W::RealVector)

"""
    performance(y_hat::IntegerVector, y::IntegerVector)

Returns the categorization performance of y_hat against y.
"""
function performance(y_hat::IntegerVector, y::IntegerVector)
    # Get the number of labels
    n_y = length(y)

    # Check lengths
    if length(y_hat) != n_y
        error("Label vectors must be the same length")
    end

    # Get the number of correct classifications
    n_correct = 0
    for ix = 1:n_y
        if y_hat[ix] == y[ix]
            n_correct += 1
        end
    end

    return n_correct/n_y
end # performance(y_hat::IntegerVector, y::IntegerVector)

"""
    get_data_shape(data::RealArray)

Returns the correct feature dimension and number of samples.
"""
function get_data_shape(data::RealArray)
    # Get the correct dimensionality and number of samples
    if ndims(data) > 1
        dim, n_samples = size(data)
    else
        # dim = 1
        # n_samples = length(data)
        dim = length(data)
        n_samples = 1
    end

    return dim, n_samples
end # get_data_shape(data::RealArray)

"""
    get_n_samples(data::RealArray)

Returns the number of samples, accounting for 1-D and 2-D arrays.
"""
function get_n_samples(data::RealArray)
    # Get the correct dimensionality and number of samples
    if ndims(data) > 1
        n_samples = size(data)[2]
    else
        # n_samples = length(data)
        n_samples = 1
    end

    return n_samples
end # get_n_samples(data::RealArray)

"""
    data_setup!(config::DataConfig, data::RealMatrix)

Sets up the data config for the ART module before training.
"""
function data_setup!(config::DataConfig, data::RealMatrix)
    if config.setup
        @warn "Data configuration already set up, overwriting config"
    else
        config.setup = true
    end

    # Get the correct dimensionality and number of samples
    config.dim, _ = get_data_shape(data)
    config.dim_comp = 2*config.dim

    # Compute the ranges of each feature
    config.mins = [minimum(data[i, :]) for i in 1:config.dim]
    config.maxs = [maximum(data[i, :]) for i in 1:config.dim]
end # data_setup!(config::DataConfig, data::RealMatrix)

"""
    data_setup!(art::ARTModule, data::RealMatrix)

Convenience method for setting up the DataConfig of an ART module in advance.
"""
function data_setup!(art::ARTModule, data::RealMatrix)
    # Modify the DataConfig of the ART module directly
    data_setup!(art.config, data)
end # data_setup!(art::ART, data::RealMatrix)

"""
    get_data_characteristics(data::RealArray ; config::DataConfig=DataConfig())

Get the characteristics of the data, taking account if a data config is passed.

If no DataConfig is passed, then the data characteristics come from the array itself. Otherwise, use the config for the statistics of the data and the data array for the number of samples.
"""
function get_data_characteristics(data::RealArray ; config::DataConfig=DataConfig())
    # If the data is setup, use the config
    if config.setup
        n_samples = get_n_samples(data)
        dim = config.dim
        mins = config.mins
        maxs = config.maxs
    else
        # Get the correct dimensionality and number of samples
        dim, n_samples = get_data_shape(data)
        # Get the ranges for each feature
        mins = [minimum(data[i, :]) for i in 1:dim]
        maxs = [maximum(data[i, :]) for i in 1:dim]
    end
    return dim, n_samples, mins, maxs
end # get_data_characteristics(data::RealArray ; config::DataConfig=DataConfig())

"""
    linear_normalization(data::RealVector ; config::DataConfig=DataConfig())

Normalize the data to the range [0, 1] along each feature.
"""
function linear_normalization(data::RealVector ; config::DataConfig=DataConfig())
    # Vector normalization requires a setup DataConfig
    if !config.setup
        error("Attempting to complement code a vector without a setup DataConfig")
    end

    # Populate a new array with normalized values.
    x_raw = zeros(config.dim)

    # Iterate over each dimension
    for i = 1:config.dim
        denominator = config.maxs[i] - config.mins[i]
        if denominator != 0
            # If the denominator is not zero, normalize
            x_raw[i] = (data[i] .- config.mins[i]) ./ denominator
        else
            # Otherwise, the feature is zeroed because it contains no useful information
            x_raw[i] = zero(Int)
        end
    end
    return x_raw
end # linear_normalization(data::RealArray ; config::DataConfig=DataConfig())

"""
    linear_normalization(data::RealMatrix ; config::DataConfig=DataConfig())

Normalize the data to the range [0, 1] along each feature.
"""
function linear_normalization(data::RealMatrix ; config::DataConfig=DataConfig())
    # Get the data characteristics
    dim, n_samples, mins, maxs = get_data_characteristics(data, config=config)

    # Populate a new array with normalized values.
    x_raw = zeros(dim, n_samples)

    # Verify that all maxs are strictly greater than mins
    if !all(mins .< maxs)
        error("Got a data max index that is smaller than the corresonding min")
    end

    # Iterate over each dimension
    for i = 1:dim
        denominator = maxs[i] - mins[i]
        if denominator != 0
            # If the denominator is not zero, normalize
            x_raw[i, :] = (data[i, :] .- mins[i]) ./ denominator
        else
            # Otherwise, the feature is zeroed because it contains no useful information
            x_raw[i, :] = zeros(length(x_raw[i, :]))
        end
    end
    return x_raw
end # linear_normalization(data::RealMatrix ; config::DataConfig=DataConfig())

"""
    complement_code(data::RealArray ; config::DataConfig=DataConfig())

Normalize the data x to [0, 1] and returns the augmented vector [x, 1 - x].
"""
function complement_code(data::RealArray ; config::DataConfig=DataConfig())
    # Normalize the data
    x_raw = linear_normalization(data, config=config)

    # Complement code the data and return a concatenated matrix
    return vcat(x_raw, 1 .- x_raw)
end # complement_code(data::RealArray ; config::DataConfig=DataConfig())

"""
    get_iterator(opts::ARTOpts, x::Array)
"""
function get_iterator(opts::ARTOpts, x::RealArray)
    # Show a progbar only if the data is 2-D and the option is on
    dim, n_samples = get_data_shape(x)
    single_sample = n_samples == 1

    # Decide if using a progress bar or not
    # Don't use one if either there is a single sample or the display option is off
    prog_bar = single_sample ? false : opts.display

    # Construct the iterator
    iter_raw = 1:n_samples
    iter = prog_bar ? ProgressBar(iter_raw) : iter_raw

    return iter
end # get_iterator(opts::ARTOpts, x::RealArray)

"""
    update_iter(art::ARTModule, iter::ARTIterator, i::Int)
"""
function update_iter(art::ARTModule, iter::ARTIterator, i::Int)
    # Check explicitly for each, as the function definition restricts the types
    if iter isa ProgressBar
        set_description(iter, string(@sprintf("Ep: %i, ID: %i, Cat: %i", art.epoch, i, art.n_categories)))
    elseif iter isa UnitRange
        return
    end
end # update_iter(art::ARTModule, iter::Union{UnitRange, ProgressBar}, i::Int)

"""
    get_sample(x::RealArray, i::Int)

Returns a sample from data array x safely, accounting for 1-D and
"""
function get_sample(x::RealArray, i::Int)
    # Get the shape of the data, irrespective of data type
    dim, n_samples = get_data_shape(x)
    # Get the type shape of the array
    x_dim = ndims(x)
    # Initialize the sample 1-D array with the original dim
    sample = zeros(dim)
    # Short-circuit error if asking for index out of bounds
    i > n_samples && error("Index of data array out of bounds.")
    # Copy the contents of the input if we got a 1-D array
    if x_dim == 1
        sample = x
    # Otherwise, take the correct slice of the 2-D array
    else
        sample = x[:, i]
    end
    return sample
end # get_sample(x::RealArray, i::Int)

"""
    init_train!(x::RealVector, art::ARTModule, preprocessed::Bool)
"""
function init_train!(x::RealVector, art::ARTModule, preprocessed::Bool)
    # If the data is not preprocessed
    if !preprocessed
        # If the data config is not setup, not enough information to preprocess
        if !art.config.setup
            error("$(typeof(art)): cannot preprocess data before being setup.")
        end
        x = complement_code(x, config=art.config)
    # If it is preprocessed and we are not setup
    elseif !art.config.setup
        # Get the dimension of the vector
        dim_comp = length(x)
        # If the complemented dimension is not even, error
        if !iseven(dim_comp)
            error("Declare that the vector is preprocessed, but it is not even")
        end
        # Half the complemented dimension and setup the DataConfig with that
        dim = Int(dim_comp/2)
        art.config = DataConfig(0, 1, dim)
    end
    return x
end # init_train!(x::RealVector, art::ARTModule, preprocessed::Bool)

"""
    init_train!(x::RealMatrix, art::ARTModule, preprocessed::Bool)
"""
function init_train!(x::RealMatrix, art::ARTModule, preprocessed::Bool)
    # If the data is not preprocessed, then complement code it
    if !preprocessed
        # Set up the data config if training for the first time
        !art.config.setup && data_setup!(art.config, x)
        x = complement_code(x, config=art.config)
    end
    return x
end # init_train!(x::RealMatrix, art::ART, preprocessed::Bool)

"""
    init_classify!(x::RealArray, art::ARTModule, preprocessed::Bool)
"""
function init_classify!(x::RealArray, art::ARTModule, preprocessed::Bool)
    # If the data is not preprocessed
    if !preprocessed
        # If the data config is not setup, not enough information to preprocess
        if !art.config.setup
            error("$(typeof(art)): cannot preprocess data before being setup.")
        end
        # Dispatch to the correct complement code method (vector or matrix)
        x = complement_code(x, config=art.config)
    end
    return x
end # init_classify!(x::RealArray, art::ART, preprocessed::Bool)

"""
    classify(art::ARTModule, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)

Predict categories of 'x' using the ART model.

Returns predicted categories 'y_hat.'

# Arguments
- `art::ARTModule`: ART or ARTMAP module to use for batch inference.
- `x::RealMatrix`: the 2-D dataset containing columns of samples with rows of features.
- `preprocessed::Bool=false`: flag, if the data has already been complement coded or not.
- `get_bmu::Bool=false`, flag, if the model should return the best-matching-unit label in the case of total mismatch.

# Examples
```julia-repl
julia> my_DDVFA = DDVFA()
DDVFA
    opts: opts_DDVFA
    ...
julia> x, y = load_data()
julia> train!(my_DDVFA, x)
julia> y_hat = classify(my_DDVFA, y)
```
"""
function classify(art::ARTModule, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)
    # Show a message if display is on
    art.opts.display && @info "Testing $(typeof(art))"

    # Preprocess the data
    x = init_classify!(x, art, preprocessed)

    # Data information and setup
    n_samples = get_n_samples(x)

    # Initialize the output vector
    y_hat = zeros(Int, n_samples)

    # Get the iterator based on the module options and data shape
    iter = get_iterator(art.opts, x)
    for ix = iter
        # Update the iterator if necessary
        update_iter(art, iter, ix)

        # Grab the sample slice
        sample = get_sample(x, ix)

        # Get the classification
        y_hat[ix] = classify(art, sample, preprocessed=true, get_bmu=get_bmu)
    end

    return y_hat
end # classify(art::ARTModule, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)

# -------------------------------------------
# Common Documentation
# -------------------------------------------

@doc raw"""
    classify(art::ARTModule, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)

Predict categories of a single sample of features 'x' using the ART model.

Returns predicted category 'y_hat.'

# Arguments
- `art::ARTModule`: ART or ARTMAP module to use for batch inference.
- `x::RealVector`: the single sample of features to classify.
- `preprocessed::Bool=false`: optional, flag if the data has already been complement coded or not.
- `get_bmu::Bool=false`: optional, flag if the model should return the best-matching-unit label in the case of total mismatch.
"""
classify(art::ARTModule, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)


text:
addprocs(4)

@everywhere using Simplices, SimplexSplitting, InvariantDistribution, TransferEntropy

 # Bogus example to trigger compilation
e_ex = SimplexSplitting.Embedding(InvariantDistribution.invariant_gaussian_embedding(npts = 12))
t_ex = SimplexSplitting.triang_from_embedding(e_ex)
mm = InvariantDistribution.mm_discrete_dense(t_ex)
@assert all(sum(mm, 2) .≈ 1)

te = TransferEntropy.te_from_embedding(e_ex.embedding, 1, n_reps = 200)
te2 = TransferEntropy.te_from_embedding(e_ex.embedding, 1, n_reps = 100)
plot(te.binsizes, median(te.TE, 2), color = "black")
plot!(te2.binsizes, median(te2.TE, 2), color = "red")

# using Plots; plotlyjs()
using ChaoticMaps.Rossler
u0 = rand(6)
rp = RosslerCoupledParams(ω₁ = 1-0.015, ω₂ = 1+0.015, ϵ₂=0.1)
s1 = solve_rossler_coupled(rp = rp, tend = 500, timestep = 2, u0 = u0)
# using ChaoticMaps.Rossler
# p1 = plot(s1, vars=(1, 2, 3), color = "black", width = 1)
# plot!(p1, s1, vars = (4,5,6), color = "blue", width = 1)
#
# p2 = plot(s1, vars=(0, 1), color = "black", width = 1, denseplot = false)
# plot!(p2, s1, vars=(0, 4), color = "blue", width = 1, denseplot = false)
#

sol1 = Array(s1)
X = vec(sol1[1, 100:end])
Y = vec(sol1[4, 100:end])
te_lag = 1
# Triangulate
embedding = hcat(X[1:end-te_lag],
                Y[1:end-te_lag],
                Y[(1 + te_lag):end])

embedding = InvariantDistribution.invariantize_embedding(embedding, max_point_remove = 10)
t = InvariantDistribution.triang_from_embedding(SimplexSplitting.Embedding(embedding))
te_rossler = TransferEntropy.te_from_embedding(embedding, 1)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct AVRFixed <: AVR
        Vf::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a AVR that returns a fixed voltage to the rotor winding

# Arguments
- `Vf::Float64`: Fixed voltage field applied to the rotor winding, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: Fixed AVR has no states
- `n_states::Int64`: Fixed AVR has no states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRFixed <: AVR
    "Fixed voltage field applied to the rotor winding"
    Vf::Float64
    ext::Dict{String, Any}
    "Fixed AVR has no states"
    states::Vector{Symbol}
    "Fixed AVR has no states"
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRFixed(Vf, ext=Dict{String, Any}(), )
    AVRFixed(Vf, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function AVRFixed(; Vf, ext=Dict{String, Any}(), )
    AVRFixed(Vf, ext, )
end

# Constructor for demo purposes; non-functional.
function AVRFixed(::Nothing)
    AVRFixed(;
        Vf=0,
        ext=Dict{String, Any}(),
    )
end

"""Get AVRFixed Vf."""
get_Vf(value::AVRFixed) = value.Vf
"""Get AVRFixed ext."""
get_ext(value::AVRFixed) = value.ext
"""Get AVRFixed states."""
get_states(value::AVRFixed) = value.states
"""Get AVRFixed n_states."""
get_n_states(value::AVRFixed) = value.n_states
"""Get AVRFixed internal."""
get_internal(value::AVRFixed) = value.internal

"""Set AVRFixed Vf."""
set_Vf!(value::AVRFixed, val::Float64) = value.Vf = val
"""Set AVRFixed ext."""
set_ext!(value::AVRFixed, val::Dict{String, Any}) = value.ext = val
"""Set AVRFixed states."""
set_states!(value::AVRFixed, val::Vector{Symbol}) = value.states = val
"""Set AVRFixed n_states."""
set_n_states!(value::AVRFixed, val::Int64) = value.n_states = val
"""Set AVRFixed internal."""
set_internal!(value::AVRFixed, val::InfrastructureSystemsInternal) = value.internal = val


text:
# Original Java version Copyright 2005 Google Inc.
# Julia translation Copyright 2022 Ted Dunning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

struct S1Angle
    radians::Float64

    S1Angle(radians::Number) = new(radians)

    md"""
    Return the angle between two points, which is also equal to the distance between these points
    on the unit sphere. The points do not need to be normalized.
    """
    S1Angle(x::Vector, y::Vector) = new(angle(x, y))
end

degrees(degrees::Number) = S1Angle(deg2rad(degrees))
radians(radians::Number) = S1Angle(radians)

md"""
Returns the angle in radians
"""
radians(x::S1Angle) = x.radians

md"""
Returns the angle in degrees
"""
degrees(x::S1Angle) = rad2deg(x.radians)

md"""
An angle larger than any finite angle
"""
const S1Angle_INFINITY = S1Angle(Inf)

const S1Angle_ZERO = S1Angle(0)

(<)(a::S1Angle, b::S1Angle) = a.radians < b.radians

(≤)(a::S1Angle, b::S1Angle) = a.radians ≤ b.radians

max(a::S1Angle, b::S1Angle) = a > b ? a : b;
min(a::S1Angle, b::S1Angle) = a < b ? a : b;

md"""
Returns the distance along the surface of a sphere of the given radius
"""
distance(a::S1Angle, radius::Number) = a.radians * radius;

(-)(a::S1Angle) = S1Angle(-a.radians)
(+)(a::S1Angle, b::S1Angle) = S1Angle(a.radians + b.radians)
(-)(a::S1Angle, b::S1Angle) = S1Angle(a.radians - b.radians)
(*)(k::Number, a::S1Angle) = S1Angle(k * a.radians)
(/)(a::S1Angle, k::Number) = S1Angle(a.radians / k)
Base.cos(a::S1Angle) = cos(a.radians)
Base.sin(a::S1Angle) = sin(a.radians)
Base.tan(a::S1Angle) = tan(a.radians)

md"""
Returns the angle normalized to the range (-180, 180] degrees.
"""
normalizeRange(a::S1Angle) = S1Angle(normalizeRange(a.radians))

function normalizeRange(r::Number)
    x = rem(r, 2π, RoundNearest)
    if x ≤ -S2_M_PI
        return S2_M_PI
    else
        return x
    end
end


text:
@testset "interface" begin
    problem = MultivariateProblems.UnconstrainedProblems.examples["Exponential"]
    f = MVP.objective(problem)
    g! = MVP.gradient(problem)
    h! = MVP.hessian(problem)
    nd = NonDifferentiable(f, fill!(similar(problem.initial_x), 0))
    od = OnceDifferentiable(f, g!, fill!(similar(problem.initial_x), 0))
    td = TwiceDifferentiable(f, g!, h!, fill!(similar(problem.initial_x), 0))
    tdref = TwiceDifferentiable(f, g!, h!, fill!(similar(problem.initial_x), 0))
    ref = optimize(tdref, problem.initial_x, Newton(), Optim.Options())
    # test AbstractObjective interface
    for obj in (nd, od, td)
        res = []
        push!(res, optimize(obj, problem.initial_x))

        push!(res, optimize(obj, problem.initial_x, Optim.Options()))

        for r in res
            @test norm(Optim.minimum(ref)-Optim.minimum(r)) < 1e-6
        end
    end
    ad_res = optimize(od, problem.initial_x, Newton())
    @test norm(Optim.minimum(ref)-Optim.minimum(ad_res)) < 1e-6
    ad_res2 = optimize(od, problem.initial_x, Newton())
    @test norm(Optim.minimum(ref)-Optim.minimum(ad_res2)) < 1e-6
    # test f, g!, h! interface
    for tup in ((f,), (f, g!), (f, g!, h!))
        fgh_res = []
        push!(fgh_res, optimize(tup..., problem.initial_x))
        for m in (NelderMead(), LBFGS(), Newton())
            push!(fgh_res, optimize(tup..., problem.initial_x; f_tol = 1e-8))
            push!(fgh_res, optimize(tup..., problem.initial_x, m))
            push!(fgh_res, optimize(tup..., problem.initial_x, m, Optim.Options()))
        end
        for r in fgh_res
            @test norm(Optim.minimum(ref)-Optim.minimum(r)) < 1e-6
        end
    end
end


text:
abstract type AbstractPotentialCalculationSetup{T, N} end

"""
    struct PotentialCalculationSetup

`PotentialCalculationSetup` holds the grid, fields and certain precalculated fixed parameters for the field calculation.
This struct will be calculated after each refinement as is depends on the grid.

`grid`: The 3-dimensional grid, either Cartesian or cylindrical, on which the field will be calculated.

The fields:
* `potential`: This is the 4-dimensional array of the extended 3D potential array.
The fourth dimensions comes from the red-black (even-odd) division in order to parallelize the field calculation.
Extended means that the grid holds one additional tick at both sides of all axes necessary for boundary handling (reflecting, fixed, ...) at the ends of the grid.
* `point_types`: Also a 4-dimensional array. Same structure as the `potential`-array.
* `volume_weights`: Also a 4-dimensional array. Same structure as the `potential`-array.
* `q_eff_imp`: Also a 4-dimensional array. Same structure as the `potential`-array.
* `imp_scale`: Also a 4-dimensional array. Same structure as the `potential`-array.
* `q_eff_fix`: Also a 4-dimensional array. Same structure as the `potential`-array.
* `ϵ_r`: Normal 3-dimensional array! In order to calculate the final 6 weights for each grid point in the SOR,
the eight values of the dielectric permittivity (of each octant) around the grid point needs to be loaded.
Here not special division is possible for the red-black (even-odd) division.

Precalculated parameters:
* `geom_weights`: The parts of the calculation of the six weights in the SOR can be precalculated. Those are stored here for each axis/dimension.
* `sor_const`: Vector holding the SOR constants. In the cartesian case only the first entry is used. As the optimal value for the SOR constant
depends on the grid, the constant is linear increased and the array holds the respective value for each radial axis tick.
* `bias_voltage`: `maximum_applied_potential - minimum_applied_potential`. Used for depletion handling, but might be obsolete by now.
* `maximum_applied_potential`: Used for depletion handling, but might be obsolete by now.
* `minimum_applied_potential`: Used for depletion handling, but might be obsolete by now.
* `grid_boundary_factors`: Used in the application of boundary conditions in the field calculation for decaying (infinite) boundary conditions
to approximate the decay of the potential (depending on the grid).
"""
struct PotentialCalculationSetup{
            T, S, N1, DATN1<:AbstractArray{T,N1},
            N2, AT,
            DATN2<:AbstractArray{T,N2},
            DATPT<:AbstractArray{PointType,N2},
            DATGW<:AbstractArray{T,2},
            DATSOR<:AbstractArray{T,1}
        } <: AbstractPotentialCalculationSetup{T, N1}
    grid::Grid{T, N1, S, AT}
    potential::DATN2 # Array{T, N2} or e.g. CuArray{T, N2}
    point_types::DATPT # Array{PointType, N2} or ...
    volume_weights::DATN2 # Array{T, N2}
    q_eff_imp::DATN2 # Array{T, N2}
    imp_scale::DATN2 # Array{T, N2}
    q_eff_fix::DATN2 # Array{T, N2}
    ϵ_r::DATN1 # Array{T, N1}
    geom_weights::NTuple{3, DATGW} # NTuple{3, Array{T, 2}}
    sor_const::DATSOR # Array{T, 1}
    bias_voltage::T
    maximum_applied_potential::T
    minimum_applied_potential::T
    grid_boundary_factors::NTuple{3, NTuple{2, T}}
end

function Adapt.adapt_structure(to, pcs::PotentialCalculationSetup{T, S, 3}) where {T, S}
    PotentialCalculationSetup(
        pcs.grid,
        adapt(to, pcs.potential),
        adapt(to, pcs.point_types),
        adapt(to, pcs.volume_weights),
        adapt(to, pcs.q_eff_imp),
        adapt(to, pcs.imp_scale),
        adapt(to, pcs.q_eff_fix),
        adapt(to, pcs.ϵ_r),
        adapt(to, pcs.geom_weights),
        adapt(to, pcs.sor_const),
        pcs.bias_voltage,
        pcs.maximum_applied_potential,
        pcs.minimum_applied_potential,
        pcs.grid_boundary_factors
    )
end

include("BoundaryConditions/BoundaryConditions.jl")
include("PotentialCalculationSetupCylindrical.jl")
include("PotentialCalculationSetupCartesian3D.jl")


text:
## CPU

function _grid2particle(p::NTuple, xi::NTuple, dxi::NTuple, F::Array)
    # check that the particle is inside the grid
    isinside(p, xi)

    # indices of lowermost-left corner of the cell
    # containing the particle
    idx = parent_cell(p, dxi)

    # normalize particle coordinates
    ti = normalize_coordinates(p, xi, dxi, idx)

    # F at the cell corners
    Fi = field_corners(F, idx)

    # Interpolate field F onto particle
    Fp = ndlinear(ti, Fi)

    return Fp
end

function grid2particle(xi, F::Array{T,N}, particle_coords) where {T,N}
    Fp = zeros(T, np)

    # cell dimensions
    dxi = grid_size(xi)

    Threads.@threads for i in eachindex(particle_coords[1])
        @inbounds Fp[i] = _grid2particle(
            ntuple(j -> particle_coords[j][i], Val(N)), xi, dxi, F
        )
    end

    return Fp
end

function grid2particle!(Fp, xi, F::Array{T,N}, particle_coords) where {T,N}
    # cell dimensions
    dxi = grid_size(xi)

    Threads.@threads for i in eachindex(particle_coords[1])
        @inbounds Fp[i] = _grid2particle(
            ntuple(j -> particle_coords[j][i], Val(N)), xi, dxi, F
        )
    end
end

## CUDA

function _grid2particle!(
    Fp::CuDeviceVector{T,1},
    p::NTuple{N,CuDeviceVector{T,1}},
    dxi::NTuple{N,T},
    xi::NTuple{N,A},
    F::CuDeviceArray{T,N},
    n::Integer,
) where {T,A,N}

    ix = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    @inbounds if ix ≤ n
        pix = particle2tuple(p, ix)

        # check that the particle is inside the grid
        isinside(pix, xi)

        # indices of lowermost-left corner of the cell
        # containing the particle
        idx = parent_cell(pix, dxi)

        # normalize particle coordinates
        ti = normalize_coordinates(pix, xi, dxi, idx)

        # F at the cell corners
        Fi = field_corners(F, idx)

        # Interpolate field F onto particle
        Fp[ix] = ndlinear(ti, Fi)
    end

    return nothing
end

function grid2particle(
    xi, Fd::CuArray{T,N}, particle_coords::NTuple{N,CuArray}; nt=512
) where {T,N}
    dxi = grid_size(xi)
    n = length(particle_coords[1])
    Fpd = CuArray{T,1}(undef, n)

    numblocks = ceil(Int, n / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _grid2particle!(
            Fpd, particle_coords, dxi, xi, Fd, n
        )
    end

    return Fpd
end

function grid2particle!(
    Fpd::CuArray{T,1}, xi, Fd::CuArray{T,N}, particle_coords::NTuple{N,CuArray}; nt=512
) where {T,N}
    dxi = grid_size(xi)
    n = length(particle_coords[1])

    numblocks = ceil(Int, n / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _grid2particle!(
            Fpd, particle_coords, dxi, xi, Fd, n
        )
    end
end


text:
using Base.Broadcast: broadcasted, instantiate
using Test
using MacroTools
using NDReducibles
using NDReducibles: AccessPattern, Index, plan, _axisfor

_indices(idxs::String) = Index.(Tuple(Symbol.(collect(idxs))))

ND(n) = zeros(ntuple(_ -> 0, n))

macro test_inferred(ex)
    ex = quote
        $Test.@test (($Test.@inferred $ex); true)
    end
    esc(ex)
end

macro test_broken_inferred(ex)
    ex = quote
        $Test.@test_broken (($Test.@inferred $ex); true)
    end
    esc(ex)
end


text:
"Replace old symbol with new (expression or symbol) in expression"
function substitute!(e::Expr, old::Symbol, new::Union{Expr, Symbol})
   for (i,a) in enumerate(e.args)
       if a==old
           e.args[i] = new
       elseif a isa Expr
           substitute!(a, old, new)
       end
       ## otherwise do nothing
   end
   e
end

function substitute(e::Expr, old::Symbol, new::Union{Expr, Symbol})
    e_copy = copy(e)
   for (i,a) in enumerate(e_copy.args)
       if a==old
           e_copy.args[i] = new
       elseif a isa Expr
           substitute!(a, old, new)
       end
       ## otherwise do nothing
   end
   e_copy
end

"For example: (x = y, y = 5*2, z = y / x) --> (y = 5*2, x = y, z = y / x)"
function get_expression_order(vargs::Expr...)
    n_expressions = length(vargs)
    exprs_left = collect(vargs)
    exprs_placed = Array{Expr, 1}()
    while length(exprs_left) > 0
        if length(exprs_left) == 1
            push!(exprs_placed, exprs_left[1])
            break
        end

        found_independent = false
        for (i, expr) in enumerate(exprs_left)

            return_symbols = Array{Symbol, 1}()
            for expr_other in filter(x -> x != expr, exprs_left)
                if expr_other.head == :(=)
                    if expr_other.args[1] isa Expr
                        # expr_other as "(x, y) = ..."
                        push!(return_symbols, expr_other.args[1].args...)
                    else
                        # expr_other as "x = ..."
                        push!(return_symbols, expr_other.args[1])
                    end
                end
            end

            expr_rightside = is_assignment(expr) ? expr.args[2] : expr # Check if like "x = ..." or "y * 2"
            if expr_rightside isa Symbol
                # expr_rightside as format "y"
                has_no_dependency = all(
                    expr_rightside == return_symbol
                    for return_symbol in return_symbols
                )
            else
                # expr_rightside as format "y * 2"
                has_no_dependency = all(
                    ~occursin(return_symbol, expr_rightside)
                    for return_symbol in return_symbols
                )
            end

            if has_no_dependency
                push!(exprs_placed, expr)
                deleteat!(exprs_left, i)
                found_independent = true
                break
            end
        end
        if ~found_independent
            throw(ArgumentError("Cross reference in expressions: $exprs_left."))
        end
    end
    exprs_placed
end

is_assignment(e::Expr) = e.head == :(=)
import Base: occursin

function Base.occursin(s::Symbol, e::Expr)

    for arg in e.args
        if s == arg
            return true
        elseif typeof(arg) == Expr
            if occursin(s, arg)
                 return true
            end
        end
    end
    return false
end


text:
function accurate_cosx(x)
    an = 1
    sqrX = x^2
    sn = 0
    m = 2
    while abs(an) > eps()
        sn += an
        an = - an * sqrX /(m-1)/m
        m += 2
    end
    return sn
end


text:
using Helloworld
using Test

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