text:
#Sparse Ensemble Kalman Inversion: specific structures and function definitions
using Convex, SCS
using SparseArrays

"""
    SparseInversion <: Process

A sparse ensemble Kalman Inversion process

# Fields

$(TYPEDFIELDS)

# Constructors

$(METHODLIST)
"""
Base.@kwdef struct SparseInversion{FT <: AbstractFloat} <: Process
    "upper limit of l1-norm"
    γ::FT
    "threshold below which the norm of parameters is pruned to zero"
    threshold_value::FT = FT(0)
    "indices of parameters included in the evaluation of l1-norm constraint"
    uc_idx::Union{AbstractVector, Colon} = Colon()
    "a small regularization value to enhance robustness of convex optimization"
    reg::FT = FT(0)
end

function SparseInversion(
    γ::FT;
    threshold_value::FT = FT(0),
    uc_idx::Union{AbstractVector, Colon} = Colon(),
    reg::FT = FT(0),
) where {FT <: AbstractFloat}
    return SparseInversion{FT}(γ, threshold_value, uc_idx, reg)
end

function FailureHandler(process::SparseInversion, method::IgnoreFailures)
    failsafe_update(ekp, u, g, y, obs_noise_cov, failed_ens) = sparse_eki_update(ekp, u, g, y, obs_noise_cov)
    return FailureHandler{SparseInversion, IgnoreFailures}(failsafe_update)
end

"""
    FailureHandler(process::SparseInversion, method::SampleSuccGauss)

Provides a failsafe update that
 - updates the successful ensemble according to the SparseEKI update,
 - updates the failed ensemble by sampling from the updated successful ensemble.
"""
function FailureHandler(process::SparseInversion, method::SampleSuccGauss)
    function failsafe_update(ekp, u, g, y, obs_noise_cov, failed_ens)
        successful_ens = filter(x -> !(x in failed_ens), collect(1:size(g, 2)))
        n_failed = length(failed_ens)
        u[:, successful_ens] =
            sparse_eki_update(ekp, u[:, successful_ens], g[:, successful_ens], y[:, successful_ens], obs_noise_cov)
        if !isempty(failed_ens)
            u[:, failed_ens] = sample_empirical_gaussian(u[:, successful_ens], n_failed)
        end
        return u
    end
    return FailureHandler{SparseInversion, SampleSuccGauss}(failsafe_update)
end

