text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -4.116064765345723 <= q <= -0.12537054293327043)

                     add_NL_constraint(m, :(log(1 + exp(-0.6620590692147021 + -0.9580027962881332*log(1 + exp(-0.38090481240177176 + -0.02939531016329644*log(1 + exp(0.2067810723024852 + 0.7469716803367321*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.3209180228281623*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + -0.9724733578803808*log(1 + exp(-0.4274287128450034 + -0.22403315015562253*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.6136110647248842*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + -0.699289524220764*log(1 + exp(0.4498619401955626 + 0.5463276179102481*log(1 + exp(0.2067810723024852 + 0.7469716803367321*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.3209180228281623*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + -0.02415027804893688*log(1 + exp(-0.4274287128450034 + -0.22403315015562253*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.6136110647248842*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))))) + log(1 + exp(-0.8107501452143855 + 0.008470632510253484*log(1 + exp(-0.38090481240177176 + -0.02939531016329644*log(1 + exp(0.2067810723024852 + 0.7469716803367321*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.3209180228281623*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + -0.9724733578803808*log(1 + exp(-0.4274287128450034 + -0.22403315015562253*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.6136110647248842*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + 0.39739497887500974*log(1 + exp(0.4498619401955626 + 0.5463276179102481*log(1 + exp(0.2067810723024852 + 0.7469716803367321*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.3209180228281623*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + -0.02415027804893688*log(1 + exp(-0.4274287128450034 + -0.22403315015562253*log(1 + exp(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3]))) + -0.6136110647248842*log(1 + exp(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@doc raw"""
    LogisticLikelihood() -> BernoulliLikelihood

Bernoulli likelihood with a logistic link for the Bernoulli likelihood
```math
    p(y|f) = \sigma(yf) = \frac{1}{1 + \exp(-yf)},
```
(for more info see : [wiki page](https://en.wikipedia.org/wiki/Logistic_function))

---

For the analytic version the likelihood, it is augmented via:
```math
    p(y|f,ω) = \exp\left(\frac{1}{2}(yf - (yf)^2 \omega)\right)
```
where ``ω \sim \mathcal{PG}(\omega | 1, 0)``, and ``\mathcal{PG}`` is the Polya-Gamma distribution.
See paper : [Efficient Gaussian Process Classification Using Polya-Gamma Data Augmentation](https://arxiv.org/abs/1802.06383).
"""
LogisticLikelihood() = BernoulliLikelihood(LogisticLink())

function implemented(
    ::BernoulliLikelihood{<:LogisticLink},
    ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling},
)
    return true
end

function Distributions.loglikelihood(
    ::BernoulliLikelihood{<:LogisticLink}, y::Real, f::Real
)
    return -log(one(f) + exp(-y * f))
end

function compute_proba(l::BernoulliLikelihood, f::Real)
    return pdf(l(f), 1)
end

### Local Updates Section ###
function local_updates!(
    local_vars,
    ::BernoulliLikelihood{<:LogisticLink},
    ::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
)
    map!(sqrt_expec_square, local_vars.c, μ, diagΣ) # √E[f^2]
    map!(local_vars.θ, local_vars.c) do c
        tanh(c / 2) / (2c)
    end # E[ω]
    return local_vars
end

function sample_local!(
    local_vars, ::BernoulliLikelihood{<:LogisticLink}, ::AbstractVector, f::AbstractVector
)
    map!(local_vars.θ, f) do f
        rand(PolyaGamma(1, abs(f)))
    end
    return local_vars
end

### Natural Gradient Section ###

function ∇E_μ(::BernoulliLikelihood{<:LogisticLink}, ::AOptimizer, y::AbstractVector, state)
    return (y / 2,)
end
function ∇E_Σ(::BernoulliLikelihood{<:LogisticLink}, ::AOptimizer, ::AbstractVector, state)
    return (state.θ / 2,)
end

### ELBO Section ###

function expec_loglikelihood(
    ::BernoulliLikelihood{<:LogisticLink},
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
    state,
)
    tot = -length(y) * logtwo / 2
    tot += (dot(μ, y) - dot(state.θ, diag_cov) - dot(state.θ, μ)) / 2
    return tot
end

AugmentedKL(l::BernoulliLikelihood{<:LogisticLink}, state, ::Any) = PolyaGammaKL(l, state)

function PolyaGammaKL(::BernoulliLikelihood{<:LogisticLink}, state)
    return sum(
        broadcast(PolyaGammaKL, ones(eltype(state.c), length(state.c)), state.c, state.θ)
    )
end

### Gradient Section ###

∇loglikehood(::BernoulliLikelihood{<:LogisticLink}, y::Real, f::Real) = y * logistic(-y * f)

function hessloglikehood(::BernoulliLikelihood{<:LogisticLink}, y::Real, f::Real)
    return -exp(y * f) / logistic(-y * f)^2
end

ext: jl


text:
############################

function compute_ordered_exp(bs::T, εₚ⁻::Float64, Δτ::Float64) where T<:Tuple{AbstractVector{Float64}, Vararg{AbstractVector{Float64}}}
	N = length(bs[1])
	M = one(SMatrix{2,2, Complex{Float128}})
	for i = N:-1:1
		M = M*Br(bs, i, εₚ⁻, Δτ)
	end
	return M
end

function compute_U_cash!(U_cash, bs::T, εₚ⁻, Δτ) where T<:Tuple{AbstractVector{Float64}, Vararg{AbstractVector{Float64}}}
	N = length(bs[1])
	M = one(SMatrix{2,2, Complex{Float128}})
	U_cash[N+1] = M
	for i′ = N:-1:1
		M = M*Br(bs, i′, εₚ⁻, Δτ)
		U_cash[i′] = M
	end
end

function compute_U_cash!(U_cash, bs::AbstractVector{Float64}, εₚ⁻, Δτ)
	N = length(bs)
	M = one(SMatrix{2,2, Complex{Float128}})
	U_cash[N+1] = M
	for i′ = N:-1:1
		M = M*Br(bs, i′, εₚ⁻, Δτ)
		U_cash[i′] = M
	end
end

############################

function compute_single_period!(∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, ∂²U::AbstractMatrix{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::NTuple{2, AbstractVector{Float64}}, εₚ⁻, Δτ)
	N = length(bs[1])
	∂₀U = @view ∂U[1:N]
	∂₁U = @view ∂U[(N+1):2N]
	∂²₀₀U = @view ∂²U[1:N, 1:N]
	∂²₁₁U = @view ∂²U[(N+1):2N, (N+1):2N]
	∂²₀₁U = @view ∂²U[1:N, (N+1):2N]
	∂²₁₀U = @view ∂²U[(N+1):2N, 1:N]
	compute_U_cash!(U_cash, bs, εₚ⁻, Δτ)
	Mb = one(SMatrix{2,2, Complex{Float128}})
	for i = 1:N
		∂²₀₀U[i,i] = U_cash[i+1]*B′′₀₀r(bs, i, εₚ⁻, Δτ)*Mb
		∂²₀₁U[i,i] = U_cash[i+1]*B′′₀₁r(bs, i, εₚ⁻, Δτ)*Mb
		∂²₁₀U[i,i] = ∂²₀₁U[i,i]
		∂²₁₁U[i,i] = U_cash[i+1]*B′′₁₁r(bs, i, εₚ⁻, Δτ)*Mb
		Mbm₀ = B′₀r(bs, i, εₚ⁻, Δτ)*Mb
		Mbm₁ = B′₁r(bs, i, εₚ⁻, Δτ)*Mb
		∂₀U[i] = U_cash[i+1]*Mbm₀
		∂₁U[i] = U_cash[i+1]*Mbm₁
		for j = i+1:N
			Me₀ = U_cash[j+1]*B′₀r(bs, j, εₚ⁻, Δτ)
			Me₁ = U_cash[j+1]*B′₁r(bs, j, εₚ⁻, Δτ)
			∂²₀₀U[j, i] = Me₀*Mbm₀
			∂²₀₀U[i, j] = ∂²₀₀U[j, i]
			∂²₁₁U[j, i] = Me₁*Mbm₁
			∂²₁₁U[i, j] = ∂²₁₁U[j, i]
			∂²₀₁U[j, i] = Me₀*Mbm₁
			∂²₀₁U[i, j] = Me₁*Mbm₀
			∂²₁₀U[i, j] = ∂²₀₁U[j, i]
			∂²₁₀U[j, i] = ∂²₀₁U[i, j]
			Bⱼr = Br(bs, j, εₚ⁻, Δτ)
			Mbm₀ = Bⱼr*Mbm₀
			Mbm₁ = Bⱼr*Mbm₁
		end
		Mb = Br(bs, i, εₚ⁻, Δτ)*Mb
	end
	return U_cash[1]
end

function compute_single_period!(∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, ∂²U::AbstractMatrix{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::Tuple{AbstractVector{Float64}}, εₚ⁻, Δτ)
	N = length(bs[1])
	compute_U_cash!(U_cash, bs, εₚ⁻, Δτ)
	Mb = one(SMatrix{2,2, Complex{Float128}})
	for i = 1:N
		∂²U[i,i] = U_cash[i+1]*B′′₀₀r(bs, i, εₚ⁻, Δτ)*Mb
		Mbm₀ = B′₀r(bs, i, εₚ⁻, Δτ)*Mb
		∂U[i] = U_cash[i+1]*Mbm₀
		for j = i+1:N
			∂²U[j, i] = U_cash[j+1]*B′₀r(bs, j, εₚ⁻, Δτ)*Mbm₀
			∂²U[i, j] = ∂²U[j, i]
			Mbm₀ = Br(bs, j, εₚ⁻, Δτ)*Mbm₀
		end
		Mb = Br(bs, i, εₚ⁻, Δτ)*Mb
	end
	return U_cash[1]
end

function compute_single_period!(∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::NTuple{3, AbstractVector{Float64}}, εₚ⁻, Δτ)
	N = length(bs[1])
	∂₀U = @view ∂U[1:N]
	∂₁U = @view ∂U[(N+1):2N]
	∂₃U = @view ∂U[(2N+1):3N]
	compute_U_cash!(U_cash, bs, εₚ⁻, Δτ)
	Mb = one(SMatrix{2,2, Complex{Float128}})
	for i = 1:N
		∂₀U[i] = U_cash[i+1]*B′₀r(bs, i, εₚ⁻, Δτ)*Mb
		∂₁U[i] = U_cash[i+1]*B′₁r(bs, i, εₚ⁻, Δτ)*Mb
		∂₃U[i] = U_cash[i+1]*B′₃r(bs, i, εₚ⁻, Δτ)*Mb
		Mb = Br(bs, i, εₚ⁻, Δτ)*Mb
	end
	return U_cash[1]
end

function compute_single_period!(∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::NTuple{2, AbstractVector{Float64}}, εₚ⁻, Δτ)
	N = length(bs[1])
	∂₀U = @view ∂U[1:N]
	∂₁U = @view ∂U[(N+1):2N]
	compute_U_cash!(U_cash, bs, εₚ⁻, Δτ)
	Mb = one(SMatrix{2,2, Complex{Float128}})
	for i = 1:N
		∂₀U[i] = U_cash[i+1]*B′₀r(bs, i, εₚ⁻, Δτ)*Mb
		∂₁U[i] = U_cash[i+1]*B′₁r(bs, i, εₚ⁻, Δτ)*Mb
		Mb = Br(bs, i, εₚ⁻, Δτ)*Mb
	end
	return U_cash[1]
end

function compute_single_period!(∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::Tuple{AbstractVector{Float64}}, εₚ⁻, Δτ)
	N = length(bs[1])
	compute_U_cash!(U_cash, bs, εₚ⁻, Δτ)
	Mb = one(SMatrix{2,2, Complex{Float128}})
	for i = 1:N
		∂U[i] = U_cash[i+1]*B′₀r(bs, i, εₚ⁻, Δτ)*Mb
		Mb = Br(bs, i, εₚ⁻, Δτ)*Mb
	end
	return U_cash[1]
end

############################

process_bs(bs, params::AbstractParamsNoB1) = (bs,)
process_bs(bs, params::AbstractParamsB1) = (view(bs, 1:params.N), view(bs,(params.N+1):2params.N))

function custom_eigen(U::SMatrix{2,2, Complex{T}}) where {T<:Real}
	eigvals, eigvecs = eigen(Matrix(U))
	if abs(eigvals[1])>=abs(eigvals[2])
		λ₁, λ₂ = real(eigvals[1]), real(eigvals[2])
		S = SMatrix{2,2}(eigvecs)
		S⁻¹ = inv(S)
	else
		λ₁, λ₂ = real(eigvals[2]), real(eigvals[1])
		S = SMatrix{2,2}([eigvecs[1,2] eigvecs[1,1]; eigvecs[2,2] eigvecs[2,1]])
		S⁻¹ = inv(S)
	end
	#print(λ₁, "\n", λ₂, "\n", S, "\n", S⁻¹, "\n\n")
	return λ₁, λ₂, S, S⁻¹
end

function compute_full_span!(∂Fₚ::AbstractVector{Float64}, ∂²Fₚ::AbstractMatrix{Float64}, ∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, ∂²U::AbstractMatrix{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::Vector{Float64}, params::AbstractParams, εₚ⁺, εₚ⁻)
	N₀, N, m, β = length(bs), params.N, params.m, params.β
	U = compute_single_period!(∂U, ∂²U, U_cash, process_bs(bs, params), εₚ⁻, β/(N*m))
	λ₁, λ₂, S, S⁻¹ = custom_eigen(U)
	mFₚ⁻¹ = m/(2*cosh(β*εₚ⁺)/λ₁^m + 1.0 + (λ₂/λ₁)^m)
	ll = 0.0
	for k=0:m-2
		ll = ll*(λ₂/λ₁) + 1.0
	end
	for i = 1:N₀
		∂Uⁱₛ = S⁻¹*∂U[i]*S
		∂Fₚ[i] = mFₚ⁻¹*(1/λ₁)*real(∂Uⁱₛ[1,1]+(λ₂/λ₁)^(m-1)*∂Uⁱₛ[2,2])
	end
	for j = 1:N₀
		∂Uʲₛ = S⁻¹*∂U[j]*S
		for i = 1:N₀
			∂²Uⁱʲₛ = S⁻¹*∂²U[i,j]*S
			∂Uⁱₛ = S⁻¹*∂U[i]*S
			Δf = (1/λ₁)*real(∂²Uⁱʲₛ[1,1] + (λ₂/λ₁)^(m-1)*∂²Uⁱʲₛ[2,2])
			if m>1
				Δf += (1/λ₁^2)*real((m-1)*(∂Uⁱₛ[1,1]*∂Uʲₛ[1,1] + (λ₂/λ₁)^(m-2)*∂Uⁱₛ[2,2]*∂Uʲₛ[2,2]) + ll*(∂Uⁱₛ[1,2]*∂Uʲₛ[2,1] + ∂Uⁱₛ[2,1]*∂Uʲₛ[1,2]))
			end
			∂²Fₚ[i,j] = mFₚ⁻¹*Δf-∂Fₚ[i]*∂Fₚ[j]
		end
	end
	β′, λ₁′, λ₂′ = Float128(β), Float128(λ₁), Float128(λ₂)
	return Float64(log((cosh(β′*εₚ⁺) + 0.5*(λ₁′^m + λ₂′^m))/(cosh(β′*εₚ⁺)+cosh(β′*εₚ⁻))))
end

function compute_full_span!(∂Fₚ::AbstractVector{Float64}, ∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::AbstractVector{Float64}, params::AbstractParams{Int64}, εₚ⁺, εₚ⁻, wₚ, a, b)
	N₀, N, m = length(bs), params.N, params.m
	U = compute_single_period!(∂U, U_cash, process_bs(bs, params), εₚ⁻, params.Δτ)
	λ₁, λ₂, S, S⁻¹ = custom_eigen(U)
    λ₁ᵐ⁻¹ = λ₁^(m-1)
    λ₂ᵐ⁻¹ = λ₂^(m-1)
    Fₚ = a + λ₁^m + λ₂^m
	Fₚ⁻¹ = 1/Fₚ
	for i = 1:N₀
		∂Uⁱₛ = S⁻¹*∂U[i]*S
        ∂Fₚ[i] += wₚ*Fₚ⁻¹*real(λ₁ᵐ⁻¹*∂Uⁱₛ[1,1]+λ₂ᵐ⁻¹*∂Uⁱₛ[2,2])
	end
    return Float64(wₚ*(log(abs(Fₚ))-b)/m)
end

function compute_full_span!(∂Fₚ::AbstractVector{Float64}, ∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::AbstractVector{Float64}, params::AbstractParams{Nothing}, εₚ⁺, εₚ⁻, wₚ, a, b)
	N₀, N, m = length(bs), params.N, params.m
	U = compute_single_period!(∂U, U_cash, process_bs(bs, params), εₚ⁻, params.Δτ)
	λ₁, λ₂, S, S⁻¹ = custom_eigen(U)
	Fₚ⁰ = log(abs(λ₁))
	if Fₚ⁰>=a
		Fₚ⁻¹ = 1/λ₁
		for i = 1:N₀
			∂Uⁱₛ = S⁻¹*∂U[i]*S
			∂Fₚ[i] += wₚ*Fₚ⁻¹*real(∂Uⁱₛ[1,1])
		end
		return Float64(wₚ*(Fₚ⁰ - b))
	else
		return Float64(wₚ*(a - b))
	end
end

function compute_full_span(bs::Vector{Float64}, params::AbstractParams{Int64}, εₚ⁺, εₚ⁻, wₚ, a, b)
	m, N = params.m, params.N
	U = compute_ordered_exp(process_bs(bs, params), εₚ⁻, params.Δτ)
	λ₁, λ₂, S, S⁻¹ = custom_eigen(U)
    return Float64(wₚ*(log(abs(a + (λ₁^m + λ₂^m)))-b)/m)
end

function compute_full_span(bs::Vector{Float64}, params::AbstractParams{Nothing}, εₚ⁺, εₚ⁻, wₚ, a, b)
	m, N = params.m, params.N
	U = compute_ordered_exp(process_bs(bs, params), εₚ⁻, params.Δτ)
	λ₁, λ₂, S, S⁻¹ = custom_eigen(U)
	return Float64(wₚ*(max(log(abs(λ₁)), a) - b))
end
############################

function process_chunk!(∂F::AbstractVector{Float64}, ∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs, params::AbstractParams, psamples::AbstractVector{Tuple{Float64, Float64, Float64, Float128, Float128}})
	fill!(∂F, 0.0)
	ΔF = zero(Float64)
	for psample in psamples
		ΔF += compute_full_span!(∂F, ∂U, U_cash, bs, params, psample...)
	end
	return ΔF
end

function process_chunk2!(ΔFs::Vector{Float64}, ∂Fs::Vector{Vector{Float64}}, ∂U::AbstractVector{SMatrix{2,2, Complex{Float128}}}, U_cash::AbstractVector{SMatrix{2,2, Complex{Float128}}}, bs::AbstractVector{Float64}, params::AbstractParams, psamples::AbstractVector{Tuple{Float64, Float64, Float64, Float128, Float128}})
    for i in eachindex(psamples)
        fill!(∂Fs[i], 0.0)
        ΔFs[i] = compute_full_span!(∂Fs[i], ∂U, U_cash, bs, params, psamples[i]...)
	end
end

function process_chunk(bs::AbstractVector{Float64}, params::AbstractParams, psamples::AbstractVector{Tuple{Float64, Float64, Float64, Float128, Float128}})
	ΔF = zero(Float64)
	for psample in psamples
		ΔF += compute_full_span(bs, params, psample...)
	end
	return ΔF
end

true_workers() = nworkers()>1 ? workers() : Int64[]

function precompute_step!(G_storage::G_Cash, bs, params::AbstractParams, psamples)
	∂Fs, ∂Us, U_cashes = G_storage.∂Fs, G_storage.∂Us, G_storage.U_cashes
	fs = Vector{Future}(undef, nprocs())
    @sync begin
        for wid in true_workers()
            fs[wid] = @spawnat wid process_chunk!(∂Fs[wid], localpart(∂Us), localpart(U_cashes), bs, params, localpart(psamples))
        end
        fs[1] = @spawnat 1 process_chunk!(∂Fs[1], localpart(∂Us), localpart(U_cashes), bs, params, localpart(psamples))
    end
	ΔF = sum(fetch.(fs))
	#print(typeof(∂Fs[2]), " ", length(∂Fs[2]), "\n")
    for wid in true_workers()
		∂Fs[1] .+= ∂Fs[wid]
	end
	return finalize!(∂Fs[1], ΔF, bs, params)
end

function compute_grad_components!(G_storage::G_Cash2, bs, params::AbstractParams, psamples)
	ΔFss, ∂Fss, ∂Us, U_cashes = G_storage.ΔFss, G_storage.∂Fss, G_storage.∂Us, G_storage.U_cashes
    @sync begin
        for wid in true_workers()
            @spawnat wid process_chunk2!(localpart(ΔFss), localpart(∂Fss), localpart(∂Us), localpart(U_cashes), bs, params, localpart(psamples))
        end
        @spawnat 1 process_chunk2!(localpart(ΔFss), localpart(∂Fss), localpart(∂Us), localpart(U_cashes), bs, params, localpart(psamples))
    end
	return G_storage
end

function precompute_step(bs::AbstractVector{Float64}, params::AbstractParams, psamples)
	fs = Vector{Future}(undef, nprocs())
    @sync begin
        for wid in true_workers()
			fs[wid] = @spawnat wid process_chunk(bs, params, localpart(psamples))
        end
		fs[1] = @spawnat 1 process_chunk(bs, params, localpart(psamples))
    end
	ΔF = sum(fetch.(fs))
   	return finalize(ΔF, bs, params)
end

function finalize!(∂F, ΔF, bs, params::AbstractParamsB1)
	N, Δτ, u, u₁ = params.N, params.Δτ, params.u, params.u₁
    ΔF = Δτ*(u*sum(x->(abs2(x)-1), bs[1:N]) + u₁*sum(abs2, bs[(N+1):2N])) - ΔF
	∂F .*= -1.0
	for i = 1:N
		∂F[i] += 2Δτ*u*bs[i]
	end
	for i = (N+1):2N
		∂F[i] += 2Δτ*u₁*bs[i]
	end
    return ΔF, pin_bs!(∂F, params)
end

function finalize(ΔF, bs, params::AbstractParamsB1)
	N, Δτ, u, u₁ = params.N, params.Δτ, params.u, params.u₁
	return Δτ*(u*sum(x->(abs2(x)-1), bs[1:N]) + u₁*sum(abs2, bs[(N+1):2N])) - ΔF
end

function finalize!(∂F, ΔF, bs, params::AbstractParamsNoB1)
	N, Δτ, u = params.N, params.Δτ, params.u
	ΔF = Δτ*u*sum(x->(abs2(x)-1), bs) - ΔF
	∂F .*= -1.0
	for i = 1:N
		∂F[i] += 2Δτ*u*bs[i]
	end
    return ΔF, pin_bs!(∂F, params)
end

function finalize(ΔF, bs, params::AbstractParamsNoB1)
	N, Δτ, u = params.N, params.Δτ, params.u
	return Δτ*u*sum(x->(abs2(x)-1), bs) - ΔF
end


text:
using Test

@testset "$(basename(@__FILE__))" begin
    @testset "1. using JustSayIt" begin
        using JustSayIt
        @test true
    end;
    @testset "2. using JustSayIt.API" begin
        using JustSayIt.API
        @test true
    end;
end;


text:
using JuliaFEM
using Lexicon

include("badges.jl")

"""
Searches recursively all the modules from packages. As documentation grows, it's a bit
troublesome to add all the new modules manually, so this function searches all the modules
automatically.

Parameters
----------
    module_: Module
        Module where we want to search modules inside
    append_list: Array{Module, 1}
        Array, where we append Modules as we find them

Returns
-------
    None. Void function, which manipulates the append_list
"""
function search_modules!(module_::Module, append_list::Array{Module, 1})
    all_names = names(module_, true)
    for each in all_names
        inner_module = module_.(each)
        if (typeof(inner_module) == Module) && !(inner_module in append_list)
            push!(append_list, inner_module)
            search_modules!(inner_module, append_list)
        end
    end
end

append_list = Array(Module, 0)
search_modules!(JuliaFEM, append_list)

const modules = append_list

cd(dirname(@__FILE__)) do
    npassed = 0
    nfailed = 0
    nskipped = 0

    for m in modules
        s = doctest(m)
        lnpassed, lnfailed, lnskipped = map(length, (passed(s), failed(s), skipped(s)))
        npassed += lnpassed
        nfailed += lnfailed
        nskipped += lnskipped
    end
    println("""DOCTEST: {"failed": $nfailed, "skipped": $nskipped, "passed": $npassed}""")
    make_badge("doctests", npassed, (npassed+nfailed+nskipped), "badges/doctests-status.svg")
end


text:
const cp = 1004.0
const R = 287.0
const p0 = 1000.0

"""
    temperature(theta, exner; celsius=false)
Returns temperature given potential temperature (`THETA` in RAMS)
and the exner function * c_p (`PI` in RAMS). Specify `celsius=true` for celsius.
"""
function temperature(theta, exner; celsius=false)
    T = theta .* exner ./ cp

    if celsius
        T = T .- 273.15
    end

    return T
end

"""
    pressure(exner)
Return the pressure given the Exner function (`PI` in RAMS).
"""
function pressure(exner)
    p = p0 .* (exner ./ cp).^(cp/R)
    return p
end

"""
    wsat(theta, exner)
Saturation mixing ratio based on `THETA` and exner (`PI`).

Note: for this function `THETA` and `PI` must be arrays
"""
function wsat(theta, exner)
    T = temperature(theta, exner; celsius=true)
    p = pressure(exner)

    c0 = 0.6105851e3
    c1 = 0.4440316e2
    c2 = 0.1430341e1
    c3 = 0.2641412e-1
    c4 = 0.2995057e-3
    c5 = 0.2031998e-5
    c6 = 0.6936113e-8
    c7 = 0.2564861e-11
    c8 = -0.3704404e-13

    T[T .< -80] .= -80
    x = T
    es = c0.+x.*(c1.+x.*(c2.+x.*(c3.+x.*(c4.+x.*(c5.+x.*(c6.+x.*(c7.+x.*c8)))))))

    ws = 0.622 .* es ./ (p.*100 .- es)
    return ws
end

"""
    rh(rv, theta, exner)
Calculates RH in percent (divide by 100 for fraction)
Note: for this function `THETA` and `PI` must be arrays
"""
function rh(rv, theta, exner)
    ws = wsat(theta, exner)
    rh = (rv ./ ws) .* 100.0

    return rh

end

"""
    mslp(temp, press, height)
Calculate Mean Sea Level Pressure from temperature, pressure, and height
"""
function mslp(temp, press, height)
    slp = press .* (1 .- (0.0065 .* height)./(temp .+ 0.0065 .* height .+ 273.15)).^-5.257
    return slp
end


text:
# Autogenerated wrapper script for GTK4_jll for aarch64-apple-darwin
export libgtk4

using Glib_jll
using Graphene_jll
using Cairo_jll
using Pango_jll
using FriBidi_jll
using FreeType2_jll
using gdk_pixbuf_jll
using Libepoxy_jll
using HarfBuzz_jll
using iso_codes_jll
using Fontconfig_jll
JLLWrappers.@generate_wrapper_header("GTK4")
JLLWrappers.@declare_library_product(libgtk4, "@rpath/libgtk-4.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, Graphene_jll, Cairo_jll, Pango_jll, FriBidi_jll, FreeType2_jll, gdk_pixbuf_jll, Libepoxy_jll, HarfBuzz_jll, iso_codes_jll, Fontconfig_jll)
    JLLWrappers.@init_library_product(
        libgtk4,
        "lib/libgtk-4.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module FEFMM

    import Base: isless
    using DataStructures

    export fefmm

    #small helper functions
    include("fefmm_nodes.jl")
    include("fefmm_neighbour_checks.jl")
    include("fefmm_piecewise_quadratic.jl")
    include("fefmm_utilities.jl")
    #core functions (initialization, main loop, node solve)
    include("factoredeikonalfmm.jl")

end


text:
"""
Graph represanted by an explicit vector of edges
"""
struct EdgeGraph <: Graph
    n::Int
    edges::Vector{Edge{Int}}
end

@inline num_nodes(g::EdgeGraph) = g.n

@inline num_edges(g::EdgeGraph) = length(g.edges)

@inline function enumerate_edges(f::F, g::EdgeGraph) where {F<:Function}
    for (i, (u, v)) in enumerate(g.edges)
        f(i, u, v)
    end
end

function Base.collect(g::Graph)::EdgeGraph
    m = num_edges(g)
    n = num_nodes(g)
    edges = Vector{Edge{Int}}(undef, m)
    enumerate_edges(g) do i::Int, u::Int, v::Int
        edges[i] = Edge((u, v))
    end
    EdgeGraph(n, edges)
end

Base.collect(g::EdgeGraph)::EdgeGraph = g

EdgeGraph(g::Graph) = collect(g)

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# test file for LittleEndianBase128.jl
using Test

include("../src/LittleEndianBase128.jl")
using .LittleEndianBase128

println("Testing LittleEndianBase128 encoding and decoding.")

print("Known results ... ")
@test encode(UInt32(624485)) ==  UInt8[0xe5,0x8e,0x26]
@test encode(UInt64(2147483647)) == UInt8[0xff,0xff,0xff,0xff,0x07]
@test encode(-1) == UInt8[0x01]
@test encode(0) == UInt8[0x00]
@test encode(2147483647) == UInt8[0xfe,0xff,0xff,0xff,0x0f]
@test encode(-2147483647) == UInt8[0xfd,0xff,0xff,0xff,0x0f]
@test decodesigned(encode(-2147483647)) == [-2147483647]
@test decode(UInt8(0x01)) == [0x01]
println("PASS")

println("Type min and max ...")
types = [UInt8, UInt16, UInt32, UInt64, UInt64, Int8, Int16, Int32, Int64, Int64, Bool]
n = 3
for t in types
  a = typemin(t)
  b = typemax(t)
  print("$t min ... ")
  u = decode(encode(a),t)[1]
  @test u == a
  @test typeof(u) == typeof(a)
  println("PASS")
  print("$t max ... ")
  u = decode(encode(b),t)[1]
  @test u == b
  @test typeof(u) == typeof(b)
  println("PASS")
  print("$t random $n x $n matrix ...")
  x = rand(a:b, n, n)
  y = reshape(decode(encode(x),t), n, n)
  @test x == y
  println("PASS")
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    vecmag(vec::Matrix{F})::Matrix{F} where {F<:AbstractFloat}

or

    vecmag(vec::Vector{F})::F where {F<:AbstractFloat}

Determines the magnitude of vectors.

# Arguments:
- `vec`: An array of column vectors represented by a matrix of size `(dim,nv)`,
  where `dim` is the dimension of the vector and `nv` the number of vectors.

# Output:
- `mvec`: Magnitude of the vectors, vector of size `(1,nv)`.

# Examples:
```julia
julia> mvec = vecmag(vec)
```

Find the magnitude of the two vectors `[0.0;2.0;1.3]` and `[1.5;3.4;2.3]`.
```julia
julia> mvec = vecmag([0.0 1.5; 2.0 3.4; 1.3 2.3])
```
"""
function vecmag(vec::Matrix{F})::Matrix{F} where {F<:AbstractFloat}
return sqrt.(sum(vec.^2,dims=1));
end # vecmag

# TODO: look if this is even necessary for this toolbox (using LinearAlgebra)
# overloaded function for vector input to get a float instead of a vector
function vecmag(vec::Vector{F})::F where {F<:AbstractFloat}
return sqrt(vec' * vec);
end # vecmag


text:
include((@__DIR__)*"/../preamble.jl")
include("reflecting_model/model_def.jl")

include("default_params.jl")

# look at time to exit of the interval [0,1]
# use the bounded fluid queue, but change the boundary to be at 1

dg_model(order) = build_discretised_model(DGMesh{t},model,1.0,order)
order1_model(order) = build_discretised_model(DGMesh{t},model,1.0/order,1)
qbdrap_model(order) = build_discretised_model(FRAPMesh{t},model,1.0,order)

d0_map_point_mass(dq) = (interior_point_mass(eps(),1,dq).coeffs)

# change d0_map here to do the make_approximations for ohter initial confitions
# remember to change write directory below too
models = make_approximations(orders,approx_types,d0_map_point_mass,(args...)->args[1])

h = 0.005
pth = mkpath((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2")
pth_coeff = mkpath((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2/approx_coeffs")
pth_cdf = mkpath((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2/approx_cdf")
x_vals = range(0,10;length=n_err_evals)
for k1 in keys(models)
    for k2 in keys(models[k1])
        models[k1][k2]["coeffs_mapped"] = Vector{Float64}(JSON.parsefile(
            (@__DIR__)*
            "/reflecting_model/transient_distribution/point_mass/data/models/order_"*string(k1)*
            "/model_"*k2*"/coeffs_mapped.json",
        ))
        @show k1
        @show k2
        if k2=="dg2"
            models[k1][k2]["coeffs_mapped"] = integrate_time(SFMDistribution(models[k1][k2]["coeffs_mapped"],models[k1][k2]["dq"]),models[k1][k2]["B"],0.1,StableRK4(h); limiter=GeneralisedMUSCL)
        elseif k2=="dg1"
            models[k1][k2]["coeffs_mapped"] = integrate_time(SFMDistribution(models[k1][k2]["coeffs_mapped"],models[k1][k2]["dq"]),models[k1][k2]["B"],0.1,StableRK4(h); limiter=NoLimiter)
        else
            models[k1][k2]["coeffs_mapped"] = integrate_time(SFMDistribution(models[k1][k2]["coeffs_mapped"],models[k1][k2]["dq"]),models[k1][k2]["B"],0.1,StableRK4(h))
        end
        CSV.write(pth_coeff*"/order_"*string(k1)*"model_"*k2*".csv",DataFrame(coeffs=models[k1][k2]["coeffs_mapped"]))
        cdf_evaluated = DataFrame(
            x=x_vals,
            phase_1 = cdf(models[k1][k2]["coeffs_mapped"]).(x,1),
            phase_2 = cdf(models[k1][k2]["coeffs_mapped"]).(x,2),
        )
        CSV.write(pth_cdf*"/order_"*string(k1)*"model_"*k2*".csv",cdf_evaluated)
    end
end


text:
module ReinforcementLearningZoo

const RLZoo = ReinforcementLearningZoo
export RLZoo

using ReinforcementLearningBase
using ReinforcementLearningCore
using Setfield: @set

include("patch.jl")
include("algorithms/algorithms.jl")
include("utils.jl")

using Requires

# dynamic loading environments
function __init__()
    @require ReinforcementLearningEnvironments = "25e41dd2-4622-11e9-1641-f1adca772921" begin
        include("experiments/rl_envs.jl")
        @require ArcadeLearningEnvironment = "b7f77d8d-088d-5e02-8ac0-89aab2acc977" include("experiments/atari.jl")
    end
end

end # module


text:
module GlobalSensitivityAnalysis

using DataStructures
using Distributions
using Statistics
using Sobol

export
    sample,
    analyze,
    SobolData,
    DeltaData

include("utils.jl")
include("sample_sobol.jl")
include("sample_delta.jl")
include("analyze_sobol.jl")
include("analyze_delta.jl")
include("test_functions/ishigami.jl")
include("test_functions/borehole.jl")

end # module

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
function scale(he::AbstractMeanfieldEquations; kwargs...)
    h = hilbert(he.hamiltonian)
    scale_aons, N = get_cluster_stuff(h)
    return scale(he, scale_aons, N; kwargs...)
end
function scale(he::MeanfieldEquations, scale_aons::Vector{<:Vector}, N::Vector;
                                                    simplify=true,
                                                    order=nothing, mix_choice=maximum,
                                                    kwargs...)
    names = get_names(he)
    equations = Symbolics.Equation[]
    for i=1:length(he.equations)
        lhs = he.operator_equations[i].lhs
        rhs = he.operator_equations[i].rhs
        push!(equations, average(lhs) ~ average(rhs))
    end

    for j=1:length(scale_aons)
        scale_aons_ = scale_aons[j]
        M = length(scale_aons_)
        for i=1:length(he.equations)
            rhs = _scale(he.states[i], equations[i].rhs, scale_aons_, N[j], M, names)
            if order===nothing
                rhs = substitute_redundants(rhs, scale_aons_, names)
            end
            equations[i] = Symbolics.Equation(equations[i].lhs, rhs)
        end
    end

    operator_equations = he.operator_equations
    for i=1:length(operator_equations)
        rhs = undo_average(equations[i].rhs)
        operator_equations[i] = Symbolics.Equation(operator_equations[i].lhs, rhs)
    end

    if order !== nothing
        for i=1:length(equations)
            rhs = cumulant_expansion(equations[i].rhs, order; simplify=simplify, mix_choice=mix_choice)
            for j=1:length(scale_aons)
                rhs = substitute_redundants(rhs, scale_aons[j], names)
            end
            equations[i] = equations[i].lhs ~ rhs
        end
    end

    he_scaled = ScaledMeanfieldEquations(equations,operator_equations,
                                    he.states, he.operators,
                                    he.hamiltonian,he.jumps,he.jumps_dagger,he.rates,
                                    he.iv, he.varmap, he.order,
                                    scale_aons,names,ones(Bool, length(he)
                                    )
    )

    if simplify
        return Symbolics.simplify(he_scaled; kwargs...)
    else
        return he_scaled
    end
end

function get_cluster_stuff(h::ClusterSpace,aon=1)
    M = h.order
    aons = [ClusterAon(aon,i) for i=1:M]
    return [aons], [h.N]
end
function get_cluster_stuff(h::ProductSpace)
    idx = findall(x->isa(x,ClusterSpace),h.spaces)
    aons = Vector{<:ClusterAon}[]
    N = []
    for i∈idx
        aon_, N_ = get_cluster_stuff(h.spaces[i],i)
        append!(aons, aon_)
        append!(N, N_)
    end
    return aons, N
end

function _scale(lhs, rhs, scale_aons, N, M, names)
    aon_lhs = acts_on(lhs)
    scale_aon_ = intersect(aon_lhs, scale_aons)
    M_ = length(scale_aon_)
    if iszero(M_)
        N_ = N / M
    elseif M_ == M # all clusters on lhs
        N_ = 1
    else
        N_ = (N - (M-M_)) / (M - M_)
    end
    rhs = _scaling_rewrite(rhs, N_, scale_aons, aon_lhs, names)
    return rhs
end

function _get_names_ops(ops)
    hs_ = hilbert(ops[1])
    if isa(hs_, ProductSpace)
        hs = hilbert(ops[1]).spaces
    else
        hs = [hs_]
    end
    names = Vector{Any}(undef, length(hs))
    for op in ops
        aon = acts_on(op)
        if isa(aon, Int)
            names[aon] = op.name
        else #ClusterAon
            h = hs[aon.i]
            order = h.order
            op_name = h.op_name[]
            cluster_names = [Symbol(op_name, :_, i) for i=1:order]
            names[aon.i] = cluster_names
        end
    end
    return names
end

function get_names(q::Union{QSym,QMul})
    ops = get_operators(q)
    unique_ops!(ops)
    _get_names_ops(ops)
end
function get_names(he)
    H = he.hamiltonian
    J = he.jumps
    ops = get_operators(H)
    for j ∈ J
        append!(ops, get_operators(j))
    end
    for l ∈ he.operators
        append!(ops, get_operators(l))
    end
    unique_ops!(ops)
    _get_names_ops(ops)
end

function unique_i_aons(aon)
    seen = eltype(aon)[]
    for a∈aon
        if !(a∈seen)
            if a isa Integer
                push!(seen, a)
            else # ClusterAon
                isone(a.j) && push!(seen, a)
            end
        end
    end
    return seen
end

get_operators(q::QSym) = [q]
function get_operators(q::QMul)
    ops = QSym[]
    seen_hashes = UInt[]
    for arg ∈ q.args_nc
        h = hash(arg)
        if !(h ∈ seen_hashes)
            push!(seen_hashes, h)
            push!(ops, arg)
        end
    end
    return ops
end
function get_operators(q::QAdd)
    ops = QSym[]
    for arg∈q.arguments
        append!(ops, get_operators(arg))
    end
    unique_ops!(ops)
    return ops
end

## Scaling terms by the correct factor

function _scaling_rewrite(rhs, N, scale_aons, aon_lhs, names)
    rw = let was_scaled=UInt[]

        SCALE_RULES = [
            SymbolicUtils.@rule(~x::should_scale(scale_aons,aon_lhs,was_scaled) => N*~x)
        ]

        rule_tree = [
            SymbolicUtils.If(
                SymbolicUtils.is_operation(sym_average), SymbolicUtils.Chain(SCALE_RULES)
                ),
        ] |> SymbolicUtils.Chain
        SymbolicUtils.Postwalk(rule_tree)
    end

    f = SymbolicUtils.Fixpoint(rw)
    return f(rhs)
end

should_scale(scale_aons,aon_lhs,was_scaled) = (x->should_scale(x, scale_aons, aon_lhs, was_scaled))
function should_scale(x, scale_aons, aon_lhs, was_scaled)
    h = hash(x)
    h ∈ was_scaled && return false # x was already scaled
    aon = [acts_on(x)...]
    intersect!(aon, scale_aons)
    all(a ∈ aon_lhs for a ∈ aon) && return false # x acts only on things contained in lhs
    should_scale_ = any(a ∈ scale_aons for a ∈ aon) # no scaling should occur if x is not part of a cluster
    should_scale_ && push!(was_scaled, h)
    return should_scale_
end

## Dealing with redundant averages

function substitute_redundants(he::MeanfieldEquations,scale_aons::Vector{<:Vector},names)
    eqs = Symbolics.Equation[]
    for i=1:length(he.equations)
        lhs = substitute_redundants(he.equations[i].lhs,scale_aons,names)
        rhs = substitute_redundants(he.equations[i].rhs,scale_aons,names)
        push!(eqs, Symbolics.Equation(lhs, rhs))
    end
    # TODO substitute jumps and hamiltonian?
    return MeanfieldEquations(lhs, rhs, he.hamiltonian, he.jumps, he.jumps_dagger, he.rates)
end

function substitute_redundants(t::SymbolicUtils.Symbolic,scale_aons::Vector{<:Vector},names)
    if SymbolicUtils.istree(t)
        f = SymbolicUtils.operation(t)
        if f === sym_average
            op = deepcopy(SymbolicUtils.arguments(t)[1])
            for j=1:length(scale_aons)
                op = substitute_redundants(op,scale_aons[j],names)
            end
            return average(op)
        else
            args = []
            for arg in SymbolicUtils.arguments(t)
                push!(args, substitute_redundants(arg,scale_aons,names))
            end
            return SymbolicUtils.similarterm(t, f, args)
        end
    else
        return t
    end
end
substitute_redundants(x::Number,args...) = x
substitute_redundants(x,scale_aons,names) = substitute_redundants(x,[scale_aons],names)

function substitute_redundants(t::QMul,scale_aons,names)
    aon = acts_on(t)
    idx_aon = findall(in(scale_aons), aon)
    isempty(idx_aon) && return t
    if is_redundant_aon(t,scale_aons)
        if length(aon)==1
            aon_subs = scale_aons[idx_aon[1]]
        else
            aon_subs = copy(aon)
            aon_subs[idx_aon] .= scale_aons[1:length(idx_aon)]
        end
        name_idx = map(get_i, aon_subs)
        return substitute_redundants(
                            _swap_aon_and_name(t, aon, aon_subs, names[name_idx]),
                            scale_aons, names
                            )
    else
        args = t.args_nc
        args_cluster = filter(x->acts_on(x)∈scale_aons, args)

        # Get the proper ordering
        p = sortperm_ref_order(args_cluster)

        if issorted(p) # arguments are already in correct order
            return t
        else
            # acts_on and relevant names
            aon_subs = copy(aon)
            names_ = names[map(get_i, aon)]

            # Permute cluster part to proper reference order
            aon_subs[idx_aon[p]] = scale_aons[1:length(idx_aon)]
            names_[idx_aon] = names_[idx_aon][p]

            # Swap and return
            return _swap_aon_and_name(t, aon, aon_subs, names_)
        end
    end
end
function substitute_redundants(x::QSym,scale_aons,names)
    if is_redundant_aon(x,scale_aons)
        aon = acts_on(x)
        if aon isa ClusterAon
            i = aon.i
            idx_aon = findfirst(x->x.i==i, scale_aons)
        else
            idx_aon = findfirst(in(scale_aons), aon)
        end
        aon_sub = scale_aons[idx_aon]
        return _swap_aon_and_name(x,aon,aon_sub,names)
    else
        return x
    end
end

function is_redundant_aon(x,scale_aons)
    # Judge whether a term is redundant from its acts_on
    aon = acts_on(x)
    if aon isa ClusterAon
        aon ∈ scale_aons || return false
        return aon != scale_aons[1]
    else
        idx_aon = findall(in(aon), scale_aons)
        isempty(idx_aon) && return false
        for i=1:length(idx_aon)-1
            (idx_aon[i]+1 == idx_aon[i+1]) || return true
        end
        idx = findfirst(in(scale_aons), aon)
        return aon[idx] != scale_aons[1]
    end
end

function sortperm_ref_order(args_cluster)
    if args_cluster[1] isa Transition
        return sortperm(args_cluster, lt=lt_reference_order)
    else
        # non-unique acts_on
        arg_aon = map(acts_on, args_cluster)

        # Separate arguments into blocks with equal acts_on
        i = 1
        blocks = []
        while i <= length(args_cluster)
            idx = findall(isequal(arg_aon[i]),arg_aon[i:end])
            n = length(idx)
            block = args_cluster[i:i+n-1]
            push!(blocks, block)
            i += n
        end

        # Sort blocks by length and count of Destroy/Create
        return sortperm(blocks, lt=_lt_num_destroy_create)
    end
end

function _lt_num_destroy_create(args1, args2)
    if length(args1) != length(args2) # Sort by length first
        return length(args1) > length(args2)
    else # Equal length is sorted by number of Destroy (lower ones to the right)
        f = x->isa(x,Destroy)
        idx_destroys1 = findall(f,args1)
        ndest1 = length(idx_destroys1)
        idx_destroys2 = findall(f,args2)
        ndest2 = length(idx_destroys2)
        return ndest1 <= ndest2
    end
end

function lt_reference_order(t1::Transition, t2::Transition)
    aon = acts_on(t1)
    isequal(aon,acts_on(t2)) && return false
    lvls = levels(t1.hilbert, aon)
    f = x->findfirst(isequal(x),lvls)
    i1, j1 = f(t1.i), f(t1.j)
    i2, j2 = f(t2.i), f(t2.j)
    d1 = abs(i1 - j1)
    d2 = abs(i2 - j2)
    if d1==d2
        if i1==i2
            return j1 <= j2
        else
            m1 = min(i1, j1)
            m2 = min(i2, j2)
            if m1==m2
                return i1 > i2
            else
                return m1 < m2
            end
        end
    else
        return d1 < d2
    end
end

_swap_aon_and_name(x::Average, aon1, aon2, names) = _average(_swap_aon_and_name(x.arguments[1], aon1, aon2, names))
function _swap_aon_and_name(t::QMul, aon1, aon2, names)
    args = QSym[]
    for arg ∈ t.args_nc
        idx = findfirst(isequal(acts_on(arg)), aon1)
        push!(args, _swap_aon_and_name(arg, aon1[idx], aon2[idx], names[idx]))
    end
    sort!(args, by=acts_on)
    return QMul(t.arg_c, args)
end
for T in (:Create, :Destroy)
    @eval _swap_aon_and_name(op::$(T), aon1, aon2, name::Symbol) = $(T)(op.hilbert, name, aon2; op.metadata)
end
function _swap_aon_and_name(op::Transition, aon1, aon2, name::Symbol)
    Transition(op.hilbert, name, op.i, op.j, aon2; op.metadata)
end
_swap_aon_and_name(op::QSym, aon1, aon2, names) = _swap_aon_and_name(op, aon1, aon2, _get_names(names, aon2))

_get_names(names, aon::Integer) = names[aon]
_get_names(names::Vector{<:Symbol}, aon::ClusterAon) = names[aon.j]
_get_names(names, aon::ClusterAon) = names[aon.i][aon.j]
function _get_names(names, aon)
    names_ = Symbol[]
    for a∈aon
        n = [_get_names(names, a1) for a1∈a]
        append!(names_, n)
    end
end

## Complete
function complete!(de::ScaledMeanfieldEquations;
                                order=de.order,
                                multithread=false,
                                filter_func=nothing,
                                mix_choice=maximum,
                                simplify=true,
                                kwargs...)
    vs = de.states
    order_lhs = maximum(get_order.(vs))
    order_rhs = 0
    for i=1:length(de.equations)
        k = get_order(de.equations[i].rhs)
        k > order_rhs && (order_rhs = k)
    end
    if order === nothing
        order_ = max(order_lhs, order_rhs)
    else
        order_ = order
    end
    maximum(order_) >= order_lhs || error("Cannot form cumulant expansion of derivative; you may want to use a higher order!")

    if order_ != de.order
        for i=1:length(de.equations)
            lhs = de.equations[i].lhs
            rhs = cumulant_expansion(de.equations[i].rhs,order_;
                                        mix_choice=mix_choice,
                                        simplify=simplify)
            de.equations[i] = Symbolics.Equation(lhs, rhs)
        end
    end

    names = de.names
    scale_aons = de.scale_aons

    vhash = map(hash, vs)
    vs′ = map(_conj, vs)
    vs′hash = map(hash, vs′)
    filter!(!in(vhash), vs′hash)
    missed = find_missing(de.equations, vhash, vs′hash; get_adjoints=false)
    filter_redundants!(missed,scale_aons,names,vhash,vs′hash)
    isnothing(filter_func) || filter!(filter_func, missed) # User-defined filter

    while !isempty(missed)
        ops_ = [SymbolicUtils.arguments(m)[1] for m in missed]
        he = meanfield(ops_,de.hamiltonian,de.jumps;
                                Jdagger=de.jumps_dagger,
                                rates=de.rates,
                                simplify=simplify,
                                multithread=multithread,
                                order=order_,
                                mix_choice=mix_choice,
                                iv=de.iv,
                                kwargs...)

        _append!(de, he)

        vhash_ = hash.(he.states)
        vs′hash_ = hash.(_conj.(he.states))
        append!(vhash, vhash_)
        for i=1:length(vhash_)
            vs′hash_[i] ∈ vhash_ || push!(vs′hash, vs′hash_[i])
        end

        missed = find_missing(he.equations, vhash, vs′hash; get_adjoints=false)
        filter_redundants!(missed,scale_aons,names,vhash,vs′hash)
        isnothing(filter_func) || filter!(filter_func, missed) # User-defined filter
    end

    if !isnothing(filter_func)
        # Find missing values that are filtered by the custom filter function,
        # but still occur on the RHS; set those to 0
        missed = find_missing(de)
        filter!(!filter_func, missed)
        filter_redundants!(missed,scale_aons,names)
        missed_adj = map(_adjoint, missed)
        subs = Dict(vcat(missed, missed_adj) .=> 0)
        for i=1:length(de.equations)
            de.equations[i] = substitute(de.equations[i], subs)
            de.states[i] = de.equations[i].lhs
        end
    end
    return de
end

function filter_redundants!(missed,scale_aons_::Vector{<:Vector},names,vhash=UInt[],vs′hash=UInt[])
    j = 1
    while j <= length(missed)
        for scale_aons∈scale_aons_
            missed[j] = substitute_redundants(missed[j],scale_aons,names)
        end
        h = hash(missed[j])
        h′ = hash(_adjoint(missed[j]))
        if h∈vhash || h′∈vhash || h∈vs′hash || h′∈vs′hash
            deleteat!(missed, j)
        else
            j += 1
            push!(vhash, h)
        end
    end
    return missed
end


text:
struct SklarDist{CT,TplMargins} <: Distributions.ContinuousMultivariateDistribution
    C::CT
    m::TplMargins
    function SklarDist(C,m)
        d = length(C)
        @assert length(m) == d
        @assert all(mᵢ isa Distributions.UnivariateDistribution for mᵢ in m)
        return new{typeof(C),typeof(m)}(C,m)
    end
end
Base.length(S::SklarDist{CT,TplMargins}) where {CT,TplMargins} = length(S.C)
Base.eltype(S::SklarDist{CT,TplMargins}) where {CT,TplMargins} = Base.eltype(S.C)

function Distributions.cdf(S::SklarDist{CT,TplMargins},x) where {CT,TplMargins}
    return cdf(S.C,cdf.(S.m,x))
end
function Distributions._rand!(rng::Distributions.AbstractRNG, S::SklarDist{CT,TplMargins}, x::AbstractVector{T}) where {CT,TplMargins,T}
    Random.rand!(rng,S.C,x)
     x .= Distributions.quantile.(S.m,x)
end
function Base.rand(rng::Distributions.AbstractRNG,S::SklarDist{CT,TplMargins}) where {CT,TplMargins}
     x = Random.rand(rng,S.C)
     x .= Distributions.quantile.(S.m,x)
     return x
end
function Distributions._logpdf(S::SklarDist{CT,TplMargins},u) where {CT,TplMargins}
    sum(Distributions.logpdf(S.m[i],u[i]) for i in 1:length(u)) + Distributions.logpdf(S.C,Distributions.cdf.(S.m,u))
end

function Distributions.fit(::Type{SklarDist{CT,TplMargins}},x) where {CT,TplMargins}
    # The first thing to do is to fit the marginals :
    @assert length(TplMargins.parameters) == size(x,1)
    m = Tuple(Distributions.fit(TplMargins.parameters[i],x[i,:]) for i in 1:size(x,1))
    u = pseudos(x)
    C = Distributions.fit(CT,u)
    return SklarDist(C,m)
end


text:
abstract type DensityFilter end
struct  LocalFilter{N <: Neighborhood , K<:Kernel} <: DensityFilter
	𝒩::N
	𝒦::K
end

struct HelmHoltzPDEFilter <: DensityFilter end

max_stars_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_issues_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_forks_repo_name: MKAbdElrahman/TopologyOptimizationTools


text:
module input
    include("../structs/Core.jl")

    function from_dataframe(source; name = "moeum.MOEUM")
        dat_sori, hol_sori = [String(item) for item in source.colindex.names], []

        datas = source.columns
        for i in eachindex(datas[1])
            item = []
            for col in datas
                append!(item, col[i])
            end

            append!(hol_sori, [item])
        end

        return Core.MOEUM(dat_sori = dat_sori, hol_sori = hol_sori, name = name)
    end

    function from_dict(source ; name::String = "moeum.MOEUM")
        convertable, result = true, nothing
        if length(source) > 0 && isa(source, Array)
            for item in source
                if !isa(item, Dict)
                    convertable = false
                    break
                end
            end

            if convertable
                dat_sori = [key for key in keys(source[1])]
                hol_sori = [[row[dat] for dat in dat_sori] for row in source]

                result = Core.MOEUM(dat_sori = dat_sori, hol_sori = hol_sori, name = name)
            else
                result = convertable
            end
        elseif isa(source, Dict)
            for value in values(source)
                if !isa(value, Array)
                    convertable = false
                    break
                elseif length(value) == 0
                    convertable = false
                    break
                end
            end

            if convertable
                dat_sori = [key for key in keys(source)]
                hol_sori = [[source[dat][i] for dat in dat_sori] for i in eachindex(source[dat_sori[1]])]

                result = Core.MOEUM(dat_sori = dat_sori, hol_sori = hol_sori, name = name)
            else
                result = convertable
            end
        end

        return result
    end

    function from_csv(source_file; name = "")
        result = nothing
        name = name == "" ? splitext(basename(source_file))[1] : name
        open(source_file, "r") do f_csv
            lines = readlines(f_csv)
            dat_sori, hol_sori = [strip(item) for item in split(lines[1], ",")], [[strip(item) for item in split(line, ",")] for line in lines[2:end]]

            result = Core.MOEUM(dat_sori = dat_sori, hol_sori = hol_sori, name = name)
        end

        return result
    end
end


text:
using Test
include("Parse.jl")
out = readinput("../src/testing/superminimal.yml")
r = out["gas"]["Reactions"][1]
s = r.products[1]

@test r.kinetics(1000.0) ≈ 25618.8 rtol=1e-4
@test getEnthalpy(s.thermo,1000.0) ≈ 20677.0 rtol=1e-4
@test getEnthalpy(out["gas"]["Species"][5].thermo,1000.0) ≈ 20677.2 rtol=1e-4


text:
import Base.==

==(a::FactNode, b::FactNode) = a.line == b.line && a.name == b.name &&a.result == b.result && a.details == b.details
==(a::FactsCollectionNode, b::FactsCollectionNode) = a.line == b.line && a.name == b.name && a.children == b.children
==(a::ContextNode, b::ContextNode) = a.line == b.line && a.name == b.name && a.children == b.children
==(a::RootNode, b::RootNode) = a.children == b.children


text:
@with_kw struct SimulatorOptions{T}
    warmstart::Bool=false
    record::Bool=true
    z_warmstart::T=0.001
    κ_warmstart::T=0.001
	failure_abort::Int=50
end

struct SimulatorStatistics{T}
    policy_time::Vector{T}
    policy_mean::Vector{T}
    policy_std::Vector{T}
    sim_time::Vector{T}
    sim_mean::Vector{T}
    sim_std::Vector{T}
end

function SimulatorStatistics(H)
    policy_time = zeros(H)
    policy_mean = zeros(1)
    policy_std = zeros(1)
    sim_time = zeros(H)
    sim_mean = zeros(1)
    sim_std = zeros(1)
    return SimulatorStatistics(policy_time, policy_mean, policy_std,
        sim_time, sim_mean, sim_std)
end

function process!(stats::SimulatorStatistics, N_sample::Int)
    # policy
    H_sim = length(stats.policy_time)
    H = Int(H_sim / N_sample)
    policy_time = sum(reshape(stats.policy_time, (H, N_sample)), dims=2)
    stats.policy_mean[1] = mean(policy_time)
    stats.policy_std[1] = sqrt(mean((policy_time .- stats.policy_mean[1]).^2.0))

    # simulation
    H_sim = length(stats.sim_time)
    H = Int(H_sim / N_sample)
    sim_time = sum(reshape(stats.sim_time, (H, N_sample)), dims=2)
    stats.sim_mean[1] = mean(sim_time)
    stats.sim_std[1] = sqrt(mean((sim_time .- stats.sim_mean[1]).^2.0))

    return nothing
end

struct Simulator{T,R,RZ,Rθ,M<:Model{T},P<:Policy{T},D<:Disturbances{T}}
    model::M
    policy::P
    dist::D
    traj::Trajectory{T}
    grad::GradientTrajectory{T}
    ip::InteriorPoint{T,R,RZ,Rθ}
    idx_z::IndicesZ
    idx_θ::Indicesθ
    f::Vector{T}
    h::T
    stats::SimulatorStatistics{T}
    opts::SimulatorOptions{T}
end

function Simulator(model, T;
        h=0.01,
        policy=empty_policy(model),
        dist=empty_disturbances(model),
        f=friction_coefficients(model),
        residual=eval(residual_expr(model)),
        jacobian_z=eval(jacobian_var_expr(model)),
        jacobian_θ=eval(jacobian_data_expr(model)),
        diff_sol=false,
        solver_opts=InteriorPointOptions(
            undercut=Inf,
            γ_reg=0.1,
            r_tol=1e-8,
            κ_tol=1e-8,
            max_ls=25,
            ϵ_min=0.25,
            diff_sol=diff_sol,
            verbose=false),
        stats=SimulatorStatistics(T),
        sim_opts=SimulatorOptions())

    idx_z = indices_z(model)
    idx_θ = indices_θ(model, nf=length(f))
    idx_opt = indices_optimization(model)

    nz = num_var(model)
    nθ = num_data(model, nf=length(f))
    z0 = zeros(nz)
    θ0 = zeros(nθ)
    q0 = nominal_configuration(model)
    initialize_z!(z0, model, idx_z, q0)
    initialize_θ!(θ0, model, idx_θ, q0, q0, zeros(model.nu), zeros(model.nw), f, h)

    ip = interior_point(
         z0,
         θ0,
         idx=idx_opt,
         r! = residual,
         rz! = jacobian_z,
         rθ! = jacobian_θ,
         rz=zeros(nz, nz),
         rθ=zeros(nz, nθ),
         opts=solver_opts)

    traj = Trajectory(model, T)
    grad = GradientTrajectory(model, T)

    Simulator(model, policy, dist, traj, grad, ip, idx_z, idx_θ, f, h, stats, sim_opts)
end

function step!(s::Simulator{T}, t::Int; verbose=false) where T
    model = s.model
    traj = s.traj
    grad = s.grad
    ip = s.ip
    idx_z = s.idx_z
    idx_θ = s.idx_θ
    f = s.f
    h = s.h

    # initialize
    initialize_z!(ip.z, model, idx_z, traj.q[t+1])
    initialize_θ!(ip.θ, model, idx_θ, traj.q[t], traj.q[t+1], traj.u[t], traj.w[t], f, h)

    # solve
    sim_time = @elapsed status = interior_point_solve!(ip)
    s.opts.record && (s.stats.sim_time[t] = sim_time)

    # status check
    if !status
        verbose && (@show norm(ip.r))
        verbose && (@warn "residual failure")
        return false
    end

    # cache solution
    solution!(traj, ip.z, idx_z, h, t)

    # cache gradients
    ip.opts.diff_sol && gradient!(grad, h, ip.δz, idx_z, idx_θ, t)

    return status
end

function step!(s::Simulator{T}, q::AbstractVector{T}, v::AbstractVector{T}, u::AbstractVector{T}, t::Int; verbose=false) where T
    # set state
    set_state!(s, q, v, t)

    # set control
    s.traj.u[t] .= u

    # step
    step!(s, t, verbose=verbose)

    return s.traj.q[t+2] #TODO: return gradients
end

function solution!(traj::Trajectory{T}, z::Vector{T}, idx_z::IndicesZ, h::T, t::Int) where T
    q = @views z[idx_z.q]
    γ = @views z[idx_z.γ]
    b = @views z[idx_z.b]

    # current configuration
    traj.q[t+2] .= q

    # current velocity
    traj.v[t+1] .= traj.q[t+2]
    traj.v[t+1] .-= traj.q[t+1]
    traj.v[t+1] ./= h

    # contact impulses
    traj.γ[t] .= γ
    traj.b[t] .= b

    return nothing
end

function gradient!(grad::GradientTrajectory{T}, h::T, δz::Matrix{T}, idx_z::IndicesZ, idx_θ::Indicesθ, t::Int) where T
    ∂q3∂q1 = @views δz[idx_z.q, idx_θ.q1]
    ∂q3∂q2 = @views δz[idx_z.q, idx_θ.q2]
    ∂q3∂u1 = @views δz[idx_z.q, idx_θ.u]
    ∂γ1∂q1 = @views δz[idx_z.γ, idx_θ.q1]
    ∂γ1∂q2 = @views δz[idx_z.γ, idx_θ.q2]
    ∂γ1∂u1 = @views δz[idx_z.γ, idx_θ.u]
    ∂b1∂q1 = @views δz[idx_z.b, idx_θ.q1]
    ∂b1∂q2 = @views δz[idx_z.b, idx_θ.q2]
    ∂b1∂u1 = @views δz[idx_z.b, idx_θ.u]

    grad.∂q3∂q1[t] .= ∂q3∂q1
    grad.∂q3∂q2[t] .= ∂q3∂q2
    grad.∂q3∂u1[t] .= ∂q3∂u1

    grad.∂γ1∂q1[t] .= ∂γ1∂q1
    grad.∂γ1∂q2[t] .= ∂γ1∂q2
    grad.∂γ1∂u1[t] .= ∂γ1∂u1

    grad.∂b1∂q1[t] .= ∂b1∂q1
    grad.∂b1∂q2[t] .= ∂b1∂q2
    grad.∂b1∂u1[t] .= ∂b1∂u1

    grad.∂q3∂v1[t] .= ∂q3∂q2
    grad.∂q3∂v1[t] .-= ∂q3∂q1
    grad.∂q3∂v1[t] ./= h

    grad.∂γ1∂v1[t] .= ∂γ1∂q2
    grad.∂γ1∂v1[t] .-= ∂γ1∂q1
    grad.∂γ1∂v1[t] ./= h

    grad.∂b1∂v1[t] .= ∂b1∂q2
    grad.∂b1∂v1[t] .-= ∂b1∂q1
    grad.∂b1∂v1[t] ./= h

    return nothing
end

function set_state!(s::Simulator{T}, q::AbstractArray, v::AbstractArray, t::Int) where T
    # initial configuration and velocity
    s.traj.q[t+1] .= q # q2
    s.traj.v[t] .= v   # v1

    # q1
    s.traj.q[t] .= v
    s.traj.q[t] .*= s.h
    s.traj.q[t] .*= -1.0
    s.traj.q[t] .+= q

    return nothing
end

function simulate!(s::Simulator{T}, q::AbstractVector{T}, v::AbstractVector{T};
        reset_traj=false, verbose=false) where T
    # reset trajectory
    reset_traj && reset!(s.traj)
    reset_traj && reset!(s.grad)

    # reset solver

    # set initial state
    set_state!(s, q, v, 1)

    # simulate
    simulate!(s, verbose=verbose)
end

function simulate!(s::Simulator{T}; verbose=false) where T
    status = false

    N = length(s.traj.u)
    p = s.policy
    w = s.dist
    traj = s.traj

    for t = 1:N
        # policy
        policy_time = @elapsed traj.u[t] .= policy(p, traj, t)
        s.opts.record && (s.stats.policy_time[t] = policy_time)

        # disturbances
        traj.w[t] .= disturbances(w, traj.q[t+1], t)

        # step
        status = step!(s, t, verbose=verbose)
        !status && break
    end
    return status
end

function visualize!(vis, s::Simulator; skip=1, fixed_camera=true)
    visualize!(vis, s.model, s.traj.q[1:skip:end], Δt=skip * s.h, fixed_camera=fixed_camera)
end


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_energy_market_dp" begin

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     |
    DP => 10h40         |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DP => 10h30         |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DP => 10h45        |
                        |

    =#
    function create_instance(ech, next_ech, ts,
                            wind_1_level,
                            prod_1_level,
                            prod_2_level;
                            fo_length=Minute(5))
        TS = [ts]
        ECH = [DateTime("2015-01-01T10:00:00"), DateTime("2015-01-01T10:40:00")]
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(20*60), Dates.Second(20*60))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                10., 100.,
                                                0., 10.,
                                                Dates.Second(3*60*60), Dates.Second(30*60))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                10., 100.,
                                                0., 15.,
                                                Dates.Second(3*60*60), Dates.Second(15*60))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ECH[1], "wind_1_1", TS[1], "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[1], "bus_1", TS[1], "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[2], "wind_1_1", TS[1], "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ECH[2], "bus_1", TS[1], "S1", 55.)
        # initial generators state : need to pay starting cost at TS[1]
        generators_init_state = SortedDict(
            "prod_1_1" => PSCOPF.OFF,
            "prod_1_2" => PSCOPF.OFF
        )
        #ManagementMode
        mode = PSCOPF.ManagementMode("test_mode", fo_length)

        context = PSCOPF.PSCOPFContext(network, TS, mode,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        firmness = PSCOPF.compute_firmness(market,
                                            ech, next_ech,
                                            TS, context)

        prod_1_definitive_state = (prod_1_level < 1e-09) ? PSCOPF.OFF : PSCOPF.ON
        prod_1_definitive_level = (PSCOPF.get_power_level_firmness(firmness, "prod_1_1", ts) == PSCOPF.DECIDED) ? prod_1_level : missing
        prod_2_definitive_state = (prod_2_level < 1e-09) ? PSCOPF.OFF : PSCOPF.ON
        prod_2_definitive_level = (PSCOPF.get_power_level_firmness(firmness, "prod_1_2", ts) == PSCOPF.DECIDED) ? prod_2_level : missing
        context.tso_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T06:00:00"), SortedDict(
                "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                    SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                    # SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing,
                    #                                                                         SortedDict("S1"=>PSCOPF.ON))),
                    SortedDict(ts => PSCOPF.UncertainValue{Float64}(missing,
                                                                    SortedDict("S1"=>wind_1_level)))
                    ),
                "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                    SortedDict(ts => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(prod_1_definitive_state,
                                                                    SortedDict("S1"=>prod_1_definitive_state))),
                    SortedDict(ts => PSCOPF.UncertainValue{Float64}(prod_1_definitive_level,
                                                                    SortedDict("S1"=>prod_1_level)))
                    ),
                "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                    SortedDict(ts => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(prod_2_definitive_state,
                                                                    SortedDict("S1"=>prod_2_definitive_state))),
                    SortedDict(ts => PSCOPF.UncertainValue{Float64}(prod_2_definitive_level,
                                                                    SortedDict("S1"=>prod_2_level)))
                    )
                )
            )

        return context, market, firmness
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     |
    DP => 10h40         |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DP => 10h30         |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DP => 10h45        |
                        |

    ECH1*         ECH2            ECH3                                 <---FO--->TS
    |             |               |
    10h           10h30           10h40          10h45                          11h
                  <------------------------------------------DP(prod1)----------->
                                  <--------------------------DP(wind1)----------->
                                                   <---------DP(prod2)----------->

    We can change the production level of a unit before reaching the DP
    ech = 10h
    TS = 11h
    DP(prod_1_1) = 30 mins => 10h30
    DP(prod_1_2) = 15 mins => 10h45

    We suppose that in the preceding step we decided :
    prod_1_1 : 0. (OFF)
    prod_1_2 : 37. (ON)

    In terms of cost, it is cheaper to use prod_1_1. But, since we are past DMO, we can no longer start prod_1_1.
    wind_1_1 : 20 MW
    remaining Demand : 35MW
    Since we haven't reached the DP yet, we can change the production level of prod_1_2 from 37 MW to 35 MW
    => the new decision :
    wind_1_1 : 20 MW
    prod_1_1 : OFF (unchanged)
    prod_1_2 : 35. (decreased from 37MW)
    =#
    @testset "energy_market_can_decrease_the_production_level_before_dp_if_ON" begin
        ech = DateTime("2015-01-01T10:00:00")
        next_ech = DateTime("2015-01-01T10:30:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    18., 0., 37.,
                                                    )
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # prod_1_1 is OFF cause it already was OFF:
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        # prod_1_2 was ON, it can change it's production : 37 => 35.
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     |
    DP => 10h40         |
      S1: 20            |
                        |
    (pilotable) prod_1_1|
    Pmin=10, Pmax=100   |
    Csta=0, Cprop=10    |
    DP => 10h30         |
                        |
    (pilotable) prod_1_2|
     Pmin=10, Pmax=100  |
     Csta=0, Cprop=15   |
     DP => 10h45        |
                        |

    ECH1*         ECH2            ECH3                                 <---FO--->TS
    |             |               |
    10h           10h30           10h40          10h45                          11h
                  <------------------------------------------DP(prod1)----------->
                                  <--------------------------DP(wind1)----------->
                                                   <---------DP(prod2)----------->

    We can change the production level of a unit before reaching the DP
    ech = 10h
    TS = 11h
    DP(prod_1_1) = 30 mins => 10h30
    DP(prod_1_2) = 15 mins => 10h45

    We suppose that in the preceding step we decided :
    prod_1_1 : 0. (OFF)
    prod_1_2 : 33. (ON)

    In terms of cost, it is cheaper to use prod_1_1. But, since we are past DMO, we can no longer start prod_1_1.
    wind_1_1 : 20 MW
    remaining Demand : 35MW
    Since we haven't reached the DP yet, we can change the production level of prod_1_2 from 33 MW to 35 MW
    => the new decision :
    wind_1_1 : 20 MW
    prod_1_1 : OFF (unchanged)
    prod_1_2 : 35. (increased from 33MW)
    =#
    @testset "energy_market_can_increase_the_production_level_before_dp_if_ON" begin
        ech = DateTime("2015-01-01T10:00:00")
        next_ech = DateTime("2015-01-01T10:30:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    22., 0., 33.,
                                                    )
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # prod_1_1 is OFF cause it already was OFF:
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        # prod_1_2 was ON, it can change it's production : 33 => 35.
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
    end

    #=
    ECH1          ECH2            ECH3*          ECH4                  <---FO--->TS
    |             |               |
    10h           10h30           10h40          10h45                          11h
                  <------------------------------------------DP(prod1)----------->
                                  <--------------------------DP(wind1)----------->
                                                   <---------DP(prod2)----------->

    We can no longer change the production level of a unit once DP is reached
    ech = 10h40
    TS = 11h
    DP(prod_1_1) = 30 mins => 10h30 => we are past DP
    DP(prod_1_2) = 15 mins => 10h45

    We suppose that in the preceding step we decided :
    prod_1_1 : 17. (ON)
    prod_1_2 : 20. (ON)

    In terms of cost, it is cheaper to use prod_1_1.
    But, since we are past DP, we can no longer increase the level of prod_1_1.

    wind_1_1 : 20 MW
        remaining Demand : 35MW
    prod_1_1 : 17. (unchanged)
    prod_1_2 : 18. (decreased from 20MW)
    =#
    @testset "energy_market_cannot_change_the_production_level_after_dp" begin
        # We are past DP[prod_1_1]
        ech = DateTime("2015-01-01T10:40:00")
        next_ech = DateTime("2015-01-01T10:45:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    18., 17., 20.,
                                                    )

        expected_firmness = PSCOPF.Firmness(
                                    SortedDict("prod_1_1" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                ),
                                    SortedDict("wind_1_1" => SortedDict(TS[1] => PSCOPF.TO_DECIDE),
                                                "prod_1_1" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.FREE),
                                                )
            )
        @test firmness == expected_firmness

        # prod_1_1 has a definitive value since we are past the DP

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # we are past DP[prod_1_] => we can't change the level of prod_1_1
        # prod_1_1 : 17.
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test 17. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        # prod_1_2 : 20 => 18.
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 18. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
    end

    #=
    at ECH2, we are at wind_1_1's DP,
    but the decision made will still be by scenario cause it depends on the uncertainties
    limitables's DP only applies to limitation decisions (=> applied in TSO)
    ECH1          ECH2            ECH3*          ECH4                  <---FO--->TS
    |             |               |
    10h           10h30           10h40          10h45                          11h
                  <------------------------------------------DP(prod1)----------->
                                  <--------------------------DP(wind1)----------->
                                                   <---------DP(prod2)----------->
    =#
    @testset "energy_market_wind_is_always_uncertain" begin
        # We are past DP[prod_1_1]
        ech = DateTime("2015-01-01T10:40:00")
        next_ech = DateTime("2015-01-01T10:45:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    18., 17., 20.,
                                                    )

        expected_firmness = PSCOPF.Firmness(
                                    SortedDict("prod_1_1" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                ),
                                    SortedDict("wind_1_1" => SortedDict(TS[1] => PSCOPF.TO_DECIDE),
                                                "prod_1_1" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.FREE),
                                                )
            )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test !PSCOPF.is_definitive(PSCOPF.get_prod_uncertain_value(context.market_schedule, "wind_1_1", TS[1]))
    end

    #=
    ECH1*         ECH2            ECH3           ECH4                  <---FO--->TS
    |             |               |
    10h           10h30           10h40          10h45                          11h
                  <------------------------------------------DP(prod1)----------->
                                  <--------------------------DP(wind1)----------->
                                                   <---------DP(prod2)----------->
    =#
    @testset "energy_market_cannot_change_the_production_level_before_dp_if_unit_is_off_and_past_DMO" begin
        ech = DateTime("2015-01-01T10:00:00")
        next_ech = DateTime("2015-01-01T10:30:00")
        TS = [DateTime("2015-01-01T11:00:00")]
        context, market, firmness = create_instance(ech, next_ech, TS[1],
                                                    35., 0., 25.,
                                                    )

        expected_firmness = PSCOPF.Firmness(
                                    SortedDict("prod_1_1" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.DECIDED),
                                                ),
                                    SortedDict("wind_1_1" => SortedDict(TS[1] => PSCOPF.FREE),
                                                "prod_1_1" => SortedDict(TS[1] => PSCOPF.FREE),
                                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.FREE),
                                                )
            )
        @test firmness == expected_firmness

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        # We are past DMO => unit state was already fixed
        # prod_1_1 : OFF
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        # prod_1_2 : 25 => 35.
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
abstract type AbstractCN0Estimator end

"""
$(SIGNATURES)

MomentsCN0Estimator to estimate the CN0
"""
struct MomentsCN0Estimator{N} <: AbstractCN0Estimator
    prompt_buffer::SVector{N, ComplexF64}
    current_index::Int
    length::Int
end

function MomentsCN0Estimator(N)
    MomentsCN0Estimator{N}(zero(SVector{N, ComplexF64}), 0, 0)
end

length(cn0_state::MomentsCN0Estimator) = cn0_state.length
get_prompt_buffer(cn0_state::MomentsCN0Estimator) = cn0_state.prompt_buffer
get_current_index(cn0_state::MomentsCN0Estimator) = cn0_state.current_index

"""
$(SIGNATURES)

Updates the `cn0_estimator` to include the information of the current prompt value.
"""
function update(cn0_estimator::MomentsCN0Estimator{N}, prompt) where N
    next_index = mod(get_current_index(cn0_estimator), N) + 1
    next_prompt_buffer = setindex(get_prompt_buffer(cn0_estimator), prompt, next_index)
    next_length = min(length(cn0_estimator) + 1, N)
    MomentsCN0Estimator{N}(next_prompt_buffer, next_index, next_length)
end

"""
$(SIGNATURES)

Estimates the CN0 based on the struct `cn0_estimator`.
"""
function estimate_cn0(cn0_estimator::MomentsCN0Estimator, integration_time)
    length(cn0_estimator) == 0 && return 0.0dBHz
    abs2_prompt_buffer = abs2.(get_prompt_buffer(cn0_estimator))
    M₂ = 1 / length(cn0_estimator) * sum(abs2_prompt_buffer)
    M₄ = 1 / length(cn0_estimator) * sum(abs2_prompt_buffer .^ 2)
    Pd = sqrt(abs(2 * M₂^2 - M₄))
    noise_power = (M₂ - Pd)
    noise_power_non_neg = noise_power - 2 * (noise_power < 0) * noise_power
    SNR = Pd / noise_power_non_neg
    dBHz(SNR / integration_time)
end


text:
module MaxHelpingHandGroupedTriggerSpikes

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/GroupedTriggerSpikesUp" GroupedTriggerSpikesUp(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, type::String="default", behindMoveBlocks::Bool=false, triggerIfSameDirection::Bool=false)
@mapdef Entity "MaxHelpingHand/GroupedTriggerSpikesDown" GroupedTriggerSpikesDown(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, type::String="default", behindMoveBlocks::Bool=false, triggerIfSameDirection::Bool=false)
@mapdef Entity "MaxHelpingHand/GroupedTriggerSpikesLeft" GroupedTriggerSpikesLeft(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, type::String="default", behindMoveBlocks::Bool=false, triggerIfSameDirection::Bool=false)
@mapdef Entity "MaxHelpingHand/GroupedTriggerSpikesRight" GroupedTriggerSpikesRight(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, type::String="default", behindMoveBlocks::Bool=false, triggerIfSameDirection::Bool=false)

const placements = Ahorn.PlacementDict()

groupedTriggerSpikes = Dict{String, Type}(
    "up" => GroupedTriggerSpikesUp,
    "down" => GroupedTriggerSpikesDown,
    "left" => GroupedTriggerSpikesLeft,
    "right" => GroupedTriggerSpikesRight
)

groupedTriggerSpikesArrows = Dict{String, String}(
    "up" => "↑",
    "down" => "↓",
    "left" => "←",
    "right" => "→"
)

const spikeTypes = String[
    "default",
    "outline",
    "cliffside",
    "dust",
    "reflection"
]

const triggerSpikeColors = [
    (242, 90, 16, 255) ./ 255,
    (255, 0, 0, 255) ./ 255,
    (242, 16, 103, 255) ./ 255
]

groupedTriggerSpikesUnion = Union{GroupedTriggerSpikesUp, GroupedTriggerSpikesDown, GroupedTriggerSpikesLeft, GroupedTriggerSpikesRight}

for variant in spikeTypes
    for (dir, entity) in groupedTriggerSpikes
        key = "Grouped Trigger Spikes ($(groupedTriggerSpikesArrows[dir]), $(uppercasefirst(variant))) (max480's Helping Hand)"
        placements[key] = Ahorn.EntityPlacement(
            entity,
            "rectangle",
            Dict{String, Any}(
                "type" => variant
            )
        )
    end
end

Ahorn.editingOptions(entity::groupedTriggerSpikesUnion) = Dict{String, Any}(
    "type" => spikeTypes
)

directions = Dict{String, String}(
    "MaxHelpingHand/GroupedTriggerSpikesUp" => "up",
    "MaxHelpingHand/GroupedTriggerSpikesDown" => "down",
    "MaxHelpingHand/GroupedTriggerSpikesLeft" => "left",
    "MaxHelpingHand/GroupedTriggerSpikesRight" => "right",
)

offsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (4, -4),
    "down" => (4, 4),
    "left" => (-4, 4),
    "right" => (4, 4),
)

groupedTriggerSpikesOffsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (0, 5),
    "down" => (0, -4),
    "left" => (5, 0),
    "right" => (-4, 0),
)

rotations = Dict{String, Number}(
    "up" => 0,
    "right" => pi / 2,
    "down" => pi,
    "left" => pi * 3 / 2
)

triggerRotationOffsets = Dict{String, Tuple{Number, Number}}(
    "up" => (3, -1),
    "right" => (4, 3),
    "down" => (5, 5),
    "left" => (-1, 4),
)

resizeDirections = Dict{String, Tuple{Bool, Bool}}(
    "up" => (true, false),
    "down" => (true, false),
    "left" => (false, true),
    "right" => (false, true),
)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::groupedTriggerSpikesUnion)
    direction = get(directions, entity.name, "up")
    theta = rotations[direction] - pi / 2

    width = Int(get(entity.data, "width", 0))
    height = Int(get(entity.data, "height", 0))

    x, y = Ahorn.position(entity)
    cx, cy = x + floor(Int, width / 2) - 8 * (direction == "left"), y + floor(Int, height / 2) - 8 * (direction == "up")

    Ahorn.drawArrow(ctx, cx, cy, cx + cos(theta) * 24, cy + sin(theta) * 24, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        x, y = Ahorn.position(entity)

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        direction = get(directions, entity.name, "up")

        ox, oy = offsets[direction]

        return Ahorn.Rectangle(x + ox - 4, y + oy - 4, width, height)
    end
end

Ahorn.minimumSize(entity::groupedTriggerSpikesUnion) = 8, 8

function Ahorn.resizable(entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        direction = get(directions, entity.name, "up")

        return resizeDirections[direction]
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "type", "default")
        direction = get(directions, entity.name, "up")
        groupedTriggerSpikesOffset = groupedTriggerSpikesOffsets[direction]

        if variant == "dust"
            rng = Ahorn.getSimpleEntityRng(entity)

            width = get(entity.data, "width", 8)
            height = get(entity.data, "height", 8)

            updown = direction == "up" || direction == "down"

            for ox in 0:8:width - 8, oy in 0:8:height - 8
                color1 = rand(rng, triggerSpikeColors)
                color2 = rand(rng, triggerSpikeColors)

                drawX = ox + triggerRotationOffsets[direction][1]
                drawY = oy + triggerRotationOffsets[direction][2]

                Ahorn.drawSprite(ctx, "danger/triggertentacle/wiggle_v06", drawX, drawY, rot=rotations[direction], tint=color1)
                Ahorn.drawSprite(ctx, "danger/triggertentacle/wiggle_v03", drawX + 3 * updown, drawY + 3 * !updown, rot=rotations[direction], tint=color2)
            end
        else

            width = get(entity.data, "width", 8)
            height = get(entity.data, "height", 8)

            for ox in 0:8:width - 8, oy in 0:8:height - 8
                drawX, drawY = (ox, oy) .+ offsets[direction] .+ groupedTriggerSpikesOffset
                Ahorn.drawSprite(ctx, "danger/spikes/$(variant)_$(direction)00", drawX, drawY)
            end
        end
    end
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
include("intervalbox.jl")
include("setdiff.jl")


text:
factor_residual_variance(fm::FactorModel) = sum(fm.factor_residuals.^2)/apply(*, size(fm.x))  # see page 201 of Bai Ng 2002
#factor_residual_variance(fm::FactorModel) = sum(mapslices(x->x'x/length(x), fm.factor_residuals, 1))/size(fm.x, 2)  # the same as above
# and var(fm.factor_residuals) is approximately the same as well but it probably corrects for small sample bias

function criterion_cumulative_variance(pca_result, threshold=0.95)   # simply use factors until a certain threshold of variance is reached
    pca_result.cumulative_variance .< threshold  # take threshold% of variance
end

# PCp criteria as defined on page 201 of Bai and Ng 2002
function criterion_PCp1(fm::FactorModel; kmax=8)  # be careful! kmax is set as in Bai, Ng 2002 which might make no sense in a different setting!
    fm_unrestricted = FactorModel(fm.x, kmax)
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    factor_residual_variance(fm) + fm.number_of_factors*factor_residual_variance(fm_unrestricted)*(N_plus_T_by_NT)*log(N_plus_T_by_NT^-1)
end
function criterion_PCp2(fm::FactorModel; kmax=8)
    fm_unrestricted = FactorModel(fm.x, kmax)  # kmax is set as in Bai, Ng 2002
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    factor_residual_variance(fm) + sum(fm.number_of_factors)*factor_residual_variance(fm_unrestricted)*(N_plus_T_by_NT)*log(minimum(size(fm.x)))
end
function criterion_PCp3(fm::FactorModel; kmax=8)
    fm_unrestricted = FactorModel(fm.x, kmax)   # kmax is set as in Bai, Ng 2002
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    factor_residual_variance(fm) + sum(fm.number_of_factors)*factor_residual_variance(fm_unrestricted)*log(minimum(size(fm.x)))/minimum(size(fm.x))
end

# ICp criteria as defined on page 201 of Bai and Ng 2002
function criterion_ICp1(fm::FactorModel)
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    log(factor_residual_variance(fm)) + sum(fm.number_of_factors)*(N_plus_T_by_NT)*log(N_plus_T_by_NT^-1)
end
function criterion_ICp2(fm::FactorModel)
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    log(factor_residual_variance(fm)) + sum(fm.number_of_factors)*(N_plus_T_by_NT)*log(minimum(size(fm.x)))
end
function criterion_ICp3(fm::FactorModel)
    N_plus_T_by_NT = apply(+, size(fm.x))/apply(*, size(fm.x))
    log(factor_residual_variance(fm)) + sum(fm.number_of_factors)*log(minimum(size(fm.x)))/minimum(size(fm.x))
end

#  default criteria which are not quite consistent at least not when the number of factors is estimated
function criterion_BIC(fm::FactorModel)
    T, N = size(fm.x)
    factor_residual_variance(fm) + fm.number_of_factors * log(T)/T
end