"""
    sparse_qp(
        ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
        v_j::Vector{FT},
        cov_vv_inv::AbstractMatrix{FT},
        H_u::SparseArrays.SparseMatrixCSC{FT},
        H_g::SparseArrays.SparseMatrixCSC{FT},
        y_j::Vector{FT};
        H_uc::SparseArrays.SparseMatrixCSC{FT} = H_u,
    ) where {FT, IT}

Solving quadratic programming problem with sparsity constraint.
"""
function sparse_qp(
    ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
    v_j::Vector{FT},
    cov_vv_inv::AbstractMatrix{FT},
    H_u::AbstractMatrix{FT},
    H_g::AbstractMatrix{FT},
    y_j::Vector{FT};
    H_uc::AbstractMatrix{FT} = H_u,
) where {FT, IT}

    P = H_g' * (ekp.obs_noise_cov \ H_g) + cov_vv_inv
    P = 0.5 * (P + P')
    q = -(cov_vv_inv * v_j + H_g' * (ekp.obs_noise_cov \ y_j))
    N_params = size(H_uc)[1]
    P1 = vcat(
        hcat(P, fill(FT(0), size(P)[1], N_params)),
        hcat(fill(FT(0), N_params, size(P)[1]), fill(FT(0), N_params, N_params)),
    )
    q1 = vcat(q, fill(FT(0), N_params, 1))
    H_uc_abs = 1.0 * I(N_params)
    G = hcat(vcat(H_uc, -1.0 * H_uc), vcat(-1.0 * H_uc_abs, -1.0 * H_uc_abs))
    h = fill(FT(0), 2 * N_params, 1)
    G1 = vcat(G, hcat(fill(FT(0), 1, size(P)[1]), fill(FT(1), 1, N_params)))
    h1 = vcat(h, ekp.process.γ)
    x = Variable(size(P1)[1])
    problem = minimize(0.5 * quadform(x, P1; assume_psd = true) + q1' * x, [G1 * x <= h1])
    solve!(problem, () -> SCS.Optimizer(verbose = false))

    return hcat(H_u, fill(FT(0), size(H_u)[1], N_params)) * evaluate(x)
end

"""
     sparse_eki_update(
        ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
        u::AbstractMatrix{FT},
        g::AbstractMatrix{FT},
        y::AbstractMatrix{FT},
        obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}},
    ) where {FT <: Real, IT}

Returns the sparse updated parameter vectors given their current values and
the corresponding forward model evaluations, using the inversion algorithm
from eqns. (3.7) and (3.14) of Schneider et al. (2021).

Localization is applied following Tong and Morzfeld (2022).
"""
function sparse_eki_update(
    ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
    u::AbstractMatrix{FT},
    g::AbstractMatrix{FT},
    y::AbstractMatrix{FT},
    obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}},
) where {FT <: Real, IT}

    cov_ug = cov(u, g, dims = 2, corrected = false) # [N_par × N_obs]
    cov_gg = cov(g, g, dims = 2, corrected = false) # [N_par × N_obs]
    cov_uu = cov(u, u, dims = 2, corrected = false) # [N_par × N_par]

    # Localization following Tong and Morzfeld (2022)
    cov_ug = cov_ug .* ekp.localizer.kernel

    v = hcat(u', g')

    # N_obs × N_obs \ [N_obs × N_ens]
    # --> tmp is [N_obs × N_ens]
    tmp = (cov_gg + obs_noise_cov) \ (y - g)
    u = u + (cov_ug * tmp) # [N_par × N_ens]

    # Sparse EKI
    cov_vv = vcat(hcat(cov_uu, cov_ug), hcat(cov_ug', cov_gg))
    H_u = hcat(1.0 * I(size(u)[1]), fill(FT(0), size(u)[1], size(g)[1]))
    H_g = hcat(fill(FT(0), size(g)[1], size(u)[1]), 1.0 * I(size(g)[1]))

    H_uc = H_u[ekp.process.uc_idx, :]

    cov_vv_inv = cov_vv \ (1.0 * I(size(cov_vv)[1]))

    # Loop over ensemble members to impose sparsity
    Threads.@threads for j in 1:size(u, 2)
        # Solve a quadratic programming problem
        u[:, j] = sparse_qp(ekp, v[j, :], cov_vv_inv, H_u, H_g, y[:, j], H_uc = H_uc)

        # Prune parameters using threshold
        u[ekp.process.uc_idx, j] =
            u[ekp.process.uc_idx, j] .* (abs.(u[ekp.process.uc_idx, j]) .> ekp.process.threshold_value)

        # Add small noise to constrained elements of u
        if isposdef(ekp.process.reg)
            len_u_sparse = length(u[ekp.process.uc_idx, j])
            u[ekp.process.uc_idx, j] += rand(ekp.rng, MvNormal(zeros(len_u_sparse), ekp.process.reg * I(len_u_sparse)))
        end
    end
    return u
end

"""
    update_ensemble!(
        ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
        g::AbstractMatrix{FT};
        cov_threshold::FT = 0.01,
        Δt_new = nothing,
        deterministic_forward_map = true,
        failed_ens = nothing,
    ) where {FT, IT}

Updates the ensemble according to a SparseInversion process.

Inputs:
 - ekp :: The EnsembleKalmanProcess to update.
 - g :: Model outputs, they need to be stored as a `N_obs × N_ens` array (i.e data are columms).
 - cov_threshold :: Threshold below which the reduction in covariance determinant results in a warning.
 - Δt_new :: Time step to be used in the current update.
 - deterministic_forward_map :: Whether output `g` comes from a deterministic model.
 - failed_ens :: Indices of failed particles. If nothing, failures are computed as columns of `g`
    with NaN entries.
"""
function update_ensemble!(
    ekp::EnsembleKalmanProcess{FT, IT, SparseInversion{FT}},
    g::AbstractMatrix{FT};
    cov_threshold::FT = 0.01,
    Δt_new = nothing,
    deterministic_forward_map = true,
    failed_ens = nothing,
) where {FT, IT}

    #catch works when g non-square
    if !(size(g)[2] == ekp.N_ens)
        throw(
            DimensionMismatch(
                "ensemble size in EnsembleKalmanProcess and g do not match, try transposing g or check ensemble size",
            ),
        )
    end

    # u: N_par × N_ens
    # g: N_obs × N_ens
    u = get_u_final(ekp)
    N_obs = size(g, 1)
    cov_init = cov(u, dims = 2)
    set_Δt!(ekp, Δt_new)
    fh = ekp.failure_handler

    # Scale noise using Δt
    scaled_obs_noise_cov = ekp.obs_noise_cov / ekp.Δt[end]
    noise = rand(ekp.rng, MvNormal(zeros(N_obs), scaled_obs_noise_cov), ekp.N_ens)

    # Add obs_mean (N_obs) to each column of noise (N_obs × N_ens) if
    # G is deterministic
    y = deterministic_forward_map ? (ekp.obs_mean .+ noise) : (ekp.obs_mean .+ zero(noise))

    if isnothing(failed_ens)
        _, failed_ens = split_indices_by_success(g)
    end
    if !isempty(failed_ens)
        @info "$(length(failed_ens)) particle failure(s) detected. Handler used: $(nameof(typeof(fh).parameters[2]))."
    end

    u = fh.failsafe_update(ekp, u, g, y, scaled_obs_noise_cov, failed_ens)

    # store new parameters (and model outputs)
    push!(ekp.u, DataContainer(u, data_are_columns = true))
    push!(ekp.g, DataContainer(g, data_are_columns = true))

    # Store error
    compute_error!(ekp)

    # Check convergence
    cov_new = cov(get_u_final(ekp), dims = 2)
    cov_ratio = det(cov_new) / det(cov_init)
    if cov_ratio < cov_threshold
        @warn string(
            "New ensemble covariance determinant is less than ",
            cov_threshold,
            " times its previous value.",
            "\nConsider reducing the EK time step.",
        )
    end
end


text:
#=
This script will create code for the Jacobian of the finite volume discretization. The
code is used to calculate Jacobians. We have four functions that get created:
density_flux!             The density part of the flux differentiated w.r.t. P_ij
density_coupling!         The density part of the flux differentiated w.r.t. T_ij
temperature_coupling!     The temperature part of the flux differentiated w.r.t. P_ij
temperature_flux!         The temperature part of the flux differentiated w.r.t. T_ij

The finite volume discretization is local, so the flux at a given point (i,j) only depends
on the variables at that point and the points surrounding it. We use a five point stencil
of the form:
            u_i,j-1
    u_i-1,j u_ij   u_i+1j
            u_ij+1
Throughout the symbolic section, we will denote points on this stencil as follows:
         u_21
    u_12 u_22 u_23
         u_32
We will have to calculate the flux of the probability distribution and the flux of the
temperature distribution. Both of these depend on the probability distribution, the
temperature and the potential. We have the following notations for discretised variables:

P_ij:    The discrete probability distribution at stencil point i,j.
T_ij:    The discrete temperature at stencil point i,j.
T_x_ij:  The x gradient of the temperature at the left edge of the i,j cell.
T_y_ij:  The y gradient of the temperature at the bottom edge of the i,j cell.
V_x_ij:  The x derivative of the potential at the left edge of the i,j cell.
V_y_ij:  The y derivative of the potential at the bottom edge of the i,j cell.
Vxs_ij:  The potential at the left edge of the i,j cell (i.e. the x shifted potential).
Vys_ij:  The potential at the bottom edge of the i,j cell (i.e. the y shifted potential).
=#

using SymPy
# Create the symbolic variables.
@syms(P_12,P_21,P_22,P_32,P_23,T_12,T_21,T_23,T_32,T_22,
      V_x_12,V_x_21,V_x_32,V_x_23,V_x_22,
      V_y_12,V_y_21,V_y_32,V_y_23,V_y_22,V_22,real=true)
@syms(Vxs_ij,Vxs_ip1j,Vys_ij,Vys_ijp1,Jx_ij,Jx_ip1j,
      Jy_ij,Jy_ip1j,Jy_ijp1,dx,dy,A,C,real=true)

Jx_ij = -((P_12+P_22)*V_x_22+(T_12+T_22)*(P_22-P_12)/dx)/2
Jx_ip1j = -((P_22+P_32)*V_x_32+(T_22+T_32)*(P_32-P_22)/dx)/2

Jy_ij = -((P_21+P_22)*V_y_22+(T_21+T_22)*(P_22-P_21)/dy)/2
Jy_ijp1 = -((P_22+P_23)*V_y_23+(T_22+T_23)*(P_23-P_22)/dy)/2

T_x_ij = (T_22-T_12)/dx
T_x_ip1j = (T_32-T_22)/dx
T_y_ij = (T_22-T_21)/dy
T_y_ijp1 = (T_23-T_22)/dy
# The definitions for the fluxes are based on the fluxes from the function
# `flux!(du,u,params...)`.
density_flux = (Jx_ij-Jx_ip1j)/dx+(Jy_ij-Jy_ijp1)/dy
energy_flux = ((Jx_ij*Vxs_ij-Jx_ip1j*Vxs_ip1j-C*(T_x_ij-T_x_ip1j))/dx
                 +((Jy_ij*Vys_ij-Jy_ijp1*Vys_ijp1-C*(T_y_ij-T_y_ijp1)))/dy)
temperature_flux = A*(energy_flux-density_flux*V_22)
# Now we want to turn the symbolic expressions into Julia code. To do this, we will
# replace the stencil variables with their corresponding array variables.
replacements = [("T_22","Tmat[i,j]"),("T_21","Tmat[i,j-1]"),("T_12","Tmat[i-1,j]"),
                ("T_23","Tmat[i,j+1]"),("T_32","Tmat[i+1,j]"),
                ("P_22","Pmat[i,j]"),("P_21","Pmat[i,j-1]"),("P_12","Pmat[i-1,j]"),
                ("P_23","Pmat[i,j+1]"),("P_32","Pmat[i+1,j]"),
                ("V_x_22","V_x[i,j]"),("V_x_21","V_x[i,j-1]"),("V_x_12","V_x[i-1,j]"),
                ("V_x_23","V_x[i,j+1]"),("V_x_32","V_x[i+1,j]"),
                ("V_y_22","V_y[i,j]"),("V_y_21","V_y[i,j-1]"),("V_y_12","V_y[i-1,j]"),
                ("V_y_23","V_y[i,j+1]"),("V_y_32","V_y[i+1,j]"),
                ("Vxs_ip1j","Vxshift[i+1,j]"),("Vxs_ij","Vxshift[i,j]"),
                ("Vys_ijp1","Vyshift[i,j+1]"),("Vys_ij","Vyshift[i,j]"),
                ("V_22","V[i,j]"),
                ("C","coupling")]
# These stencils correspond to the functions:
#    density_flux!
#    density_coupling!
#    temperature_coupling!
#    temperature_flux!
density_stencil = Array{String}(undef,5)
density_coupling_stencil = Array{String}(undef,5)
temperature_coupling_stencil = Array{String}(undef,5)
temperature_stencil = Array{String}(undef,5)
@syms dP_P dT_P
for (index,var) in enumerate([P_22,P_12,P_21,P_32,P_23])
    dP_P = diff(density_flux,var) |> simplify |> string
    dT_P = diff(temperature_flux,var) |> simplify |> string
    for replacement in replacements
        dP_P = replace(dP_P,replacement[1]=>replacement[2])
        dT_P = replace(dT_P,replacement[1]=>replacement[2])
    end
    dP_P = replace(dP_P," "=>"")
    dT_P = replace(dT_P," "=>"")
    density_stencil[index] = dP_P
    temperature_coupling_stencil[index] = dT_P
end
@syms dP_T dT_T
for (index,var) in enumerate([T_22,T_12,T_21,T_32,T_23])
    dP_T = diff(density_flux,var) |> simplify |> string
    dT_T = diff(temperature_flux,var) |> simplify |> string
    for replacement in replacements
        dP_T = replace(dP_T,replacement[1]=>replacement[2])
        dT_T = replace(dT_T,replacement[1]=>replacement[2])
    end
    dP_T = replace(dP_T," "=>"")
    dT_T = replace(dT_T," "=>"")
    density_coupling_stencil[index] = dP_T
    temperature_stencil[index] = dT_T
end

function_names = ["density_flux!","density_coupling!",
                  "temperature_coupling!","temperature_flux!"]
stencils = [density_stencil,density_coupling_stencil,
            temperature_coupling_stencil,temperature_stencil]
for i in 1:4
    function_name = Meta.parse(function_names[i])
    stencil = Meta.parse.(stencils[i])

    quote
        function ($function_name)(::Type{Val{:jac}},jac,u,params)
            @unpack Pmat,Tmat,jac_,Jx,Jy,V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling = params
            nx = length(axes(Tmat)[1])-2
            ny = length(axes(Tmat)[2])-2
            Pmat[1:nx,1:ny] .= reshape(u[1:nx*ny],nx,ny)
            # Periodicity in the x direction.
            Pmat[0,1:ny] .= Pmat[nx-1,1:ny]
            Pmat[nx+1,1:ny] .= Pmat[2,1:ny]

            Tmat[1:nx,1:ny] .= reshape(u[(nx*ny+1):end],nx,ny)
            # Periodicity in the x direction.
            Tmat[0,1:ny] .= Tmat[nx-1,1:ny]
            Tmat[nx+1,1:ny] .= Tmat[2,1:ny]

            inds = Array{Int64}(undef,5)
            i2s = CartesianIndices((nx,ny))
            linear_indices = LinearIndices(jac)
            diag0 = Vector{eltype(jac)}(undef,size(jac)[1])
            diagm1 = Vector{eltype(jac)}(undef,size(jac)[1])
            diagm1_inds = Vector{Int}(undef,size(jac)[1])
            diagmnx = Vector{eltype(jac)}(undef,size(jac)[1])
            diagp1 = Vector{eltype(jac)}(undef,size(jac)[1])
            diagp1_inds = Vector{Int}(undef,size(jac)[1])
            diagpnx = Vector{eltype(jac)}(undef,size(jac)[1])
            for row in 1:nx*ny
                i_j = i2s[row]
                i,j = i_j[1],i_j[2]
                stencil_indices!(inds,(nx,ny),row)
                # jac[row,inds] .= [$(stencil[1]),$(stencil[2]),$(stencil[3]),
                #                   $(stencil[4]),$(stencil[5])]
                diag0[row] = $(stencil[1])
                diagm1[row] = $(stencil[2])
                diagm1_inds[row] = linear_indices[row,inds[2]]
                diagmnx[row] = $(stencil[3])
                diagp1[row] = $(stencil[4])
                diagp1_inds[row] = linear_indices[row,inds[4]]
                diagpnx[row] = $(stencil[5])
            end
            jac[diagind(jac,0)] .= diag0
            jac[diagm1_inds] .= diagm1
            jac[diagp1_inds] .= diagp1
            jac[[diagind(jac,nx*ny-nx);diagind(jac,-nx)]] .= diagmnx
            jac[[diagind(jac,nx);diagind(jac,-nx*ny+nx)]] .= diagpnx
            jac[diagind(jac,nx*ny-1)] .= zero(eltype(jac))
            jac[diagind(jac,-nx*ny+1)] .= zero(eltype(jac))
            jac[diagind(jac,nx*ny-nx)] .= zero(eltype(jac))
            jac[diagind(jac,-nx*ny+nx)] .= zero(eltype(jac))
            jac[diagind(jac,nx*ny-1)] .= zero(eltype(jac))
            jac[diagind(jac,-nx*ny+1)] .= zero(eltype(jac))
            jac[diagind(jac,nx*ny-nx)] .= zero(eltype(jac))
            jac[diagind(jac,-nx*ny+nx)] .= zero(eltype(jac))
            jac
        end
    end |> eval

    eval(:(($function_name)(::Type{Val{:jac}},jac,u,params,t) = ($function_name)(Val{:jac},jac,u,params)))
end

ext: jl


text:
export AbnormalTerminationNtf, AgentLifecycleNtf, AgentStartNtf, AgentTerminationNtf
export FailureNtf, RefuseRsp, CapabilityListRsp, CapabilityReq, ParamChangeNtf
export DatagramReq, DatagramNtf, ClearReq, DatagramCancelReq
export DatagramDeliveryNtf, DatagramFailureNtf, DatagramProgressNtf
export TxFrameReq, TxRawFrameReq, TxFrameStartNtf, TxFrameNtf, RxFrameStartNtf, RxFrameNtf
export BadFrameNtf, CollisionNtf, TxJanusFrameReq, RxJanusFrameNtf, FecDecodeReq
export TxBasebandSignalReq, RxBasebandSignalNtf, RecordBasebandSignalReq, GetPreambleSignalReq
export AddressAllocReq, AddressAllocRsp, AddressResolutionReq, AddressResolutionRsp
export RangeReq, BeaconReq, RangeNtf, InterrogationNtf, RespondReq
export ReservationReq, ReservationCancelReq, ReservationAcceptReq, TxAckReq
export ReservationRsp, ReservationStatusNtf, LinkStatusNtf
export EditRouteReq, GetRouteReq, RouteRsp, RouteChangeNtf, RouteTraceReq, RouteTraceNtf
export DatagramTraceReq, RouteDiscoveryReq, RouteDiscoveryNtf
export RemoteTextReq, RemoteFileGetReq, RemoteFilePutReq, RemoteExecReq
export RemoteTextNtf, RemoteFileNtf, RemoteSuccessNtf, RemoteFailureNtf
export AddScheduledSleepReq, RemoveScheduledSleepReq, GetSleepScheduleReq, SleepScheduleRsp, WakeFromSleepNtf
export ClearStateReq, SaveStateReq
export TestReportNtf

global AbnormalTerminationNtf, AgentLifecycleNtf, AgentStartNtf, AgentTerminationNtf
global FailureNtf, RefuseRsp, CapabilityListRsp, CapabilityReq, ParamChangeNtf
global DatagramReq, DatagramNtf, ClearReq, DatagramCancelReq
global DatagramDeliveryNtf, DatagramFailureNtf, DatagramProgressNtf
global TxFrameReq, TxRawFrameReq, TxFrameStartNtf, TxFrameNtf, RxFrameStartNtf, RxFrameNtf
global BadFrameNtf, CollisionNtf, TxJanusFrameReq, RxJanusFrameNtf, FecDecodeReq
global TxBasebandSignalReq, RxBasebandSignalNtf, RecordBasebandSignalReq, GetPreambleSignalReq
global AddressAllocReq, AddressAllocRsp, AddressResolutionReq, AddressResolutionRsp
global RangeReq, BeaconReq, RangeNtf, InterrogationNtf, RespondReq
global ReservationReq, ReservationCancelReq, ReservationAcceptReq, TxAckReq
global ReservationRsp, ReservationStatusNtf, LinkStatusNtf
global EditRouteReq, GetRouteReq, RouteRsp, RouteChangeNtf, RouteTraceReq, RouteTraceNtf
global DatagramTraceReq, RouteDiscoveryReq, RouteDiscoveryNtf
global RemoteTextReq, RemoteFileGetReq, RemoteFilePutReq, RemoteExecReq
global RemoteTextNtf, RemoteFileNtf, RemoteSuccessNtf, RemoteFailureNtf
global AddScheduledSleepReq, RemoveScheduledSleepReq, GetSleepScheduleReq, SleepScheduleRsp, WakeFromSleepNtf
global ClearStateReq, SaveStateReq
global TestReportNtf

function __init__()
  global AbnormalTerminationNtf = MessageClass(@__MODULE__, "org.arl.unet.AbnormalTerminationNtf")
  global AgentLifecycleNtf = MessageClass(@__MODULE__, "org.arl.unet.AgentLifecycleNtf")
  global AgentStartNtf = MessageClass(@__MODULE__, "org.arl.unet.AgentStartNtf")
  global AgentTerminationNtf = MessageClass(@__MODULE__, "org.arl.unet.AgentTerminationNtf")
  global FailureNtf = MessageClass(@__MODULE__, "org.arl.unet.FailureNtf")
  global RefuseRsp = MessageClass(@__MODULE__, "org.arl.unet.RefuseRsp")
  global CapabilityListRsp = MessageClass(@__MODULE__, "org.arl.unet.CapabilityListRsp")
  global CapabilityReq = MessageClass(@__MODULE__, "org.arl.unet.CapabilityReq")
  global ParamChangeNtf = MessageClass(@__MODULE__, "org.arl.unet.ParamChangeNtf")
  global DatagramReq = AbstractMessageClass(@__MODULE__, "org.arl.unet.DatagramReq")
  global DatagramNtf = AbstractMessageClass(@__MODULE__, "org.arl.unet.DatagramNtf")
  global ClearReq = MessageClass(@__MODULE__, "org.arl.unet.ClearReq")
  global DatagramCancelReq = MessageClass(@__MODULE__, "org.arl.unet.DatagramCancelReq")
  global DatagramDeliveryNtf = MessageClass(@__MODULE__, "org.arl.unet.DatagramDeliveryNtf")
  global DatagramFailureNtf = MessageClass(@__MODULE__, "org.arl.unet.DatagramFailureNtf")
  global DatagramProgressNtf = MessageClass(@__MODULE__, "org.arl.unet.DatagramProgressNtf")
  global TxFrameReq = MessageClass(@__MODULE__, "org.arl.unet.phy.TxFrameReq", DatagramReq)
  global TxRawFrameReq = MessageClass(@__MODULE__, "org.arl.unet.phy.TxRawFrameReq")
  global TxFrameStartNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.TxFrameStartNtf")
  global TxFrameNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.TxFrameNtf")
  global RxFrameStartNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.RxFrameStartNtf")
  global RxFrameNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.RxFrameNtf", DatagramNtf)
  global BadFrameNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.BadFrameNtf")
  global CollisionNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.CollisionNtf")
  global TxJanusFrameReq = MessageClass(@__MODULE__, "org.arl.unet.phy.TxJanusFrameReq")
  global RxJanusFrameNtf = MessageClass(@__MODULE__, "org.arl.unet.phy.RxJanusFrameNtf")
  global FecDecodeReq = MessageClass(@__MODULE__, "org.arl.unet.phy.FecDecodeReq")
  BasebandSignal = AbstractMessageClass(@__MODULE__, "org.arl.unet.bb.BasebandSignal")
  global TxBasebandSignalReq = MessageClass(@__MODULE__, "org.arl.unet.bb.TxBasebandSignalReq", BasebandSignal)
  global RxBasebandSignalNtf = MessageClass(@__MODULE__, "org.arl.unet.bb.RxBasebandSignalNtf", BasebandSignal)
  global RecordBasebandSignalReq = MessageClass(@__MODULE__, "org.arl.unet.bb.RecordBasebandSignalReq")
  global GetPreambleSignalReq = MessageClass(@__MODULE__, "org.arl.unet.bb.GetPreambleSignalReq")
  global RangeReq = MessageClass(@__MODULE__, "org.arl.unet.localization.RangeReq")
  global BeaconReq = MessageClass(@__MODULE__, "org.arl.unet.localization.BeaconReq")
  global RangeNtf = MessageClass(@__MODULE__, "org.arl.unet.localization.RangeNtf")
  global InterrogationNtf = MessageClass(@__MODULE__, "org.arl.unet.localization.InterrogationNtf")
  global RespondReq = MessageClass(@__MODULE__, "org.arl.unet.localization.RespondReq")
  global AddressAllocReq = MessageClass(@__MODULE__, "org.arl.unet.addr.AddressAllocReq")
  global AddressAllocRsp = MessageClass(@__MODULE__, "org.arl.unet.addr.AddressAllocRsp")
  global AddressResolutionReq = MessageClass(@__MODULE__, "org.arl.unet.addr.AddressResolutionReq")
  global AddressResolutionRsp = MessageClass(@__MODULE__, "org.arl.unet.addr.AddressResolutionRsp")
  global ReservationReq = MessageClass(@__MODULE__, "org.arl.unet.mac.ReservationReq")
  global ReservationAcceptReq = MessageClass(@__MODULE__, "org.arl.unet.mac.ReservationAcceptReq")
  global ReservationCancelReq = MessageClass(@__MODULE__, "org.arl.unet.mac.ReservationCancelReq")
  global TxAckReq = MessageClass(@__MODULE__, "org.arl.unet.mac.TxAckReq")
  global ReservationRsp = MessageClass(@__MODULE__, "org.arl.unet.mac.ReservationRsp")
  global ReservationStatusNtf = MessageClass(@__MODULE__, "org.arl.unet.mac.ReservationStatusNtf")
  global LinkStatusNtf = MessageClass(@__MODULE__, "org.arl.unet.link.LinkStatusNtf")
  global EditRouteReq = MessageClass(@__MODULE__, "org.arl.unet.net.EditRouteReq")
  global GetRouteReq = MessageClass(@__MODULE__, "org.arl.unet.net.GetRouteReq")
  global RouteRsp = MessageClass(@__MODULE__, "org.arl.unet.net.RouteRsp")
  global RouteChangeNtf = MessageClass(@__MODULE__, "org.arl.unet.net.RouteChangeNtf")
  global RouteTraceReq = MessageClass(@__MODULE__, "org.arl.unet.net.RouteTraceReq")
  global RouteTraceNtf = MessageClass(@__MODULE__, "org.arl.unet.net.RouteTraceNtf")
  global DatagramTraceReq = MessageClass(@__MODULE__, "org.arl.unet.net.DatagramTraceReq")
  global RouteDiscoveryReq = MessageClass(@__MODULE__, "org.arl.unet.net.RouteDiscoveryReq")
  global RouteDiscoveryNtf = MessageClass(@__MODULE__, "org.arl.unet.net.RouteDiscoveryNtf")
  global RemoteTextReq = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteTextReq")
  global RemoteFileGetReq = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteFileGetReq")
  global RemoteFilePutReq = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteFilePutReq")
  global RemoteExecReq = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteExecReq")
  global RemoteTextNtf = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteTextNtf")
  global RemoteFileNtf = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteFileNtf")
  global RemoteSuccessNtf = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteSuccessNtf")
  global RemoteFailureNtf = MessageClass(@__MODULE__, "org.arl.unet.remote.RemoteFailureNtf")
  global AddScheduledSleepReq = MessageClass(@__MODULE__, "org.arl.unet.scheduler.AddScheduledSleepReq")
  global RemoveScheduledSleepReq = MessageClass(@__MODULE__, "org.arl.unet.scheduler.RemoveScheduledSleepReq")
  global GetSleepScheduleReq = MessageClass(@__MODULE__, "org.arl.unet.scheduler.GetSleepScheduleReq")
  global SleepScheduleRsp = MessageClass(@__MODULE__, "org.arl.unet.scheduler.SleepScheduleRsp")
  global WakeFromSleepNtf = MessageClass(@__MODULE__, "org.arl.unet.scheduler.WakeFromSleepNtf")
  global SaveStateReq = MessageClass(@__MODULE__, "org.arl.unet.state.SaveStateReq")
  global ClearStateReq = MessageClass(@__MODULE__, "org.arl.unet.state.ClearStateReq")
  global TestReportNtf = MessageClass(@__MODULE__, "org.arl.unet.TestReportNtf")
end


text:
import Observables
import Makie

export Cylinder
export update

"""
    Represents a cylinder.

fields: q, radius, height, segments, color and observable.
"""
mutable struct Cylinder
    q::Biquaternion
    height::Float64
    radius::Float64
    segments::Int
    color::Observables.Observable{Array{Makie.ColorTypes.RGBA{Float32},2}}
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Cylinder(q, scene)

Construct a Cylinder with the given configuration `q` and `scene`, and the optional
arguments: `height`, `radius`, `segments`, `color` and `transparency`.
"""
function Cylinder(q::Biquaternion,
                  scene::Makie.Scene;
                  height::Float64 = 1.0,
                  radius::Float64 = 0.1,
                  segments::Int = 36,
                  color::Makie.RGBAf = Makie.RGBAf(0.7529,
                                                                           0.7529,
                                                                           0.7529,
                                                                           0.5),
                  transparency::Bool = false)
    cylinder = constructcylinder(q, height, radius, segments = segments)
    colorarray = Observables.Observable(fill(color, segments, segments))
    observable = buildsurface(scene, cylinder, colorarray, transparency = transparency)
    Cylinder(q, height, radius, segments, colorarray, observable)
end

"""
    update(cylinder, q)

Update a Cylinder by changing its observable with the given `cylinder` and configuration `q`
.
"""
function update(cylinder::Cylinder, q::Biquaternion)
    cylinder.q = q
    value = constructcylinder(cylinder.q,
                              cylinder.height,
                              cylinder.radius,
                              segments = cylinder.segments)
    updatesurface(value, cylinder.observable)
end

"""
    update(cylinder, color)

Update a Cylinder by changing its observable with the given `cylinder` and `color`.
"""
function update(cylinder::Cylinder, color::Makie.RGBAf)
    cylinder.color[] = fill(color, cylinder.segments, cylinder.segments)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
abstract AbstractStateSpaceModel
typealias AbstractSSM AbstractStateSpaceModel
abstract AbstractGaussianSSM <: AbstractStateSpaceModel

# LinearGaussianSSM     LinearKalmanFilter
#                       NonlinearKalmanFilter
#                       NonlinearFilter
#
# NonlinearGaussianSSM  NonlinearKalmanFilter
#                       NonlinearFilter
#
# NonlinearSSM          NonlinearFilter

###########################################################################
# Linear Gaussian
###########################################################################

immutable LinearGaussianSSM <: AbstractGaussianSSM
    # Process transition matrix, control matrix, and noise covariance
    F::Function
    B::Function
    V::Function

    # Observation matrix and noise covariance
    G::Function
    W::Function

    # Model dimensions
    nx::Int
    ny::Int
    nu::Int

    function LinearGaussianSSM(F::Function, B::Function, V::Function,
                                G::Function, W::Function)
        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))

        nx, ny, nu = confirm_matrix_sizes(F(1), B(1), V(1), G(1), W(1))
        new(F, B, V, G, W, nx, ny, nu)
    end
end

# Time-dependent constructor
LinearGaussianSSM{T}(F::Function, V::Function, G::Function, W::Function;
                          B::Function=_->zeros(T, size(V(1), 1), 1)) =
	  LinearGaussianSSM(F, B, V, G, W)

# Time-independent constructor
LinearGaussianSSM{T}(F::Matrix{T}, V::Matrix{T}, G::Matrix{T}, W::Matrix{T};
                          B::Matrix{T}=zeros(T, size(F, 1), 1)) =
	  LinearGaussianSSM(_->F, _->B, _->V, _->G, _->W)

function show(io::IO, mod::LinearGaussianSSM)
    dx, dy = mod.nx, mod.ny
    println("LinearGaussianSSM, $dx-D process x $dy-D observations")
    println("Process evolution matrix F:")
    show(mod.F(1))
    println("\n\nGontrol input matrix B:")
    show(mod.B(1))
    println("\n\nProcess error covariance V:")
    show(mod.V(1))
    println("\n\nObservation matrix G:")
    show(mod.G(1))
    println("\n\nObseration error covariance W:")
    show(mod.W(1))
end

## Core methods
process_matrix(m::LinearGaussianSSM, state::Vector, t::Real=0.0) = m.F(t)
process_matrix(m::LinearGaussianSSM, state::AbstractMvNormal, t::Real=0.0) = m.F(t)

observation_matrix(m::LinearGaussianSSM, state::Vector, t::Real=0.0) = m.G(t)
observation_matrix(m::LinearGaussianSSM, state::AbstractMvNormal, t::Real=0.0) = m.G(t)

control_input(m::LinearGaussianSSM, u, t::Real=0.0) = m.B(t) * u

###########################################################################
# Nonlinear Gaussian
###########################################################################

immutable NonlinearGaussianSSM <: AbstractGaussianSSM
    # Process transition function, jacobian, and noise covariance matrix
    f::Function
    V::Function

    # Control input function
    b::Function

    # Observation function,
    g::Function # actual observation function
    W::Function

    nx::Int
    ny::Int
    nu::Int

    function NonlinearGaussianSSM(f::Function, V::Function,
            b::Function, g::Function, W::Function, nx::Int, ny::Int, nu::Int)
        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))
        @assert (nx, nx) == size(V(1))
        @assert (ny, ny) == size(W(1))

        new(f, V, b, g, W, nx, ny, nu)
    end
end

function NonlinearGaussianSSM{T}(f::Function, V::Matrix{T}, g::Function,
        W::Matrix{T}; b::Function=_ -> 0, nu::Int=0)
    nx = size(V, 1)
    ny = size(W, 1)

    return NonlinearGaussianSSM(f, _-> V, b, g, _-> W, nx, ny, nu)
end

## Core methods
process_matrix(m::NonlinearGaussianSSM, x::Vector, t::Real=0.0) = ForwardDiff.jacobian(m.f, x)
function process_matrix(m::NonlinearGaussianSSM, x::AbstractMvNormal, t::Real=0.0)
    return process_matrix(m, mean(x), t)
end

observation_matrix(m::NonlinearGaussianSSM, x::Vector, t::Real=0.0) = ForwardDiff.jacobian(m.g, x)
function observation_matrix(m::NonlinearGaussianSSM, x::AbstractMvNormal, t::Real=0.0)
    return process_matrix(m, mean(x), t)
end

control_input(m::NonlinearGaussianSSM, u, t::Real=0.0) = m.b(u)


text:
#push!(LOAD_PATH, "/home/caseykneale/Desktop/ChemometricsTools/ChemometricsTools.jl/");
using ChemometricsTools
using Plots

#View the data in the package space
ChemometricsToolsDatasets()
#Load the data from package-space: remove the time column while we're at it
DF = ChemometricsToolsDataset("actuator.csv");
actuator = convert(Matrix,  DF)[:,2:end];
#Split into a training set based on normal behaviour
(Train, Test) = SplitByProportion(actuator, 0.40)

#Plot the data
plot(Train, legend = false, title = "Train Set - Raw Data", ylabel = "Measurements", xlabel = "time")
#alot of DC offsets let's column center each variable.
clean = Center(Train);
pca = PCA(clean(Train); Factors = 5);
cumsum( ExplainedVariance( pca ) )
#Old Hotelling Interface
#(Lambdas, limits) = Hotelling(clean(Train), pca; Quantile = 0.015, Variance = 0.95)
T = Hotelling(clean(Train), pca; Quantile = 0.015, Variance = 0.95)
scatter(T(clean(Train)), title = "Train - T2 vs Time", xlabel = "time", ylabel = "T2", label = "T2");
hline!([T.UpperLimit], label = "Control Limit")

#Old Q Interface
#(qs, limits) = Q(clean(Train), pca; Quantile = 1.0 - 1e-16, Variance = 0.99)
q_residuals = Q(clean(Train), pca; Quantile = 1.0 - 1e-16, Variance = 0.99)
scatter(q_residuals(Train), title = "Train - Q vs Time", xlabel = "time", ylabel = "Q", label = "Q");
hline!([T.UpperLimit], label = "Control Limit")

#Let's apply the control limits obtained from the train set to the test set...
scatter(T(clean(Test)), title = "Test - T2 vs Time", xlabel = "time", ylabel = "T2", label = "T2");
hline!([T.UpperLimit], label = "Control Limit");
vline!([120], label = "Known Fault")

scatter(q_residuals(clean(Test)), title = "Test - Q vs Time", xlabel = "time", ylabel = "Q", label = "Q");
hline!([q_residuals.UpperLimit], label = "Control Limit");
vline!([120], label = "Known Fault")

#Hotellings T2 is slower to detect the failure then the Q plot. Not a surprise but good to see.

#Let's try an EWMA on a single variable...
plot(Train[:,3], title = "Variable #1", xlabel = "time", ylabel = "Magnitude")
#Use the "Burn in" constructor for EWMA
EWM = EWMA(Train[:,3], 0.1)
(UpperEWMA, LowerEWMA) = Limits(EWM; k = 3)

#Let's get to predicting already...
ewmascore = []
for s in 1:length( Test[:,3] )
    push!( ewmascore, EWM( Test[s,3] ; train = false) )
end
plot(ewmascore, label = "EWMA", title = "EWMA Fault detection", ylabel = "EWMA", xlabel = "time");
hline!([UpperEWMA], label = "Upper Limit");
hline!([LowerEWMA], label = "Lower Limit");
vline!([120], label = "Known Fault")

score = OneClassJKNN( Train, Test; J = 1, K = 1,  DistanceType = "euclidean" );
scatter(log10.(score));
hline!([1.0], label = "Threshold of 1");
vline!([120], label = "Known Fault")


text:
include("ex2_multiprony.jl")
include("ex3_decompose.jl")
include("ex4_newton.jl")


text:
@testset "Integrals" begin
    @reset

    @set {
        printstyle none
        jkfit auto
        rifit auto
        basis cc-pvtz
    }

    jkfit = Fermi.Integrals.JKFIT()
    rifit = Fermi.Integrals.RIFIT()

    # Test auto-complete
    @test jkfit.basisset.basis_name == "cc-pvtz-jkfit"
    @test rifit.basisset.basis_name == "cc-pvtz-rifit"

    # Test displays
    x = Fermi.Integrals.string_repr(jkfit)
    @test occursin(r"JKFIT:\s+?cc-pvtz-jkfit",x)
    x = Fermi.Integrals.string_repr(rifit)
    @test occursin(r"RIFIT:\s+?cc-pvtz-rifit",x)

    ints = Fermi.Integrals.IntegralHelper()

    # Try invalid precision
    @set precision invalid
    @test_throws Fermi.Options.FermiException Fermi.Integrals.IntegralHelper()

    @set precision single
    ints32 = Fermi.Integrals.IntegralHelper()

    # Test if single and double are approx equal
    @test isapprox(ints["S"], ints32["S"])
    @test isapprox(ints["T"], ints32["T"])
    @test isapprox(ints["V"], ints32["V"])
    @test isapprox(ints["ERI"], ints32["ERI"])

    # Test string_repr
    x = Fermi.Integrals.string_repr(ints)
    @test begin
        occursin(r"Fermi\s+?IntegralHelper",x) &&
        occursin(r"Data\s+?Type:\s+?Float64",x) &&
        occursin(r"Basis:\s+?cc-pvtz",x) &&
        occursin(r"ERI:\s+?SparseERI",x) &&
        occursin(r"Orbitals:\s+?AtomicOrbitals",x)
    end

    # Test cache
    @test begin
        "S" in keys(ints.cache) &&
        "T" in keys(ints.cache) &&
        "V" in keys(ints.cache) &&
        "ERI" in keys(ints.cache)
    end

    # Test two basis set input
    bf = ints.orbitals.basisset
    S = Fermi.Integrals.ao_1e(bf, bf, "overlap")
    T = Fermi.Integrals.ao_1e(bf, bf, "kinetic")
    V = Fermi.Integrals.ao_1e(bf, bf, "nuclear")

    @test S ≈ ints["S"]
    @test T ≈ ints["T"]
    @test V ≈ ints["V"]

    # Test delete one key
    delete!(ints, "V")
    @test begin
        "S" in keys(ints.cache) &&
        "T" in keys(ints.cache) &&
      !("V" in keys(ints.cache)) &&
        "ERI" in keys(ints.cache)
    end
    # Test deleting everything
    delete!(ints)
    @test isempty(ints.cache)

    # Test invalid key
    @test_throws Fermi.Options.FermiException ints["invalid"]

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: condenser.mso 555 2008-07-18 19:01:13Z rafael $
#*--------------------------------------------------------------------
type condenser2
	condenser2()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Initial_Temperature),
			:(OutletLiquid.z(1:NComp-1) = Initial_Composition(1:NComp-1)/sum(Initial_Composition)),
		],
		[
			"Initial level Percent","Initial Temperature","Initial Composition",
		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Low Flow",
				:Default=>1E-6,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"No Flow",
				:Default=>0,
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Condenser type",
				:Valid=>["partial", "total"],
				:Default=>"partial"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Vapour Flow",
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"K factor for pressure drop",
				:Lower=>1E-8,
				:Default=>1E-3
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Temperature of Condenser"
			)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Composition",
				:Lower=>1E-6
			)),(NComp)),
			VesselVolume(Dict{Symbol,Any}(
				:Brief=>"Vessel Geometry",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour inlet stream",
				:PosX=>0.13,
				:PosY=>0,
				:Symbol=>"_{in}^{Vapour}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.35,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:Hidden=>true,
				:Symbol=>"_{out}^{Vapour}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour or Liquid product stream",
				:PosX=>0.54,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat supplied",
				:Protected=>true,
				:PosX=>1,
				:PosY=>0.08,
				:Symbol=>"Q_{in}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Reflux Ratio",
				:Default=>10,
				:Lower=>0.05
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.33,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.43,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.25,
				:PosY=>0
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray",
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Inlet Vapour Mass Density",
				:Hidden=>true,
				:Symbol=>"\\rho"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P",
				:Protected=>true
			)),
			[
				:(Vapour.F = zero_flow),
				:(RefluxRatio*Product.F = OutletLiquid.F),
				:(Product.v = Vapour.v),
				:(Product.h = Vapour.h),
				:(Product.z = Vapour.z),
				:(Product.v = OutletLiquid.v),
				:(Product.h = OutletLiquid.h),
				:(Product.z = OutletLiquid.z),
				:(InletVapour.F*vV = Kfactor *sqrt(Pdrop/rho)*"m^2"),
				:(InletVapour.F = zero_flow),
				:(diff(M) = InletVapour.F*InletVapour.z - OutletLiquid.F*OutletLiquid.z- Product.F*Product.z),
				:(diff(E) = InletVapour.F*InletVapour.h - OutletLiquid.F*OutletLiquid.h- Product.F*Product.h + InletQ),
				:(M = ML*OutletLiquid.z + MV*Vapour.z),
				:(E = ML*OutletLiquid.h + MV*Vapour.h - Vapour.P*Geometry.Vtotal),
				:(sum(OutletLiquid.z)=1.0),
				:(sum(OutletLiquid.z)=sum(Vapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(Vapour.T, Vapour.P, Vapour.z)),
				:(rho = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(Vapour.T, Vapour.P, Vapour.z)*Vapour.z),
				:(OutletLiquid.T = Vapour.T),
				:(OutletLiquid.T = Product.T),
				:(Vapour.P = OutletLiquid.P),
				:(Vapour.P = Product.P),
				:(OutletLiquid.P = InletVapour.P - Pdrop),
				:(Geometry.Vtotal = ML*vL + MV*vV),
				:(ML * vL = Geometry.Vfilled),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
				:(LI*Geometry.Vtotal= Geometry.Vfilled),
			],
			[
				"","Reflux ratio","","","","","","","","","Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol fraction normalisation","Mol fraction Constraint","Liquid Volume","Vapour Volume","Inlet Vapour Density","Chemical Equilibrium","Thermal Equilibrium","","Mechanical Equilibrium","","Pressure Drop","Geometry Constraint","Liquid Level","Temperature indicator","Pressure indicator","Level indicator",
			],
			[:PP,:NComp,:Mw,:low_flow,:zero_flow,:CondenserType,:VapourFlow,:Kfactor,:Levelpercent_Initial,:Initial_Temperature,:Initial_Composition,],
			[:Geometry,:InletVapour,:OutletLiquid,:Vapour,:Product,:InletQ,:RefluxRatio,:TI,:LI,:PI,:M,:ML,:MV,:E,:vL,:vV,:rho,:Pdrop,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	low_flow::flow_mol
	zero_flow::flow_mol
	CondenserType::DanaSwitcher
	VapourFlow::DanaSwitcher
	Kfactor::positive
	Levelpercent_Initial::positive
	Initial_Temperature::temperature
	Initial_Composition::Array{positive}
	Geometry::VesselVolume
	InletVapour::stream
	OutletLiquid::liquid_stream
	Vapour::vapour_stream
	Product::stream
	InletQ::power
	RefluxRatio::positive
	TI::control_signal
	LI::control_signal
	PI::control_signal
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	rho::dens_mass
	Pdrop::press_delta
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export condenser2
function set(in::condenser2)
	Mw = PP.MolecularWeight()
	 low_flow = 1E-6 * "kmol/h"
	 zero_flow = 0 * "kmol/h"

end
function initial(in::condenser2)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::condenser2)
	addEquation(1)
	addEquation(2)
	let switch=CondenserType
		if switch=="partial"
			addEquation(3)
			addEquation(4)
			addEquation(5)
		elseif switch=="total"
			addEquation(6)
			addEquation(7)
			addEquation(8)
		end
	end
	let switch=VapourFlow
		if VapourFlow==InletVapour.F < low_flow
			set(switch,"off")
		end
		if VapourFlow==InletVapour.P > OutletLiquid.P
			set(switch,"on")
		end
		if switch=="on"
			addEquation(9)
		elseif switch=="off"
			addEquation(10)
		end
	end
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
	addEquation(27)
	addEquation(28)
	addEquation(29)
	addEquation(30)
end
function atributes(in::condenser2,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Condenser"
	fields[:Brief]="Model of a  dynamic condenser with control."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SPECIFY ==
* the InletVapour stream;
* the outlet flows: Product.F and Reflux.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the condenser model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator of Condenser;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the condenser temperature (OutletLiquid.T);
* Levelpercent_Initial : the condenser liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
condenser2(_::Dict{Symbol,Any})=begin
	newModel=condenser2()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#
# To plot all cutoff functions, uncomment the code at the end of th e
# file and then execute the entire file.
#

using Test, NBodyIPs, BenchmarkTools
using NBodyIPs.Cutoffs
using NBodyIPs: fcut, fcut_d
import JuLIP

r0 = 0.712
rn = 1.234
rc = 2.890

cutoffs = [ CosCut(rc-1.0, rc),
            CosCut2s(r0, r0+0.4, rc-0.9, rc),
            PolyCut(2, rc),
            PolyCutSym(2, rc),
            PolyCut2sA(r0, rn, rc) ]
cnames = ["CosCut", "CosCut2s", "PolyCut", "PolyCutSym", "PolyCut2sA"]

##
@info("Testing Correctness of Derivative Implementations")
xt = [range(r0-0.5, rc+0.5, length=50); r0; rc]
for (C, name) in zip(cutoffs, cnames)
   println("Testing Derivatives of $name")
   f = r -> fcut(C, r)
   df = r -> fcut_d(C, r)
   println(@test JuLIP.Testing.fdtest_R2R(f, df, xt))
   Cd = Dict(C)
   C1 = JuLIP.decode_dict(Cd)
   Cd1 = Dict(Cd)
   println((@test C1 == C))
   println((@test Cd == Cd1))
end

##
@info("Testing Performance of fcut and fcut_d")
@info("One test includes 1_000 fcut/fcut_d calls")
function runN(f, x, dx, N)
   s = 0.0
   for n = 1:N
      x += dx
      s += f(x)
   end
   return s
end

N = 1_000
x0 = r0
dx = (rc+0.3-r0)/N

for (C, name) in zip(cutoffs, cnames)
   println(name)
   print("   fcut  : ")
   @btime runN($(r->fcut(C,r)), $x0, $dx, $N)
   print("   fcut_d: ")
   @btime runN($(r->fcut_d(C,r)), $x0, $dx, $N)
end

##
# using Plots
# xx = range(r0-0.5, rc+0.5, length=200)
# P = plot(; ylims = [-0.2, 1.2])
# for (C, name) in zip(cutoffs, cnames)
#    plot!(xx, fcut.(Ref(C), xx), label=name)
# end
# P = scatter!([r0, rc], [0, 0], label = "r0, rcut")
# P = vline!([rn], c=:black, ls=:dash, label = "rn")
# display(P)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
export @fixture
export yield_fixture

using Base.Meta
import MetaTools

function yield_fixture(args...)
  error("yield_fixture must be in the top level of scope within @fixture function")
end

const EXTRA_ARG = MetaTools.ParsedArgument(:(ecbf47d557eb469c9fc755f8e07f11f7::Function))

macro fixture(ex::Expr)
  local pfunc = MetaTools.ParsedFunction(ex)
  pfunc.args = [EXTRA_ARG, pfunc.args...]
  if pfunc.body.head == :block
    const i = findfirst(v -> isexpr(v, :call) && v.args[1]==:yield_fixture, pfunc.body.args)
    if i>0
      # Copy the call so that we get any arguments
      const call_expr = copy(pfunc.body.args[i])
      call_expr.args[1] = :ecbf47d557eb469c9fc755f8e07f11f7

      pfunc.body.args = Any[pfunc.body.args[1:(i-1)]...,
                              quote
                                try
                                  return $call_expr
                                finally
                                  $(Expr(:block, pfunc.body.args[(i+1):end]...))
                                end
                              end
                           ]
    else
      pfunc.body.args = Any[pfunc.body.args..., :(ecbf47d557eb469c9fc755f8e07f11f7())]
    end
  else
    error("Expected a :block got a $(pfunc.body.head)")
  end

  return esc(MetaTools.emit(pfunc))
end


text:
@doc raw"""
    RecordAction

A `RecordAction` is a small functor to record values.
The usual call is given by `(p,o,i) -> s` that performs the record based on
a [`Problem`](@ref) `p`, [`Options`](@ref) `o` and the current iterate `i`.

By convention `i<=0` is interpreted as "For Initialization only", i.e. only
initialize internal values, but not trigger any record, the same holds for
`i=typemin(Inf)` which is used to indicate `stop`, i.e. that the record is
called from within [`stop_solver!`](@ref) which returns true afterwards.

# Fields (assumed by subtypes to exist)
* `recorded_values` an `Array` of the recorded values.
"""
abstract type RecordAction <: AbstractOptionsAction end

@doc raw"""
    RecordOptions <: Options

append to any [`Options`](@ref) the decorator with record functionality,
Internally a `Dict`ionary is kept that stores a [`RecordAction`](@ref) for
several concurrent modes using a `Symbol` as reference.
The default mode is `:Iteration`, which is used to store information that is recorded during
the iterations. RecordActions might be added to `:Start` or `:Stop` to record values at the
beginning or for the stopping time point, respectively

The original options can still be accessed using the [`get_options`](@ref) function.

# Fields
* `options` – the options that are extended by debug information
* `recordDictionary` – a `Dict{Symbol,RecordAction}` to keep track of all
  different recorded values

# Constructors
    RecordOptions(o,dR)

construct record decorated [`Options`](@ref), where `dR` can be

* a [`RecordAction`](@ref), then it is stored within the dictionary at `:Iteration`
* an `Array` of [`RecordAction`](@ref)s, then it is stored as a
  `recordDictionary`(@ref) within the dictionary at `:All`.
* a `Dict{Symbol,RecordAction}`.
"""
mutable struct RecordOptions{O<:Options,TRD<:NamedTuple} <: Options
    options::O
    recordDictionary::TRD
    function RecordOptions{O}(o::O; kwargs...) where {O<:Options}
        return new{O,typeof(values(kwargs))}(o, values(kwargs))
    end
end
function RecordOptions(o::O, dR::RecordAction) where {O<:Options}
    return RecordOptions{O}(o; Iteration=dR)
end
function RecordOptions(o::O, dR::Dict{Symbol,<:RecordAction}) where {O<:Options}
    return RecordOptions{O}(o; dR...)
end
function RecordOptions(o::O, format::Vector{<:Any}) where {O<:Options}
    return RecordOptions{O}(o; RecordFactory(get_options(o), format)...)
end
function RecordOptions(o::O, s::Symbol) where {O<:Options}
    return RecordOptions{O}(o; Iteration=RecordFactory(get_options(o), s))
end

dispatch_options_decorator(::RecordOptions) = Val(true)

@doc """
    has_record(o::Options)

check whether the [`Options`](@ref)` o` are decorated with
[`RecordOptions`](@ref)
"""
has_record(::RecordOptions) = true
has_record(o::Options) = has_record(o, dispatch_options_decorator(o))
has_record(o::Options, ::Val{true}) = has_record(o.options)
has_record(::Options, ::Val{false}) = false

@doc """
    get_record_options(o::Options)

return the [`RecordOptions`](@ref) among the decorators from the [`Options`](@ref) `o`
"""
get_record_options(o::Options) = get_record_options(o, dispatch_options_decorator(o))
get_record_options(o::Options, ::Val{true}) = get_record_options(o.options)
get_record_options(::Options, ::Val{false}) = error("No Record decoration found")
get_record_options(o::RecordOptions) = o

@doc raw"""
    get_record_action(o::Options, s::Symbol)

return the action contained in the (first) [`RecordOptions`](@ref) decorator within the [`Options`](@ref) `o`.

"""
function get_record_action(o::Options, s::Symbol=:Iteration)
    if haskey(o.recordDictionary, s)
        return o.recordDictionary[s]
    else
        error("No record known for key :$s found")
    end
end
@doc raw"""
    get_record(o::Options, [,s=:Iteration])
    get_record(o::RecordOptions, [,s=:Iteration])

return the recorded values from within the [`RecordOptions`](@ref) `o` that where
recorded with respect to the `Symbol s` as an `Array`. The default refers to
any recordings during an `:Iteration`.

When called with arbitrary [`Options`](@ref), this method looks for the
[`RecordOptions`](@ref) decorator and calls `get_record` on the decorator.
"""
function get_record(o::RecordOptions, s::Symbol=:Iteration)
    return get_record(get_record_action(o, s))
end
function get_record(o::RecordOptions, s::Symbol, i...)
    return get_record(get_record_action(o, s), i...)
end
get_record(o::Options, s::Symbol=:Iteration) = get_record(get_record_options(o), s)

@doc raw"""
    get_record(r::RecordAction)

return the recorded values stored within a [`RecordAction`](@ref) `r`.
"""
get_record(r::RecordAction, i) = r.recorded_values
get_record(r::RecordAction) = r.recorded_values

"""
    get_index(ro::RecordOptions, s::Symbol)
    ro[s]

Get the recorded values for reording type `s`, see [`get_record`](@ref) for details.

    get_index(ro::RecordOptions, s::Symbol, i...)
    ro[s, i...]

Acces the recording type of type `s` and call its [`RecordAction`](@ref) with `[i...]`.
"""
getindex(ro::RecordOptions, s::Symbol) = get_record(ro, s)
getindex(ro::RecordOptions, s::Symbol, i...) = get_record_action(ro, s)[i...]

"""
    record_or_reset!(r,v,i)

either record (`i>0` and not `Inf`) the value `v` within the [`RecordAction`](@ref) `r`
or reset (`i<0`) the internal storage, where `v` has to match the internal
value type of the corresponding Recordaction.
"""
function record_or_reset!(r::RecordAction, v, i::Int)
    if i > 0
        push!(r.recorded_values, deepcopy(v))
    elseif i < 0 # reset if negative
        r.recorded_values = Array{typeof(v),1}()
    end
end

"""
    RecordGroup <: RecordAction

group a set of [`RecordAction`](@ref)s into one action, where the internal [`RecordAction`](@ref)s
act independently, but the results can be collected in a grouped fashion, i.e. tuples per calls of this group.
The enries can be later addressed either by index or semantic Symbols

# Constructors
    RecordGroup(g::Array{<:RecordAction, 1})

construct a group consisting of an Array of [`RecordAction`](@ref)s `g`,

    RecordGroup(g, symbols)

# Examples
    r = RecordGroup([RecordIteration(), RecordCost()])

A RecordGroup to record the current iteration and the cost. The cost can then be accessed using `get_record(r,2)` or `r[2]`.

    r = RecordGroup([RecordIteration(), RecordCost()], Dict(:Cost => 2))

A RecordGroup to record the current iteration and the cost, wich can then be accesed using `get_record(:Cost)` or `r[:Cost]`.

    r = RecordGroup([RecordIteration(), :Cost => RecordCost()])

A RecordGroup identical to the previous constructor, just a little easier to use.
"""
mutable struct RecordGroup <: RecordAction
    group::Array{RecordAction,1}
    indexSymbols::Dict{Symbol,Int}
    function RecordGroup(
        g::Array{<:RecordAction,1}, symbols::Dict{Symbol,Int}=Dict{Symbol,Int}()
    )
        if length(symbols) > 0
            if maximum(values(symbols)) > length(g)
                error(
                    "Index $(maximum(values(symbols))) must not be larger than number of elements ($(length(g))) in this RecordGroup.",
                )
            end
            if minimum(values(symbols)) < 1
                error("Index $(minimum(values(symbols))) nonpositive not allowed.")
            end
        end
        return new(g, symbols)
    end
    function RecordGroup(
        records::Vector{<:Union{<:RecordAction,Pair{Symbol,<:RecordAction}}}
    )
        g = Array{RecordAction,1}()
        si = Dict{Symbol,Int}()
        for i in 1:length(records)
            if records[i] isa RecordAction
                push!(g, records[i])
            else
                push!(g, records[i].second)
                push!(si, records[i].first => i)
            end
        end
        return RecordGroup(g, si)
    end
    RecordGroup() = new(Array{RecordAction,1}(), Dict{Symbol,Int}())
end

function (d::RecordGroup)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    for ri in d.group
        ri(p, o, i)
    end
end
@doc raw"""
    get_record(r::RecordGroup)

return an array of tuples, where each tuple is a recorded set, e.g. per iteration / record call.

    get_record(r::RecordGruop, i::Int)

return an array of values corresponding to the `i`th entry in this record group

    get_record(r::RecordGruop, s::Symbol)

return an array of recorded values with respect to the `s`, see [`RecordGroup`](@ref).

    get_record(r::RecordGroup, s1::Symbol, s2::Symbol,...)

return an array of tuples, where each tuple is a recorded set corresponding to the symbols `s1, s2,...` per iteration / record call.
"""
get_record(r::RecordGroup) = [zip(get_record.(r.group)...)...]
get_record(r::RecordGroup, i) = get_record(r.group[i])
get_record(r::RecordGroup, s::Symbol) = get_record(r.group[r.indexSymbols[s]])
function get_record(r::RecordGroup, s::NTuple{N,Symbol}) where {N}
    inds = getindex.(Ref(r.indexSymbols), s)
    return [zip(get_record.([r.group[i] for i in inds])...)...]
end

@doc raw"""
    getindex(r::RecordGroup, s::Symbol)
    r[s]
    getindex(r::RecordGroup, sT::NTuple{N,Symbol})
    r[sT]
    getindex(r::RecordGroup, i)
    r[i]

return an array of recorded values with respect to the `s`, the symbols from the tuple `sT` or the index `i`.
See [`get_record`](@ref get_record(r::RecordGroup)) for details.
"""
getindex(::RecordGroup, ::Any...)
getindex(r::RecordGroup, s::Symbol) = get_record(r, s)
getindex(r::RecordGroup, s::NTuple{N,Symbol}) where {N} = get_record(r, s)
getindex(r::RecordGroup, i) = get_record(r, i)

@doc raw"""
    RecordEvery <: RecordAction

record only every $i$th iteration.
Otherwise (optionally, but activated by default) just update internal tracking
values.

This method does not perform any record itself but relies on it's childrens methods
"""
mutable struct RecordEvery <: RecordAction
    record::RecordAction
    every::Int
    alwaysUpdate::Bool
    function RecordEvery(r::RecordAction, every::Int=1, alwaysUpdate::Bool=true)
        return new(r, every, alwaysUpdate)
    end
end
function (d::RecordEvery)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if i <= 0
        d.record(p, o, i)
    elseif (rem(i, d.every) == 0)
        d.record(p, o, i)
    elseif d.alwaysUpdate
        d.record(p, o, 0)
    end
end
get_record(r::RecordEvery) = get_record(r.record)
get_record(r::RecordEvery, i) = get_record(r.record, i)
getindex(r::RecordEvery, i) = get_record(r, i)

#
# Special single ones
#
@doc raw"""
    RecordChange <: RecordAction

debug for the amount of change of the iterate (stored in `o.x` of the [`Options`](@ref))
during the last iteration.

# Additional Fields
* `storage` a [`StoreOptionsAction`](@ref) to store (at least) `o.x` to use this
  as the last value (to compute the change)
"""
mutable struct RecordChange <: RecordAction
    recorded_values::Array{Float64,1}
    storage::StoreOptionsAction
    function RecordChange(a::StoreOptionsAction=StoreOptionsAction((:x,)))
        return new(Array{Float64,1}(), a)
    end
    function RecordChange(x0, a::StoreOptionsAction=StoreOptionsAction((:x,)))
        update_storage!(a, Dict(:x => x0))
        return new(Array{Float64,1}(), a)
    end
end
function (r::RecordChange)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    record_or_reset!(
        r,
        has_storage(r.storage, :x) ? distance(p.M, o.x, get_storage(r.storage, :x)) : 0.0,
        i,
    )
    r.storage(p, o, i)
    return r.recorded_values
end

@doc raw"""
    RecordEntry{T} <: RecordAction

record a certain fields entry of type {T} during the iterates

# Fields
* `recorded_values` – the recorded Iterates
* `field` – Symbol the entry can be accessed with within [`Options`](@ref)

"""
mutable struct RecordEntry{T} <: RecordAction
    recorded_values::Array{T,1}
    field::Symbol
    RecordEntry{T}(f::Symbol) where {T} = new(Array{T,1}(), f)
end
RecordEntry(::T, f::Symbol) where {T} = RecordEntry{T}(f)
RecordEntry(d::DataType, f::Symbol) = RecordEntry{d}(f)
function (r::RecordEntry{T})(::Problem, o::Options, i) where {T}
    return record_or_reset!(r, getfield(o, r.field), i)
end

@doc raw"""
    RecordEntryChange{T} <: RecordAction

record a certain entries change during iterates

# Additional Fields
* `recorded_values` – the recorded Iterates
* `field` – Symbol the field can be accessed with within [`Options`](@ref)
* `distance` – function (p,o,x1,x2) to compute the change/distance between two values of the entry
* `storage` – a [`StoreOptionsAction`](@ref) to store (at least) `getproperty(o, d.field)`
"""
mutable struct RecordEntryChange <: RecordAction
    recorded_values::Vector{Float64}
    field::Symbol
    distance::Any
    storage::StoreOptionsAction
    function RecordEntryChange(f::Symbol, d, a::StoreOptionsAction=StoreOptionsAction((f,)))
        return new(Float64[], f, d, a)
    end
    function RecordEntryChange(
        v::T where {T}, f::Symbol, d, a::StoreOptionsAction=StoreOptionsAction((f,))
    )
        update_storage!(a, Dict(f => v))
        return new(Float64[], f, d, a)
    end
end
function (r::RecordEntryChange)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    record_or_reset!(
        r,
        if has_storage(r.storage, r.field)
            r.distance(p, o, getfield(o, r.field), get_storage(r.storage, r.field))
        else
            0.0
        end,
        i,
    )
    return r.storage(p, o, i)
end

@doc raw"""
    RecordIterate <: RecordAction

record the iterate

# Constructors
    RecordIterate(x0)

initialize the iterate record array to the type of `x0`, e.g. your initial data.

    RecordIterate(P)

initialize the iterate record array to the data type `T`.
"""
mutable struct RecordIterate{T} <: RecordAction
    recorded_values::Array{T,1}
    RecordIterate{T}() where {T} = new(Array{T,1}())
end
RecordIterate(::T) where {T} = RecordIterate{T}()
function RecordIterate()
    return throw(
        ErrorException(
            "The iterate's data type has to be provided, i.e. RecordIterate(x0)."
        ),
    )
end

function (r::RecordIterate{T})(::Problem, o::Options, i) where {T}
    return record_or_reset!(r, o.x, i)
end

@doc raw"""
    RecordIteration <: RecordAction

record the current iteration
"""
mutable struct RecordIteration <: RecordAction
    recorded_values::Array{Int,1}
    RecordIteration() = new(Array{Int,1}())
end
function (r::RecordIteration)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    return record_or_reset!(r, i, i)
end

@doc raw"""
    RecordCost <: RecordAction

record the current cost function value, see [`get_cost`](@ref).
"""
mutable struct RecordCost <: RecordAction
    recorded_values::Array{Float64,1}
    RecordCost() = new(Array{Float64,1}())
end
function (r::RecordCost)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    return record_or_reset!(r, get_cost(p, o.x), i)
end

@doc raw"""
    RecordFactory(o::Options, a)

given an array of `Symbol`s and [`RecordAction`](@ref)s and `Ints`

* The symbol `:Cost` creates a [`RecordCost`](@ref)
* The symbol `:iteration` creates a [`RecordIteration`](@ref)
* The symbol `:Change` creates a [`RecordChange`](@ref)
* any other symbol creates a [`RecordEntry`](@ref) of the corresponding field in [`Options`](@ref)
* any [`RecordAction`](@ref) is directly included
* an semantic pair `:symbol => RecordAction` is directly included
* an Integer `k` introduces that record is only performed every `k`th iteration
"""
function RecordFactory(o::Options, a::Array{<:Any,1})
    # filter out every
    group = Array{Union{<:RecordAction,Pair{Symbol,<:RecordAction}},1}()
    for s in filter(x -> !isa(x, Int), a) # for all that are not integers or stopping criteria
        if s isa Symbol
            push!(group, s => RecordActionFactory(o, s))
        elseif s isa Pair{<:Symbol,<:RecordAction}
            push!(group, s)
        else
            push!(group, RecordActionFactory(o, s))
        end
    end
    record = RecordGroup(group)
    # filter ints
    e = filter(x -> isa(x, Int), a)
    if length(e) > 0
        record = RecordEvery(record, last(e))
    end
    return (; Iteration=record)
end
RecordFactory(o::Options, s::Symbol) = RecordActionFactory(o, s)

@doc raw"""
    RecordActionFactory(s)

create a [`RecordAction`](@ref) where

* a [`RecordAction`](@ref) is passed through
* a [`Symbol`] creates [`RecordEntry`](@ref) of that symbol, with the exceptions
  of `:Change`, `:Iterate`, `:Iteration`, and `:Cost`.
"""
RecordActionFactory(::Options, a::RecordAction) = a
RecordActionFactory(::Options, sa::Pair{Symbol,<:RecordAction}) = sa
function RecordActionFactory(o::Options, s::Symbol)
    if (s == :Change)
        return RecordChange()
    elseif (s == :Iteration)
        return RecordIteration()
    elseif (s == :Iterate)
        return RecordIterate(o.x)
    elseif (s == :Cost)
        return RecordCost()
    end
    return RecordEntry(getfield(o, s), s)
end


text:
#=
Triangle, pentagonal, and hexagonal numbers are generated by the following formulae:

Triangle	 	Tn=n(n+1)/2	 	1, 3, 6, 10, 15, ...
Pentagonal	 	Pn=n(3n−1)/2	 	1, 5, 12, 22, 35, ...
Hexagonal	 	Hn=n(2n−1)	 	1, 6, 15, 28, 45, ...
It can be verified that T285 = P165 = H143 = 40755.

Find the next triangle number that is also pentagonal and hexagonal.
=#
function ispentagonalnum(x)
  n = (sqrt(24x+1)+1)/6
  return n == round(Int, n) && n > 0
end

function ishexagonalnum(x)
  n = (sqrt(8x+1) + 1) / 4
  return n == round(Int,n) && n > 0
end

function calc()
  i = 286
  while true
    triangle = BigInt(i)*(i+1)/2
    if ispentagonalnum(triangle) && ishexagonalnum(triangle)
      return BigInt(triangle)
    end
    i += 1
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# Dense matrix linear operator

"""
  Dense array linear operator.  Serial only.

  **Public Fields**

   * A: an Array{Float64, 1}
"""
mutable struct DenseLO <: AbstractDenseLO
  A::Array{Float64, 2}
  ipiv::Array{BlasInt, 1}
  is_setup::Array{Bool, 1} # true if LO is already setup (ie. factored), false otherwise
  is_shared::Bool
  nfactorizations::Int
  nsolves::Int  # regular solves
  ntsolves::Int  # transpose solves

  # MPI stuff
  comm::MPI.Comm
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for DenseLO

  **Inputs**

   * pc: a PCNone
   * mesh
   * sbp
   * eqn
   * opts
"""
function DenseLO(pc::PCNone, mesh::AbstractMesh, sbp::AbstractOperator,
                 eqn::AbstractSolutionData, opts::Dict)

  @assert mesh.commsize == 1  # serial only
  A = zeros(mesh.numDof, mesh.numDof)
  ipiv = zeros(BlasInt, mesh.numDof)
  is_setup = Bool[false]
  is_shared = false
  nfactorizations = 0
  nsolves = 0
  ntsolves = 0
  comm = eqn.comm
  myrank = eqn.myrank
  commsize = eqn.commsize

  return DenseLO(A, ipiv, is_setup, is_shared, nfactorizations, nsolves,
                 ntsolves, comm, myrank, commsize)
end

function free(lo::DenseLO)
  # nothing to do here

  return nothing
end

function calcLinearOperator(lo::DenseLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

#  physicsJac(lo, mesh, sbp, eqn, opts, lo.A, ctx_residual, t)
  setIsSetup(lo, false)

  return nothing
end

# note: this may be faster if lo is not yet set up
function applyLinearOperator(lo::AbstractDenseLO, mesh::AbstractMesh,
                             sbp::AbstractOperator, eqn::AbstractSolutionData,
                             opts::Dict, ctx_residual, t, x::AbstractVector,
                             b::AbstractVector)

  lo2 = getBaseLO(lo)
  # BLAS operates in-place
  if getIsSetup(lo2)  # A is already factored
    throw(ErrorException("Multiplying factored matrix not supported"))
    # this doesn't work because the permutation is not correct for the vector
    # I think the problem is that ipiv is not a permuatation matrix, so really
    # what should happen is laswp(L)*U*b rather than laswp(L*U*b)
    #=
    println("multiplying in factored form")
    # A = P*L*U
    copy!(b, x)
    trmv!('U', 'N', 'N', lo2.A, b)
    trmv!('L', 'N', 'U', lo2.A, b)
    println("b = \n", b)
    println("ipiv = \n", lo2.ipiv)
    laswp!(b, 1, length(b), lo2.ipiv)
    println("b = \n", b)
    =#
  else
    gemv!('N', 1.0, lo2.A, x, 0.0, b)
  end

  return nothing
end

function applyLinearOperatorTranspose(lo::AbstractDenseLO,
                             mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::AbstractSolutionData, opts::Dict,
                             ctx_residual, t, x::AbstractVector,
                             b::AbstractVector)
  lo2 = getBaseLO(lo)
  if getIsSetup(lo2)  # A is already factored
    # A.'x = (P*L*U).'x = (U.'*L.'*P.')*x
    # there is no way to apply P.' in lapack (ipiv is not a true permutation)

    throw(ErrorException("applyLinearOperatorTranspose not supported for DenseLO in factored state: recompute the linear operator"))
  else
    gemv!('T', 1.0, lo2.A, x, 0.0, b)
  end

  return nothing
end

function getBaseLO(lo::DenseLO)

  # this is the bottom of the recursion tree
  return lo
end


text:
function get_lib()
    default_lib = "default"

    lib = @load_preference("shared_lib", default_lib)
    if lib != "default" && lib !== nothing
        if !isfile(lib)
            @error("User provided shared library: $lib is not a file.")
            lib = default_lib
        end
    else
        lib = default_lib
    end
    return lib
end

const artifact_or_path = get_lib()

"""
    set_lib!(path; export_prefs::Bool = false)

Set the shared library path for SuiteSparse:GraphBLAS. Set to "default" to use the provided artifact.
"""
function set_lib!(path; export_prefs::Bool = false)
    if path != "default" && path !== nothing
        if !isfile(path)
            @error("User provided shared library: $path is not a file.")
        end
    end
    set_preferences!(@__MODULE__, "shared_lib" => path; export_prefs, force = true)
    if path != artifact_or_path
        path = get_lib()
        @info("SuiteSparse:GraphBLAS shared library changed; restart Julia for this change to take effect", path)
    end
end


text:
#### thermodynamics

export new_thermo_state, recover_thermo_state

"""
    new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

Create a new thermodynamic state, based on the `state`, and _not_
the `aux` state.

!!! note
    This method calls the iterative saturation adjustment
    procedure for EquilMoist models.
"""
new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(atmos, atmos.moisture, state, aux)

function new_thermo_state(
    atmos::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseDry(atmos.param_set, e_int, state.ρ)
end

function new_thermo_state(
    atmos::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseEquil(
        atmos.param_set,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        moist.maxiter,
        moist.tolerance,
    )
end

function new_thermo_state(
    atmos::AtmosModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        q,
    )
end

"""
    recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

An atmospheric thermodynamic state.

!!! warn
    While recover_thermo_state is an ideal long-term solution,
    right now we are directly calling new_thermo_state to avoid
    inconsistent aux states in kernels where the aux states are
    out of sync with the boundary state.

# TODO:
    - Allow a safe way to call
    `recover_thermo_state(state, moist::EquilMoist, ...)`
"""
recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(atmos, atmos.moisture, state, aux)

function recover_thermo_state(
    atmos::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseDry(atmos.param_set, e_int, state.ρ)
end

function recover_thermo_state(
    atmos::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        aux.moisture.temperature,
    )
end

function recover_thermo_state(
    atmos::AtmosModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        q,
    )
end

max_forks_repo_name: jm-c/CLIMA


text:
module TuringPatterns

include("util.jl")
include("blur.jl")
include("export.jl")

# optionally load GeoStats.jl API
using Requires
function __init__()
    @require Meshes="eacbb407-ea5a-433e-ab97-5258b1ca43fa" begin
        @require GeoStatsBase="323cb8eb-fbf6-51c0-afd0-f8fba70507b2" include("geostats.jl")
    end
end

export Pattern, Params, SimplePattern, CompositePattern
export Sim, step!, simulate
export Zero, Clamp, BoxBlur, IteratedBoxBlur
export saveframe, saveframes, scale01

"An N-dimensional Turing pattern"
abstract type Pattern{N} end

"The difference between small and large scales in a pattern"
function difference end

"The effect of the pattern"
function effect end

"Allow a pattern to batch work; called once per frame"
function prepare! end

"Represents the parameters for a single Turing pattern"
struct Params
    r::Int
    R::Int
    effectsize::Float64
end

"Represents the simulation of a single Turing pattern"
struct SimplePattern{N} <: Pattern{N}
    params::Params
    blur_r::Array{Float64, N}
    blur_R::Array{Float64, N}
    difference::Array{Float64, N}
end

SimplePattern(params, sz) = SimplePattern(params, zeros(sz), zeros(sz), zeros(sz))

difference(p::SimplePattern, ind) = p.difference[ind]
effect(p::SimplePattern, ind) = p.params.effectsize * sign(difference(p, ind))
function prepare!(p::SimplePattern, blur::Blur, edge::EdgeCondition)
    # Compute small and large-scale averages and differences for all pixels
    blur!(blur, p.blur_r, p.params.r, edge)
    blur!(blur, p.blur_R, p.params.R, edge)
    @. p.difference = p.blur_r - p.blur_R
end

"Represents a pattern made up of multiple simple patterns"
struct CompositePattern{N} <: Pattern{N}
    patterns::Vector{SimplePattern{N}}
    couplings::Matrix{Float64}
end
CompositePattern(patterns) = CompositePattern(patterns, eye(length(patterns)))

difference(p::CompositePattern, ind) = sum(p.couplings * difference.(p.patterns, ind))
effect(p::CompositePattern, ind) = sum(p.couplings * effect.(p.patterns, ind))
prepare!(p::CompositePattern, blur, edge) =
    for p in p.patterns
        prepare!(p, blur, edge)
    end

# Promote heterogenous collections of `Pattern`s to a common type
Base.promote_rule(::Type{CompositePattern{N}}, ::Type{SimplePattern{N}}) where {N} = CompositePattern{N}
Base.convert(::Type{CompositePattern{N}}, p::SimplePattern{N}) where {N} = CompositePattern([p])

"Represents simulation state for multiple patterns acting on a 1D or 2D grid"
struct Sim{N, T<:Pattern}
    fluid::Array{Float64, N}
    patterns::Vector{T}
    edge::EdgeCondition
    blur::Blur
end

"Performs a single simulation step"
function step!(sim::Sim)
    # Prepare the blur and all of the patterns
    prepare!(sim.blur, sim.fluid)
    for p in sim.patterns
        prepare!(p, sim.blur, sim.edge)
    end

    # Select and apply the weakest pattern for every pixel
    for ind in eachindex(sim.fluid)
        p = let ind = ind # hack for better lambda optimization
            findmin(p -> abs(difference(p, ind)), sim.patterns)
        end
        sim.fluid[ind] += effect(p, ind)
    end

    # Re-normalize the fluid to [-1, 1]
    signed_scale01!(sim.fluid)
end

"Convenience method to simulate a set of patterns from some initial conditions"
function simulate(initial, patterns, blur, iters)
    sim = Sim(initial, patterns, Clamp(), blur)
    @time for i in 1:iters
        step!(sim)
    end
    scale01(sim.fluid)
end

end # module


text:
x = sort(rand(123))
ys = [sort(rand(456)), rand()]

@testset "mergesorted" for y in ys
    z, inds = mergesorted(x, y)

    @test z == sort(vcat(x, y))
    @test z[inds] == y
    @test z[setdiff(1:length(z), inds)] == x
end

n  = 1000
x  = range(0, 1, length = n)
atol = 1e-6

f(x) = x <= 0.5 ? 2x : 2x - 1

@testset "auc trapezoidal" begin
    @test auc_trapezoidal(x, x) ≈ 0.5  atol = atol
    @test auc_trapezoidal(x, x./2) ≈ 0.25 atol = atol
    @test auc_trapezoidal(x, f.(x)) ≈ 0.5  atol = atol
end


text:
using DelimitedFiles

# Converting V_0 and wfc into formats easily plotted by gnuplot
function convert_to_grid(xDim, filebase, start, incr, final)
    arr = zeros(xDim, div((final-start),incr)+1)
    count = 1
    for i = start:incr:final
        wfc = readdlm(filebase * "_" * string(i)) +
              readdlm(filebase * "i_" * string(i))*im
        arr[:,count] = abs2.(wfc)
        count += 1

    end
    writedlm("wfcgrid.dat", arr)
end

#convert_to_grid(512, "data/wfc_0_const", 0, 10, 5000)
convert_to_grid(512, "data/wfc_ev", 0, 10, 5000)

max_stars_repo_name: leios/thesis

max_issues_repo_name: leios/thesis

max_forks_repo_name: leios/thesis