text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

function update_parameter_dictionary(parameter_dictionary::Dict{String,Any}, varargs...)::Union{Nothing, Dict{String,Any}}

    # Default: return the unmodified dictionary -
    @info "Default implementation of update_parameter_dictionary returns unmodified dictionary"

    # return parameter dictionary -
    return parameter_dictionary
end

function generate_default_parameter_dictionary()::Dict{String,Any}

    # initialize -
    parameter_dictionary = Dict{String,Any}()
    system_type_flag = 	:CF_TXTL

    try

        # load the stoichiometric_matrix (SM) and degradation_dilution_matrix (DM) -
        SM = readdlm(joinpath(_PATH_TO_NETWORK,"Network.dat"))
        DM = readdlm(joinpath(_PATH_TO_NETWORK,"Degradation.dat"))

        # build the species initial condition array -
        initial_condition_array = [
			5.0	;	#	1	gene_venus	units: nM
			0.0	;	#	2	mRNA_venus	units: nM
			0.0	;	#	3	P_venus	units: nM

			# translation capacity -
			100.0	;	#	4	translation capacity	units: NA
		]

        # setup the system dimension -

		number_of_transcription_processes = 1
		number_of_translation_processes = 1

        # build the system species concentration array -
        system_concentration_array = [
			0.07	;	#	RNAP	units: µM
			2.3	;	#	RIBOSOME	units: µM
			0.035	;	#	P_σ70	units: µM
		]

        # get the biophysical parameters for this system type -
        biophysical_parameters_dictionary = get_biophysical_parameters_dictionary_for_system(system_type_flag)

        # setup the model parameter array -
        model_parameter_array = [

			# control parameters: gene_venus
			0.001	;	#	1	W_gene_venus
			1.0		;	#	2	W_gene_venus_P_σ70
			1.0		;	#	3	K_gene_venus_P_σ70
			1.0		;	#	4	n_gene_venus_P_σ70

			# transcription kinetic limit parameters: gene_venus
			1.0		;	#	5	K_gene_venus
			1.0		;	#	6	𝛕_gene_venus

			# translation kinetic limit parameters: P_venus
			1.0		;	#	7	K_P_venus
			1.0		;	#	8	𝛕_P_venus

			# species degradation constants -
			0.001		;	#	9	𝛳_mRNA_venus
			0.0001		;	#	10	𝛳_P_venus
		]

        # setup the parameter symbol index map -
        model_parameter_symbol_index_map = Dict{Symbol,Int}()
		model_parameter_symbol_index_map[:W_gene_venus] = 1
		model_parameter_symbol_index_map[:W_gene_venus_P_σ70] = 2
		model_parameter_symbol_index_map[:K_gene_venus_P_σ70] = 3
		model_parameter_symbol_index_map[:n_gene_venus_P_σ70] = 4
		model_parameter_symbol_index_map[:K_gene_venus] = 5
		model_parameter_symbol_index_map[:𝛕_gene_venus] = 6
		model_parameter_symbol_index_map[:K_P_venus] = 7
		model_parameter_symbol_index_map[:𝛕_P_venus] = 8
		model_parameter_symbol_index_map[:𝛳_mRNA_venus] = 9
		model_parameter_symbol_index_map[:𝛳_P_venus] = 10

        # setup the inverse parameter symbol index map -
        inverse_model_parameter_symbol_index_map = Dict{Int, Symbol}()
		inverse_model_parameter_symbol_index_map[1] = :W_gene_venus
		inverse_model_parameter_symbol_index_map[2] = :W_gene_venus_P_σ70
		inverse_model_parameter_symbol_index_map[3] = :K_gene_venus_P_σ70
		inverse_model_parameter_symbol_index_map[4] = :n_gene_venus_P_σ70
		inverse_model_parameter_symbol_index_map[5] = :K_gene_venus
		inverse_model_parameter_symbol_index_map[6] = :𝛕_gene_venus
		inverse_model_parameter_symbol_index_map[7] = :K_P_venus
		inverse_model_parameter_symbol_index_map[8] = :𝛕_P_venus
		inverse_model_parameter_symbol_index_map[9] = :𝛳_mRNA_venus
		inverse_model_parameter_symbol_index_map[10] = :𝛳_P_venus

        # specific growth rate (default units: h^-1)
        μ = 0.0 # default units: h^-1

        # translation capacity half-life (default units: h^-1) -
		half_life_translation_capacity=value(biophysical_parameters_dictionary,string(:half_life_translation_capacity))	# units: h^-1
		half_life_translation_capacity = half_life_translation_capacity*(3600)

        # == DO NOT EDIT BELOW THIS LINE ========================================================== #
        parameter_dictionary["initial_condition_array"] = initial_condition_array
        parameter_dictionary["number_of_states"] = length(initial_condition_array)
        parameter_dictionary["number_of_transcription_processes"] = number_of_transcription_processes
        parameter_dictionary["number_of_translation_processes"] = number_of_translation_processes
        parameter_dictionary["system_concentration_array"] = system_concentration_array
        parameter_dictionary["biophysical_parameters_dictionary"] = biophysical_parameters_dictionary
        parameter_dictionary["model_parameter_array"] = model_parameter_array
        parameter_dictionary["model_parameter_symbol_index_map"] = model_parameter_symbol_index_map
        parameter_dictionary["inverse_model_parameter_symbol_index_map"] = inverse_model_parameter_symbol_index_map
        parameter_dictionary["stoichiometric_matrix"] = SM
        parameter_dictionary["dilution_degradation_matrix"] = DM
        parameter_dictionary["specific_growth_rate"] = μ
        parameter_dictionary["half_life_translation_capacity"] = half_life_translation_capacity

        # return -
        return parameter_dictionary
        # ========================================================================================= #
    catch error
        throw(error)
    end
end


text:
"Abstract type for all flash types"
abstract type AbstractFlash end
"Abstract type for all flash types that use Newton in some form"
abstract type AbstractNewtonFlash <: AbstractFlash end

"""
Flash method that uses successive subtition.

Unconditionally convergent, does not require derivatives, but is very slow around critical regions.

See also: [`flash_2ph!`](@ref), [`NewtonFlash`](@ref) [`SSINewtonFlash`](@ref)
"""
struct SSIFlash <: AbstractFlash end

"""
    NewtonFlash([dMax = 0.2])

Flash using Newton's method for zero solve.

Only conditionally convergent, but has better convergence rate than SSI.

# Arguments
- `dMax`: dampening factor for the newton iteration

See also: [`flash_2ph!`](@ref), [`SSIFlash`](@ref) [`SSINewtonFlash`](@ref)
"""
@Base.kwdef struct NewtonFlash <: AbstractNewtonFlash
    dMax::Float64 = 0.2
end

"""
    SSINewtonFlash([swap_iter = 5,dMax = 0.2])

Perform a number of SSI iterations, followed by Newton until convergence.

See also: [`flash_2ph!`](@ref), [`SSIFlash`](@ref) [`NewtonFlash`](@ref)
"""
@Base.kwdef struct SSINewtonFlash <: AbstractNewtonFlash
    swap_iter::Int = 5
    dMax::Float64 = 0.2
end

"""
    flash_2ph(eos, c, [K], [V]; <keyword arguments>)

Perform two-phase flash with a given EOS under a set of specific conditions. Returns vapor fraction. Modifies K in-place.

Given a mixture with pressure, temperature and mole fractions, this routine performs a vapor-liquid
equilibrium calculation after a stability test.

Two outcomes are possible:
1. A single-phase condition prevails (returned vapor fraction is `NaN`) and a single phase (liquid or vapor) is stable.
2. A two-phase condition is possible. The routine produces K-values and vapor fraction so that the following holds:
    1. Isofugacity constraint for all components (``f_{li} = f_{vi}``)
    2. Molar balance for all components (``(1-V) x_i - V y_i - z_i``)
    3. Unity condition (``\\sum_i (x_i - y_i) = 0 ``)

# Arguments
- `eos`: the equation-of-state to be used for the flash
- `c`: conditions to flash the mixture at on the form `(p = 10e5, T = 303.15, z = [0.5, 0.3, 0.2])`
- `K`: optionally a buffer of length `number_of_components(eos)` used to hold K-values. Modified in-place.
- `V`: optionally the initial guess for V. If this value is not `NaN`, the stability check will be skipped.

# Keyword arguments
- `method = SSIFlash()`: Flash method to use. Can be `SSIFlash()`, `NewtonFlash()` or `SSINewtonFlash()`.
- `tolerance = 1e-8`: Tolerance for the convergence criterion. Relative to ``\\|1-R_i\\|_\\infty`` where ``R_i = f_{il}/f_{iv}``
- `maxiter = 10000`: Maximum nubmer of iterations for both stability tests and the main flash.
- `verbose = false`: Emit extra information during solve.
- `extra_out = false`: Return `(V, K, conv)` where `conv` contains iterations and oncergence status instead of just `V`.

See also: [`flash_2ph!`](@ref), [`single_phase_label`](@ref)
"""
function flash_2ph(eos, c::T, K = initial_guess_K(eos, c); method = SSIFlash(), kwarg...) where T
    nc = number_of_components(eos)
    @assert hasfield(T, :p)
    @assert hasfield(T, :T)
    @assert hasfield(T, :z)
    @assert length(c.z) == nc
    @assert length(K) == nc
    method::AbstractFlash

    S = flash_storage(eos, c, method = method)
    flash_2ph!(S, K, eos, c, update_forces = false; method = method, kwarg...)
end

"""
    flash_2ph!(storage, K, eos, c, [V]; <keyword arguments>)

Non-allocating version of [`flash_2ph`](@ref) where storage is pre-allocated.

Useful if you are performing many flashes of the same system with varying conditions.

# Arguments
- `storage`: Should be output from `flash_storage(eos, c, method = method)`. Preallocated storage.
Remaining arguments documented in [`flash_2ph`](@ref).

# Keyword arguments
- `update_forces = true`: Update the `p, T` dependent forces in storage initially.

"""
function flash_2ph!(storage, K, eos, c, V = NaN;
                                method = SSIFlash(), verbose = false, maxiter = 10000,
                                tolerance = 1e-8, extra_out = false, update_forces = true, kwarg...)
    z_min = MINIMUM_COMPOSITION
    z = c.z
    for i in eachindex(z)
        @inbounds z[i] = max(z[i], z_min)
    end
    @assert all(isfinite, K)
    forces = storage.forces
    if update_forces
        force_coefficients!(forces, eos, c)
    end
    single_phase_init = isnan(V) || V == 1.0 || V == 0.0
    if single_phase_init
        stable = stability_2ph!(storage, K, eos, c, maxiter = maxiter, verbose = verbose; kwarg...)
    else
        stable = false
    end
    converged = false
    if stable
        i = 0
    else
        i = 1
        V = solve_rachford_rice(K, z, V)
        while !converged
            V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i)
            converged = ϵ ≤ tolerance
            if converged
                if verbose
                    @info "Flash done in $i iterations." V K
                end
                break
            end
            if i == maxiter
                @warn "Flash did not converge in $i iterations. Final ϵ = $ϵ > $tolerance = tolerance"
            end
            i += 1
        end
    end
    if extra_out
        return (V, K, (its = i, converged = converged))
    else
        return V
    end
end

"""
    flash_storage(eos, [c]; <keyword arguments>)

Pre-allocate storage for `flash_2ph!`.

# Arguments
- `eos`: the equation-of-state to be used for the flash
- `c`: conditions for the mixture. Types in conditions must match later usage (e.g. for use with ForwardDiff).

# Keyword arguments
- `method = SSIFlash()`: Flash method to use. Can be `SSIFlash()`, `NewtonFlash()` or `SSINewtonFlash()`.
- `static_size = false`: Use `SArrays` and `MArrays` for fast flash, but slower compile times.
- `inc_jac`: Allocate storage for Newton/Jacobian. Required for Newton (and defaults to `true` for that method) or for `diff_externals`.
- `diff_externals = false`: Allocate storage for matrix inversion required to produce partial derivatives of flash using `set_partials`.

See also: [`flash_2ph!`](@ref) [`set_partials`](@ref)
"""
function flash_storage(eos, cond = (p = 10e5, T = 273.15, z = zeros(number_of_components(eos)));  method = SSIFlash(), kwarg...)
    out = Dict{Symbol,Any}()
    d = flash_storage_internal!(out, eos, cond, method; kwarg...)
    # Convert to named tuple
    return NamedTuple(pairs(d))
end

function flash_storage_internal!(out, eos, cond, method; inc_jac = isa(method, AbstractNewtonFlash), static_size = false, kwarg...)
    n = number_of_components(eos)
    out[:forces] = force_coefficients(eos, cond, static_size = static_size)
    if static_size
        alloc_vec = () -> @MVector zeros(n)
    else
        alloc_vec = () -> zeros(n)
    end
    out[:x] = alloc_vec()
    out[:y] = alloc_vec()

    out[:buffer1] = alloc_vec()
    out[:buffer2] = alloc_vec()
    if inc_jac
        flash_storage_internal_newton!(out, eos, cond, method, static_size = static_size; kwarg...)
    end
    return out
end

function flash_storage_internal_newton!(out, eos, cond, method; static_size = false, diff_externals = false, kwarg...)
    n = number_of_components(eos)
    np = 2*n + 1
    primary_ad(ix) = get_ad(0.0, np, :Flash, ix)
    V_ad = primary_ad(np)
    T = typeof(V_ad)
    if static_size
        x_ad = @MVector zeros(T, n)
        y_ad = @MVector zeros(T, n)
        r = @MVector zeros(np)
        J = @MMatrix zeros(np, np)
    else
        x_ad = zeros(T, n)
        y_ad = zeros(T, n)
        r = zeros(np)
        J = zeros(np, np)
    end
    out[:r] = r
    out[:J] = J

    for i = 1:n
        x_ad[i] = primary_ad(i)
        y_ad[i] = primary_ad(i+n)
    end
    out[:AD] = (x = x_ad, y = y_ad, V = V_ad)
    if diff_externals
        flash_storage_internal_inverse!(out, eos, cond, method, static_size = static_size; kwarg...)
    end
    return out
end

function flash_storage_internal_inverse!(out, eos, cond, method; static_size = false, npartials = nothing)
    n = number_of_components(eos)
    np = length(out[:r])
    external_partials = n + 2 # p, T, z_1, ... z_n
    secondary_ad(ix) = get_ad(0.0, external_partials, :InverseFlash, ix)
    p_ad = secondary_ad(1)
    T_ad = secondary_ad(2)
    T_cond = typeof(p_ad)

    if static_size
        z_ad = @MVector zeros(T_cond, n)
        J_inv = @MMatrix zeros(np, external_partials)
    else
        z_ad = zeros(T_cond, n)
        J_inv = zeros(np, external_partials)
    end
    out[:J_inv] = J_inv
    for i = 1:n
        z_ad[i] = secondary_ad(i+2)
    end
    out[:AD_cond] = (p = p_ad, T = T_ad, z = z_ad)
    if !isnothing(npartials)
        out[:buf_inv] = zeros(npartials)
    end
end

function get_ad(v::T, npartials, tag, diag_pos = nothing) where {T<:Real}
    D = ntuple(x -> T.(x == diag_pos), npartials)
    partials = ForwardDiff.Partials{npartials, T}(D)
    v = ForwardDiff.Dual{tag, T, npartials}(v, partials)
    return v
end

update_value(v::Real, newv::Real) = newv
function update_value(v::T, newv::Real) where T<:ForwardDiff.Dual
    P = v.partials
    return T(newv, P)
end

function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration)
    z = cond.z
    x, y = storage.x, storage.y
    p, T = cond.p, cond.T
    return ssi!(K, p, T, x, y, z, V, eos, forces)
end

function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
    # Initialize conditions for vapor and liquid phases based on K-values
    x = liquid_mole_fraction!(x, z, K, V)
    y = vapor_mole_fraction!(y, x, K)

    liquid = (p = p, T = T, z = x)
    vapor = (p = p, T = T, z = y)

    Z_l, s_l = prep(eos, liquid, forces)
    Z_v, s_v = prep(eos, vapor, forces)

    ϵ = zero(F)
    @inbounds for c in eachindex(K)
        f_l = component_fugacity(eos, liquid, c, Z_l, forces, s_l)
        f_v = component_fugacity(eos, vapor, c, Z_v, forces, s_v)
        r = f_l/f_v
        K[c] *= r
        ϵ = max(ϵ, abs(1-r))
    end
    V = solve_rachford_rice(K, z, V)
    return (V, ϵ)::Tuple{F, F}
end

cap_z(z) = min(max(z, MINIMUM_COMPOSITION), one(z))
cap_unit(v) = min(max(v, zero(z)), one(z))
cap_VL(v) = min(max(v, MINIMUM_COMPOSITION), 1 - MINIMUM_COMPOSITION)

function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration)
    x, y = storage.x, storage.y
    z = cond.z
    x = liquid_mole_fraction!(x, z, K, V)
    y = vapor_mole_fraction!(y, x, K)
    # Newton part
    Δ = update_and_solve!(storage, eos, cond, forces, x, y, V)
    newton_dampen!(type.dMax, Δ)
    V, ϵ = update_newton_from_increment!(K, x, y, V, Δ)
    return (V, ϵ)
end

function update_newton_from_increment!(K, x, y, V, Δ)
    ϵ = zero(eltype(K))
    n = length(x)
    @inbounds for i in 1:n
        x_n = cap_z(x[i] - Δ[i])
        y_n = cap_z(y[i] - Δ[i+n])
        K_next = y_n/x_n
        r = K_next/K[i]
        ϵ = max(ϵ, abs(1-r))
        # Assign new value, overwriting old
        K[i] = K_next
    end
    V = cap_VL(V - Δ[end])
    return (V, ϵ)
end

function get_newton_ad_values(AD, x, y, V)
    xAD = AD.x
    yAD = AD.y
    @inbounds for i in eachindex(x)
        xAD[i] = update_value(xAD[i], x[i])
        yAD[i] = update_value(yAD[i], y[i])
    end
    VAD = update_value(AD.V, V)
    return (xAD, yAD, VAD)
end

function update_primary_jacobian!(storage, eos, cond, forces, x, y, V)
    J = storage.J
    r = storage.r
    # Transfer to AD buffers
    xAD, yAD, VAD = get_newton_ad_values(storage.AD, x, y, V)
    update_flash_jacobian!(J, r, eos, cond.p, cond.T, cond.z, xAD, yAD, VAD, forces)
    return (J, r)
end

function update_secondary_jacobian!(storage, eos, c, x, y, V)
    c_ad = storage.AD_cond
    J = storage.J_inv
    p, T, z = get_inverse_ad_values(c_ad, c)
    forces = force_coefficients(eos, (p = p, T = T, z = z))
    update_flash_jacobian!(J, nothing, eos, p, T, z, x, y, V, forces)
    return (J, storage.r)
end

function update_and_solve!(storage, eos, cond, forces, x, y, V)
    J, r = update_primary_jacobian!(storage, eos, cond, forces, x, y, V)
    return solve_jacobian!(J, r)
end

function solve_jacobian!(J, r)
    # Factorize in-place, overwriting Jacobian with its factorization
    F = lu!(J)
    ldiv!(F, r)
end

function solve_jacobian!(J::MMatrix, r)
    M = length(r)
    tmp = SMatrix{M, M}(J)\SVector{M}(r)
    r .= tmp
end

function get_inverse_ad_values(cond_ad, cond)
    p, T, z = cond.p, cond.T, cond.z
    z_ad = cond_ad.z
    for i in eachindex(z)
        @inbounds z_ad[i] = update_value(z_ad[i], z[i])
    end
    p_ad = update_value(cond_ad.p, p)
    T_ad = update_value(cond_ad.T, T)
    return (p_ad, T_ad, z_ad)
end

newton_dampen!(::Nothing, Δ) = nothing

function newton_dampen!(dMax, Δ)
    ω = minimum((x) -> min(abs(dMax/x), 1.0), Δ)
    @. Δ *= ω
end

function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration)
    if iteration >= type.swap_iter
        flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration)
    else
        flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration)
    end
end

function prep(eos, cond, forces)
    s = force_scalars(eos, cond, forces)
    Z = mixture_compressibility_factor(eos, cond, forces, s)
    return (Z, s)
end

function update_flash_jacobian!(J, r, eos, p, T, z, x, y, V, forces)
    has_r = !isnothing(r)
    n = number_of_components(eos)
    liquid = (p = p, T = T, z = x)
    vapor = (p = p, T = T, z = y)

    Z_l, s_l = prep(eos, liquid, forces)
    Z_v, s_v = prep(eos, vapor, forces)

    if isa(V, ForwardDiff.Dual)
        np = length(V.partials)
    else
        np = length(p.partials)
    end
    # Isofugacity constraint
    # f_li - f_vi ∀ i
    @inbounds for c in 1:n
        f_l = component_fugacity(eos, liquid, c, Z_l, forces, s_l)
        f_v = component_fugacity(eos, vapor, c, Z_v, forces, s_v)
        Δf = f_l - f_v
        if has_r
            r[c+n] = Δf.value
        end
        for i = 1:np
            J[c+n, i] = Δf.partials[i]
        end
    end
    # x_i*(1-V) - V*y_i - z_i = 0 ∀ i
    # Σ_i x_i - y_i = 0
    T = Base.promote_type(eltype(x), eltype(y))
    L = 1 - V
    Σxy = zero(T)
    @inbounds for c in 1:n
        xc, yc = x[c], y[c]
        Σxy += (xc - yc)
        M = L*xc + V*yc - z[c]
        if has_r
            r[c] = M.value
        end
        @inbounds for i = 1:np
            J[c, i] = ∂(M, i)
        end
    end
    if has_r
        @inbounds r[end] = Σxy.value
    end
    if isa(Σxy, ForwardDiff.Dual)
        @inbounds for i = 1:np
            J[end, i] = ∂(Σxy, i)
        end
    else
        @. J[end, :] = 0
    end
end

∂(D, i) = D.partials[i]

"Compute liquid mole fraction from overall mole fraction, K-value and V vapor fraction"
@inline liquid_mole_fraction(z, K, V) = z/(1 - V + V*K)
"Compute vapor mole fraction from overall mole fraction, K-value and V vapor fraction"
@inline vapor_mole_fraction(z, K, V) = K*liquid_mole_fraction(z, K, V)
"Compute vapor mole fraction from liquid mole fraction and K-value"
@inline vapor_mole_fraction(x, K) = x*K

liquid_mole_fraction!(x, z, K, V) = begin x .= z ./ (1 .- V .+ V .* K);x end
vapor_mole_fraction!(y, x, K) = begin y .= x .* K;y end
function vapor_mole_fraction!(y,z, K, V)
    x = y
    liquid_mole_fraction!(x, z, K, V)
    vapor_mole_fraction!(y, x, K)
end

ext: jl


text:
# Code to test basic functionality and plotting
using RobustControl, PyPlot

A = [0 1;-1 -1];
B = reshape([0;1],2,1);
C = RobustControl.eye(2);
D = reshape([0;0],2,1);;
s = RobustControl.StateSpace(A,B,C);
rr,ii,om,d = RobustControl.rifd(s);
z = RobustControl.TransmissionZeros(s,2,1);
p = RobustControl.Poles(s);
d = RobustControl.Damping(s);
om = RobustControl.NaturalFrequencies(s);

# Frequency response with auto frequency grid.
G,om = RobustControl.FrequencyResponse(s,1,1);

# Frequency response with specified frequency grid.
om1 = RobustControl.logspace(1E-2,1E2,200);
G1 = RobustControl.FrequencyResponse(s,2,1,om1);

# Test plots
figure(1);clf(); pygui(true);
subplot(2,1,1);
RobustControl.BodeMagPlot(G,om); title("Bode Mag Plot");
RobustControl.BodeMagPlot(G1,om1);
grid(which="both",color=[1,1,1]*0.9,linewidth=0.35)
legend(["G1","G2"]);
tight_layout();

# Another way to get Bode plots -- just magnitude.
subplot(2,1,2);
RobustControl.BodeMagPlot(s,1,1);
RobustControl.BodeMagPlot(s,2,1);
grid(which="both",color=[1,1,1]*0.9,linewidth=0.35)
legend(["G1","G2"]);
tight_layout();

# 3rd way -- get both mag and phase in one plot
figure(2); clf();
RobustControl.BodePlot(s,1,1);
RobustControl.BodePlot(s,2,1);
subplot(211);grid(which="both",color=[1,1,1]*0.9,linewidth=0.35)
subplot(212);grid(which="both",color=[1,1,1]*0.9,linewidth=0.35)
subplot(211);legend(["G1","G2"]);
title("Bode Plot");
tight_layout();

# Nyquist Plot
figure(3);clf();
subplot(2,1,1);RobustControl.NyquistPlot(s,1,1); title("Nyquist Plot"); # With auto frequency grid
subplot(2,1,2);RobustControl.NyquistPlot(s,1,1,om1); # With specified frequency grid

# System Responses
figure(4); clf();
subplot(2,2,1);
t,y = RobustControl.ImpulseResponse(s,1,1);
plot(t,y); grid();
xlabel("Time (sec)");
#ylabel("Impulse Response");
title("Impulse Response");tight_layout();

subplot(2,2,2);
t,y = RobustControl.StepResponse(s,1,1);
plot(t,y); grid();
xlabel("Time (sec)");
#ylabel("Step Response");
title("Step Response");tight_layout();

subplot(2,2,3);
t = Vector(0:0.01:50);
u = ones(length(t)) + 0.1*randn(Float64,length(t)); # Step + Disturbance
x0 = [0.0;0.0];
y = RobustControl.ForcedResponse(s,1,1,t,u,x0);
plot(t,y); grid();
xlabel("Time (sec)");
#ylabel("Time Response to Step+Disturbance");
title("Forced Response");tight_layout();

subplot(2,2,4);
t = Vector(0:0.01:10);
x0 = [1.0;0.0];
t1,y1 = RobustControl.InitialResponse(s,x0); # Auto time grid
y2 = RobustControl.InitialResponse(s,x0,t); # Specified time grid
l1=plot(t1,y1,"b",linewidth=1);
l2=plot(t,y2,"r",linestyle=":",linewidth=2);
xlabel("Time (sec)");
#ylabel("Initial Response");
legend([l1[1],l2[1]],["(auto time grid)","(specified time grid)"]);
title("Initial Response");
grid(); tight_layout();

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
export print_sudoku, lines, fill
#Modified from: https://github.com/scheinerman/Sudoku.jl

""""
    print_sudoku(grid::Matrix)

Print a sudoku grid with a nice format
"""
function print_sudoku(grid::Matrix) ## does this function used in the script? if not maybe remove it to make it look tidy :)
    for string in lines(grid)
        println(string)
    end
    nothing
end

"""
    lines(grid::Matrix)

Draw the edges of a Sudoku grid
"""
function lines(grid::Matrix) ## does this function used in the script? if not maybe remove it to make it look tidy :)
    horiz_piece = "---------+"
    sep_line = "+" * (horiz_piece)^3

    result = Vector{String}()
    push!(result, sep_line)

    for i=1:9
        push!(result, "|" *fill(grid[i,1:3])*"|" * fill(grid[i,4:6])*"|" * fill(grid[i, 7:9])*"|" )
        if mod(i,3) == 0
            push!(result,sep_line)
        end
    end
    return result
end

"""
    fill(r::Vector)

Fills the sudoku grid, using " - " where zero is found,
takes as input a row vector
"""
function fill(r::Vector) ## does this function used in the script? if not maybe remove it to make it look tidy :)
    result = ""
    for k=1:3
        a = r[k]
        result *=  a > 0 ? " $a " : " - "
    end
    result *= ""
    return result
end


text:
"""
	This File is used to calculate the Heuristics.
    Also the init results
"""

include("../src/Graph.jl")
include("../src/Tools.jl")
include("../src/Algorithm.jl")
include("../src/Sampling.jl")
include("../src/Methods.jl")

using StatsBase
using JSON
using JLD2

"""
Set Seed Nodes & change innate opinions
Set parameters.
k = 50
epsilon : s + epsilon
delta : edge probability parameter
beta: threshold
json struct js
"""

epsilon = 0.1
delta = 0.5
beta = 0.5
tags = []

jsdict = Dict{String, Any}()

# Read Graph
buf = split(ARGS[1], ',')
fileName = string("../data/all/", buf[1], ".jld2")

(name, n, m, V, E) = load(fileName, "name", "n", "m", "V", "E")
G = Graph(n, m, V, E)

merge!(jsdict, Dict("name"=>name, "n"=>n, "m"=>m))

if ARGS[2] == "1"
    #Weighted Cascade
    push!(tags, "WC")
    E_p = load(fileName, "E_WC")
elseif ARGS[2] == "2"
    #Trivalency
    push!(tags, "TM") # Trivalency Model
    E_p = load(fileName, "E_TM")
end
G2 = Graph(n, m, V, E_p)

# Before Making changes
if ARGS[3] == "1"
    # Uniform distribution
    push!(tags, "Uniform")
    s = load(fileName, "s_Uni")
elseif ARGS[3] == "2"
    # Exponential distribution
    push!(tags, "Exponential")
    s = load(fileName, "s_Exp")
elseif ARGS[3] == "3"
    # Power-law distribution
    push!(tags, "Power-law")
    s = load(fileName, "s_Pow")
end

# Different parameters to select k
if ARGS[4] == "1"
    push!(tags, "0.005")
    k = load(fileName, "k1")
elseif ARGS[4] == "2"
    push!(tags, "0.015")
    k = load(fileName, "k2")
elseif ARGS[4] == "3"
    push!(tags, "0.02")
    k = load(fileName, "k3")
end

merge!(jsdict, Dict("k"=>k))
merge!(jsdict, Dict("model"=>tags[1], "s_dist"=>tags[2], "k_par"=>tags[3]))

name = ["z_sum", "aci", "ad", "ap", "aidc"]

##################### Above are configs
#####################
# The only difference for marketing setting and backfire setting is that Δs is different

if ARGS[5] == "1"
    Δs = delta_s_marketing(s, epsilon)
elseif ARGS[5] == "2"
    Δs = delta_s_backfire(s, epsilon, beta)
end

##################### Below are related to Heuristics
######################

#println("s is ", Δs)
Δaci, Δad, Δap, Δaidc, Δaidc_upp = assign_node_weights(G, s, Δs)
#println("aci is ", Δaci, "ad is ", Δad, "ap is ", Δap, "aidc is ", Δaidc)

lin = [Δs, Δaci, Δad, Δap, Δaidc, Δaidc_upp]
lin_labels = ["max_z_", "max_aci_", "max_ad_", "max_ap_", "max_aidc_", "max_aidc_upp_"]

ris_me, p_me, T_me, R_me = [], [], [], []

# The initial of different measures
_, init_aci, init_ad, init_ap, init_aidc = Approx(G, s)
init_result = [sum(s), init_aci, init_ad, init_ap, init_aidc, init_aidc] # the last one is a upper bound

init_labels = "init_".*name
for i in 1:length(name)
    merge!(jsdict, Dict(init_labels[i]=>init_result[i]))
end

#println("init aci is ", init_aci)

# Calculation
for i in 1:length(lin)

    # if i <= 5
    #     println("now is ", name[i])
    # else
    #     println("now is aidc_upp")
    # end

    flag = 1

    #if (i == 2) || (i == 3) || (i == 4)
    #    flag = 1
    #end

    (ris, seeds_lin, R_size, T_lin) = ris_measure_linear(G2, k, lin[i], flag, delta) ## RR to select seed nodes
    #if i == 2
    #    println("aci gain linear is ", ris)
    #end
    (_, p) = cascadeIC(G2, seeds_lin, round(Int, 2*(10*10)*log(2*n)), delta) ## direct calculation, just to check
    push!(ris_me, ris)
    #println("ris_lower_bound is: ", ris)
    push!(p_me, p)
    push!(T_me, T_lin)
    push!(R_me, R_size)

end

# Only in order to check the z_sum
#println("z_sum is, if we calculated by delta s ", sum(Δs .* p_me[1]))

# Output
for i in 1:length(lin)
    nn = lin_labels[i].*name
    lin_nn = lin_labels[i]*"lin" # only the linear part
    lin_t = lin_labels[i]*"time"
    lin_size = lin_labels[i]*"r_size"

    s2 = Δs .* p_me[i] .+ s

    _, aci, ad, ap, aidc = Approx(G, s2)
    z_sum = sum(s2)

    #if lin_labels[i] == "max_aci_"
    #    println("seed aci gain is, ", aci - init_aci)
    #end

    #println("z_sum gain is ", z_sum - sum(s))

    result = [z_sum, aci, ad, ap, aidc]
    merge!(jsdict, Dict(nn[1]=>result[1], nn[2]=>result[2], nn[3]=>result[3], nn[4]=>result[4], nn[5]=>result[5]))
    merge!(jsdict, Dict(lin_t=>T_me[i]))
    merge!(jsdict, Dict(lin_size=>R_me[i]))
    merge!(jsdict, Dict(lin_nn=>ris_me[i] + init_result[i]))
end

#println("sum of s is ", sum(s))

######################
####################### IO

js = JSON.json(jsdict)

println(jsdict)

if ARGS[5] == "1"
    open("out/method_marketing_small_par.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
elseif ARGS[5] == "2"
    open("out/method_backfire_small_par.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
end


text:
macro define_show(type)
  return quote
      Base.show(io::IO, x::$(esc(type))) = print(io, pystr(x.o))
      Base.show(io::IO, ::MIME"text/plain", x::$(esc(type))) = print(io, pystr(x.o))
      Base.show(io::IO, ::MIME"text/latex", x::$(esc(type))) = print(io, "\\begin{align*}" * galgebra.printer.latex(x.o) * "\\end{align*}")
  end
end

macro delegate_properties(type, obj_field)
  return quote
      Base.convert(::Type{$(esc(type))}, o::PyCall.PyObject) = $(esc(type))(o)
      PyCall.PyObject(o::$(esc(type))) = PyCall.PyObject(getfield(o, $obj_field))

      function Base.getproperty(o::$(esc(type)), s::AbstractString)
          if s == String($obj_field)
              return getfield(o, $obj_field)
          else
              return getproperty(getfield(o, $obj_field), s)
          end
      end

      Base.getproperty(o::$(esc(type)), s::Symbol) = getproperty(o, String(s))

      Base.propertynames(o::$(esc(type))) = map(x->Symbol(first(x)),
                                      pycall(inspect."getmembers", PyObject, getfield(o, $obj_field)))

      # avoiding method ambiguity
      Base.setproperty!(o::$(esc(type)), s::Symbol, v) = _setproperty!(o,s,v)
      Base.setproperty!(o::$(esc(type)), s::AbstractString, v) = _setproperty!(o,s,v)

      function _setproperty!(o::$(esc(type)), s::Union{Symbol,AbstractString}, v)
          obj = getfield(o, $obj_field)
          setproperty!(obj, s, v)
          o
      end

      hasproperty(o::$(esc(type)), s::Symbol) = hasproperty(getfield(o, $obj_field), s)
      hasproperty(o::$(esc(type)), s::AbstractString) = hasproperty(getfield(o, $obj_field), s)
  end
end

macro delegate_doc(type)
  return quote
      # Expose Python docstrings to the Julia doc system
      Docs.getdoc(x::$(esc(type))) = Text(convert(String, x."__doc__"))
      Docs.Binding(x::$(esc(type)), s::Symbol) = getproperty(x, s)
  end
end

macro define_op(type, op, method)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(type)), y::$(esc(type))) = x.$method(y)
  end
end

macro define_op_with_impl(type, op, impl)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(type)), y::$(esc(type))) = $impl
  end
end

macro define_lop(type, rtype, op, lmethod)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(type)), y::$(esc(rtype))) = x.$lmethod(y)
  end
end

macro define_rop(type, ltype, op, rmethod)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(ltype)), y::$(esc(type))) = y.$rmethod(x)
  end
end

macro define_unary_op(type, op, method)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(type))) = x.$method()
  end
end

macro define_unary_op_with_impl(type, op, impl)
  return quote
    Core.@__doc__ @pure $(esc(op))(x::$(esc(type))) = $impl
  end
end

macro define_postfix_symbol(super_script)
  return quote
      const $(esc(super_script)) = () -> String(:super_script)
      export $(esc(super_script))
  end
end

macro define_postfix_op(type, super_script, func)
  return quote
      Base.:(*)(x::$(esc(type)),::typeof($(esc(super_script)))) = $func(x)
      Base.:(^)(x::$(esc(type)),::typeof($(esc(super_script)))) = $func(x)
  end
end


text:
using LatinSquares, Test

A,B = ortho_latin(4)
@test check_ortho(A,B)
A,B = ortho_latin(5)
@test check_ortho(A,B)

A,B = ortho_latin(4,true)
@test check_ortho(A,B)
@test A==B'


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)
module Foo
import Qux
this = "Foo2"
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# # Shear instability of a free-surface flow
#
# This script simulates the instability of a sheared, free-surface
# flow using `ClimateMachine.Ocean.HydrostaticBoussinesqSuperModel`.

using Printf
using Plots
using ClimateMachine

ClimateMachine.init()

ClimateMachine.Settings.array_type = Array

using ClimateMachine.Ocean
using ClimateMachine.Ocean.Domains
using ClimateMachine.Ocean.Fields

using ClimateMachine.GenericCallbacks: EveryXSimulationTime
using ClimateMachine.Ocean: current_step, Δt, current_time
using CLIMAParameters: AbstractEarthParameterSet, Planet

# We begin by specifying the domain and mesh,

domain = RectangularDomain(
    Ne = (24, 24, 1),
    Np = 4,
    x = (-3π, 3π),
    y = (-3π, 3π),
    z = (0, 1),
    periodicity = (true, false, false),
)

# Note that the default solid-wall boundary conditions are free-slip and
# insulating on tracers. Next, we specify model parameters and the sheared
# initial conditions

struct NonDimensionalParameters <: AbstractEarthParameterSet end
Planet.grav(::NonDimensionalParameters) = 1

initial_conditions = InitialConditions(
    u = (x, y, z) -> tanh(y) + 0.1 * cos(x / 3) + 0.01 * randn(),
    v = (x, y, z) -> 0.1 * sin(y / 3),
    θ = (x, y, z) -> x,
)

model = Ocean.HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = 0.05,
    initial_conditions = initial_conditions,
    parameters = NonDimensionalParameters(),
    turbulence_closure = (νʰ = 1e-2, κʰ = 1e-2, νᶻ = 1e-2, κᶻ = 1e-2),
    rusanov_wave_speeds = (cʰ = 0.1, cᶻ = 1),
    boundary_tags = ((0, 0), (1, 1), (1, 2)),
    boundary_conditions = (
        OceanBC(Impenetrable(FreeSlip()), Insulating()),
        OceanBC(Penetrable(FreeSlip()), Insulating()),
    ),
)

# We prepare a callback that periodically fetches the horizontal velocity and
# tracer concentration for later animation,

u, v, η, θ = model.fields
fetched_states = []

start_time = time_ns()

data_fetcher = EveryXSimulationTime(1) do
    step = @sprintf("Step: %d", current_step(model))
    time = @sprintf("time: %.2f min", current_time(model) / 60)
    max_u = @sprintf("max|u|: %.6f", maximum(abs, u))

    elapsed = (time_ns() - start_time) * 1e-9
    wall_time = @sprintf("elapsed wall time: %.2f min", elapsed / 60)

    @info "$step, $time, $max_u, $wall_time"

    push!(
        fetched_states,
        (u = assemble(u), θ = assemble(θ), time = current_time(model)),
    )
end

# and then run the simulation.

model.solver_configuration.timeend = 100.0

result = ClimateMachine.invoke!(
    model.solver_configuration;
    user_callbacks = [data_fetcher],
)

# Finally, we make an animation of the evolving shear instability.

animation = @animate for (i, state) in enumerate(fetched_states)
    local u
    local θ

    @info "Plotting frame $i of $(length(fetched_states))..."

    kwargs =
        (xlim = domain.x, ylim = domain.y, linewidth = 0, aspectratio = 1)

    x, y = state.u.x[:, 1, 1], state.u.y[1, :, 1]

    u = state.u.data[:, :, 1]
    θ = state.θ.data[:, :, 1]

    ulim = 1
    θlim = 8

    ulevels = range(-ulim, ulim, length = 31)
    θlevels = range(-θlim, θlim, length = 31)

    u_plot = contourf(
        x,
        y,
        clamp.(u, -ulim, ulim)';
        levels = ulevels,
        color = :balance,
        kwargs...,
    )
    θ_plot = contourf(
        x,
        y,
        clamp.(θ, -θlim, θlim)';
        levels = θlevels,
        color = :thermal,
        kwargs...,
    )

    u_title = @sprintf("u at t = %.2f", state.time)
    θ_title = @sprintf("θ at t = %.2f", state.time)

    plot(u_plot, θ_plot, title = [u_title θ_title], size = (1200, 500))
end

gif(animation, "shear_instability.gif", fps = 8)

max_forks_repo_name: jm-c/CLIMA


text:
# Expected output:
#  calling callback with: Base1
#  callback with: Base1
#  calling callback with: Base2
#  callback with: Base2

using Clang.cindex, Clang.wrap_c, Clang.wrap_cpp

hpath = "test_classparents.h"
clsname = "Derived"

  idx = cindex.idx_create(0,1)
  tu = cindex.tu_parse(idx, hpath,
    ["-x", "c++",
     "-I./",
     "-I$JULIA_HOME/../../deps/llvm-3.2/build/Release/lib/clang/3.3/include",
     "-v", "-c"], 0)

  topcu = cindex.getTranslationUnitCursor(tu)
  topcl = children(topcu)

  clscu = cindex.CXCursor()
  found = false
  for i=1:topcl.size
    clscu = topcl[i]
    if(cu_kind(clscu) == cindex.CurKind.CLASSDECL && name(clscu) == clsname)
      found = true
      println("found at idx: $i")
      break
    end
  end

function cxxclass_cb(cursor::Ptr{Uint8}, data::Ptr{Void})
         cu = CXCursor()
         ccall(:memcpy, Void, (Ptr{Void},Ptr{Void}, Uint), cu.data, cursor, cindex.CXCursor_size)
         println("callback with: ", name(cu))
         return int(0)
       end

_cb = cfunction(cxxclass_cb, Int, (Ptr{Uint8},Ptr{Void}))

ccall( (:wci_getCXXClassParents, :libwrapclang), Int, (Ptr{Uint8}, Ptr{Void}), clscu.data, _cb)


text:
export
    InstantModel

"""
    abstract type InstantModel{I,O} <: Model{I,Nothing,O}

A model for a variable with no time dependencies.
Since this model has no transitions, it can only called with
`make_initial` and not `make_transition`
(i.e. `make_transition` = `make_initial`)
"""
abstract type InstantModel{I,O} <: Model{I,Nothing,O} end

# Instant models can be freely used dynamically, with no time dependencies
make_transition(m::InstantModel, parenttimes, time) = make_initial(m)


text:
# ---
# title: JuliaRL\_BasicDQN\_EmptyRoom
# cover: assets/JuliaRL_BasicDQN_EmptyRoom.png
# description: A simple example to demonstrate how to use environments in GridWorlds.jl
# date: 2021-05-22
# author: "[Siddharth Bhatia](https://github.com/Sid-Bhatia-0)"
# ---

#+ tangle=true
using ReinforcementLearning
using GridWorlds
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BasicDQN},
    ::Val{:EmptyRoom},
    ::Nothing;
    seed=123,
)
    rng = StableRNG(seed)

    inner_env = GridWorlds.EmptyRoomDirected(rng=rng)
    action_space_mapping = x -> Base.OneTo(length(RLBase.action_space(inner_env)))
    action_mapping = i -> RLBase.action_space(inner_env)[i]
    env = RLEnvs.ActionTransformedEnv(
        inner_env,
        action_space_mapping=action_space_mapping,
        action_mapping=action_mapping,
    )
    env = RLEnvs.StateTransformedEnv(env;state_mapping=x -> vec(Float32.(x)))
    env = RewardOverriddenEnv(env, x -> x - convert(typeof(x), 0.01))
    env = MaxTimeoutEnv(env, 240)

    ns, na = length(state(env)), length(action_space(env))
    agent = Agent(
        policy=QBasedPolicy(
            learner=BasicDQNLearner(
                approximator=NeuralNetworkApproximator(
                    model=Chain(
                        Dense(ns, 128, relu; init=glorot_uniform(rng)),
                        Dense(128, 128, relu; init=glorot_uniform(rng)),
                        Dense(128, na; init=glorot_uniform(rng)),
                    ) |> cpu,
                    optimizer=ADAM(),
                ),
                batch_size=32,
                min_replay_history=100,
                loss_func=huber_loss,
                rng=rng,
            ),
            explorer=EpsilonGreedyExplorer(
                kind=:exp,
                ϵ_stable=0.01,
                decay_steps=500,
                rng=rng,
            ),
        ),
        trajectory=CircularArraySARTTrajectory(
            capacity=1000,
            state=Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=false)
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "")
end

#+ tangle=false
using Plots
pyplot() #hide
ex = E`JuliaRL_BasicDQN_EmptyRoom`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_BasicDQN_EmptyRoom.png") #hide

# ![](assets/JuliaRL_BasicDQN_EmptyRoom.png)


text:
import PlotlyJS
import Colors
import DelimitedFiles

cm = mapslices(
    r->Colors.RGB(r[1],r[2],r[3]),
    DelimitedFiles.readdlm(joinpath(dirname(pathof(CompScienceMeshes)),"..","examples","cm.csv"),',',Float64),
    # readcsv(Pkg.dir("CompScienceMeshes","examples","cm.csv")),
    dims=[2])

function patch(Γ, fcr=nothing)

    v = vertexarray(Γ)
    c = cellarray(Γ)

    x = v[:,1]; y = v[:,2]; z = v[:,3]
    i = c[:,1].-1; j = c[:,2].-1; k = c[:,3].-1

    if fcr == nothing
        #a = [barytocart(chart(Γ,cells(Γ,i)),[1,1]/3)[3] for i in 1:numcells(Γ)]
        a = [cartesian(center(chart(Γ,cell)))[3] for cell in cells(Γ)]
    else
        a = fcr
    end

    m, M = extrema(a)
    if isapprox(m, M)
        n = ones(Integer, a)
    else
        n = floor.(Integer, (a.-m)/(M-m)*(length(cm)-1)).+1
    end
    fc = [cm[i] for i in n]

    s = PlotlyJS.mesh3d(;
        x=x,
        y=y,
        z=z,
        i=i,
        j=j,
        k=k,
        facecolor=fc,
        colorbar=PlotlyJS.attr(title="z"),
    )

    #PlotlyJS.plot([s])
    s
end


text:
mutable struct Conf
    MIP_SOLVER
    MIP_SUCCESS
end

"""
    conf

A struct with fields for global configuration of the `Allocations` module.

# Fields

    MIP_SOLVER :: Any

The (factory for) the JuMP optimizer to be used (by default) for mixed-integer
programming. Initially set to `Cbc.Optimizer`, with a `logLevel` of `0`. This
can be overridden either by setting `MIP_SOLVER` to another value (e.g., using
the JuMP function `optimizer_with_attributes`) or by passing the solver
directly to the appropriate allocation functions.

    MIP_SUCCESS :: Any

Container of acceptable MIP statuses. By default, has the value
`[MOI.OPTIMAL]`.
"""
const conf = Conf(nothing, nothing)

function __init__()
    conf.MIP_SOLVER = optimizer_with_attributes(
        Cbc.Optimizer, "logLevel" => 0
    )
    conf.MIP_SUCCESS = [MOI.OPTIMAL]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailabilitySetListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{AvailabilitySet} } # spec name: value

    function AvailabilitySetListResult(;value=nothing)
        o = new()
        validate_property(AvailabilitySetListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type AvailabilitySetListResult

const _property_map_AvailabilitySetListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_AvailabilitySetListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{AvailabilitySet}")
Base.propertynames(::Type{ AvailabilitySetListResult }) = collect(keys(_property_map_AvailabilitySetListResult))
Swagger.property_type(::Type{ AvailabilitySetListResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AvailabilitySetListResult[name]))}
Swagger.field_name(::Type{ AvailabilitySetListResult }, property_name::Symbol) =  _property_map_AvailabilitySetListResult[property_name]

function check_required(o::AvailabilitySetListResult)
    true
end

function validate_property(::Type{ AvailabilitySetListResult }, name::Symbol, val)
end


text:
abstract type Sampler{T} end

function set_x0!(S::Sampler{T}, x0::AbstractVector{T}) where T
    S.x = x0
    return S
end

mutable struct HMCSampler{T,
                          LL <: LogLikelihood{T},
                          V<:AbstractVector{T},
                          S<:AbstractMatrix{T},
                          SS<:AbstractMatrix{T}
                          } <: Sampler{T}
    logf::LL
    x::V
    ϵ::T
    M::Int
    ΣL::S
    Σ2::SS
end

function HMCSampler(logf::LogLikelihood{T}, x::AbstractVector{T}, ϵ, M, ΣL) where T
    Σ2 = ΣL * ΣL'
    HMCSampler{eltype(x), typeof(logf), typeof(x), typeof(ΣL), typeof(Σ2)}(logf, x, ϵ, M, ΣL, Σ2)
end

function Base.copy(S::HMCSampler)
    HMCSampler(copy(S.logf), copy(S.x), copy(S.ϵ), copy(S.M), copy(S.ΣL), copy(S.Σ2))
end

mutable struct LAHMCSampler{T,
                            LL<:LogLikelihood{T},
                            V<:AbstractVector{T},
                            S<:AbstractMatrix{T},
                            SS<:AbstractMatrix{T}} <: Sampler{T}
    logf::LL
    x::V
    ϵ::T
    M::Int
    K::Int
    ΣL::S
    Σ2::SS
    β::Float64
    p::V
    k::Vector{Int}
end

function LAHMCSampler(logf::LogLikelihood{T}, x::AbstractVector{T}, ϵ, M, K, ΣL, β) where T
    p = ΣL*randn(T, length(x))
    Σ2 = ΣL * ΣL'
    LAHMCSampler{eltype(x), typeof(logf), typeof(x), typeof(ΣL), typeof(Σ2)}(logf, x, ϵ, M, K, ΣL, Σ2, β, p, zeros(Int, K+1))
end

function Base.copy(S::LAHMCSampler)
    LAHMCSampler(copy(S.logf), copy(S.x), copy(S.ϵ), S.M, S.K, copy(S.ΣL), copy(S.Σ2), S.β, copy(S.p), copy(S.k))
end

function momentum_resample!(S::LAHMCSampler{T}; β = S.β) where T
    pp1 = (S.ΣL * randn(T, length(S.x))) .* sqrt(β)
    pp2 = S.p .* sqrt(1 - β)
    S.p = pp1 + pp2
    return S
end

function sample!(S::HMCSampler{T}) where T
    ΣL = S.ΣL
    Σ2 = S.Σ2

    ϵ = S.ϵ
    M = S.M
    logf = S.logf

    x0 = S.x
    p0 = ΣL * randn(T, length(x0))

    x1, p1, logf1, grad1, logf0, grad0 = leapfrog(x0, p0, ϵ, M, Σ2, logf)
    p1 = -p1

    current_U  = -logf0
    current_K  = KEngry(ΣL, p0)
    proposed_U = -logf1
    proposed_K = KEngry(ΣL, p1)

    prob = exp(current_U - proposed_U + current_K - proposed_K)
    if rand() < prob
        S.x = x1
    end

    return S
end

KEngry(ΣL, p) = (half = ΣL\p; 0.5*dot(half,half))
∇KEngry(Σ2, p) = Σ2\p

function leapfrog(x0::AbstractVector{T}, p0::AbstractVector{T}, ϵ, M, Σ2, logf::LogLikelihood{T}) where T
    logf0, grad0 = logf(x0, Val(:vg))
    x1 = copy(x0)
    p1 = copy(p0)
    p1 += T(0.5) * ϵ * grad0
    logf1 = T(0.0)
    grad1 = similar(grad0)
    for i in 1:M
        x1          += ϵ * ∇KEngry(Σ2, p1)
        logf1, grad1 = logf(x1, Val(:vg))
        p1          += ϵ * grad1
    end
    p1 -= T(0.5) * ϵ * grad1
    return x1, p1, logf1, grad1, logf0, grad0
end

function sample!(S::LAHMCSampler{T}) where T
    ΣL = S.ΣL
    Σ2 = ΣL * ΣL'

    x0 = S.x[:]
    x1 = S.x[:]

    p0 = S.p[:]
    p1 = S.p[:]

    K = S.K
    M = S.M
    ϵ = S.ϵ
    β = S.β
    logf = S.logf

    x_candidate = [x0 for k in 1:K+1]
    p_candidate = [p0 for k in 1:K+1]
    rand_comparison = rand()
    p_cum = 0.0
    C = zeros(Float64, K+1, K+1)
    fill!(C, NaN)

    for k in 1:K
        x_new, p_new= leapfrog( x_candidate[k], p_candidate[k],
                                ϵ, M, Σ2, logf)

        x_candidate[k+1] = x_new
        p_candidate[k+1] = p_new
        p_cum, Cl = leap_prob_recurse( x_candidate[1:(k+1)], p_candidate[1:(k+1)],
                                       C[1:(k+1),1:(k+1)], ΣL, logf)
        C[1:(k+1),1:(k+1)] = Cl

        if p_cum >= rand_comparison
            S.k[k] += 1
            x1 = x_new
            p1 = p_new
            break
        end
    end
    if p_cum < rand_comparison
        S.k[K+1] += 1
        x1 = copy(x0)
        p1 = -p0
    end

    pp1 = (ΣL * randn(T, length(x0))) .* sqrt(β)
    pp2 = p1 .* sqrt(1 - β)
    S.p = pp1 + pp2

    S.x = x1
    return S
end

function leap_prob(θ₁::AbstractVector{T}, r₁::AbstractVector{T},
                   θ₂::AbstractVector{T}, r₂::AbstractVector{T},
                   ΣL, logfgrad::LogLikelihood{T}) where T
    logf1 = logfgrad(θ₁, Val(:v))
    logf2 = logfgrad(θ₂, Val(:v))

    Ediff = (-logf1 + KEngry(ΣL, r₁))-(-logf2 + KEngry(ΣL, r₂))
    return min(1.0, exp(Ediff))
end

function leap_prob_recurse(θ_chain, r_chain, C, ΣL, logfgrad::LogLikelihood{T}) where T
    if isfinite(C[1,end])
        # we've already visited this leaf
        cumu = C[1,end]
        return cumu, C
    end

    if length(θ_chain) == 2
        p_acc = leap_prob(θ_chain[1], r_chain[1], θ_chain[2], r_chain[2], ΣL, logfgrad)
        C[1,end] = p_acc
        return p_acc, C
    end

    cum_forward, Cl = leap_prob_recurse(θ_chain[1:end-1],
                                        r_chain[1:end-1],
                                        C[1:end-1, 1:end-1], ΣL, logfgrad)
    C[1:end-1,1:end-1] = Cl
    cum_reverse, Cl = leap_prob_recurse(θ_chain[end:-1:2],
                                        r_chain[end:-1:2],
                                        C[end:-1:2, end:-1:2], ΣL, logfgrad)
    C[end:-1:2, end:-1:2] = Cl

    _H0 = logfgrad(θ_chain[1], Val(:v))
    _H1 = logfgrad(θ_chain[end], Val(:v))

    H0 = -_H0 + KEngry(ΣL, r_chain[1])
    H1 = -_H1 + KEngry(ΣL, r_chain[end])
    Ediff = H0 - H1

    start_state_ratio = exp(Ediff)

    if isnan(start_state_ratio*(1 - cum_reverse))
        prob = 1 - cum_forward
    else
        prob = min(1 - cum_forward, start_state_ratio*(1 - cum_reverse))
    end

    cumu = cum_forward + prob
    C[1, end] = cumu

    return cumu, C
end


text:
struct LeastSquares{T}
    N::Int
    fz::Vector{T}
    fη::Vector{Vector{T}}
    η::Vector{Vector{T}}
    cost::Function
    grad::Function
    hess::Function
    val_cost::Vector{T}
    val_grad::Vector{T}
    val_hess::Matrix{T}
    cache_cost::Vector{T}
    cache_grad::Vector{T}
    cache_hess::Matrix{T}
    θ::Vector{T}
    Δ::Vector{T}
    solver::RoboDojo.LinearSolver
end

function eval_cost!(ls::LeastSquares)
    fill!(ls.val_cost, 0.0)
    for i = 1:ls.N
        ls.cost(ls.cache_cost, ls.fz, ls.fη[i], ls.η[i], ls.θ)
        ls.val_cost .+= ls.cache_cost
    end
end

function eval_grad!(ls::LeastSquares)
    fill!(ls.val_grad, 0.0)
    for i = 1:ls.N
        ls.grad(ls.cache_grad, ls.fz, ls.fη[i], ls.η[i], ls.θ)
        ls.val_grad .+= ls.cache_grad
    end
end

function eval_hess!(ls::LeastSquares)
    fill!(ls.val_hess, 0.0)
    for i = 1:ls.N
        ls.hess(ls.cache_hess, ls.fz, ls.fη[i], ls.η[i], ls.θ)
        ls.val_hess .+= ls.cache_hess
    end
end

function update!(ls::LeastSquares; tol=1.0e-8, verbose=true)
    iter = 0
    eval_grad!(ls)
    res = norm(ls.val_grad, Inf)
    verbose && (@show res)
    while res > tol && iter < 100
        eval_hess!(ls)
        fill!(ls.Δ, 0.0)
        RoboDojo.linear_solve!(ls.solver, ls.Δ, ls.val_hess, ls.val_grad)
        ls.θ .-= ls.Δ
        eval_grad!(ls)
        res = norm(ls.val_grad, Inf)
        verbose && (@show res)
        verbose && (@show ls.θ)
        iter += 1
    end
end

# ## test
# function f(z)
#     x = z[1:2]
#     u = z[2 .+ (1:1)]
#     A = [1.0 1.0; 0.0 1.0]
#     B = [0.0; 1.0]

#     A * x + B * u[1]
# end

# nx = 2
# nu = 1
# nz = nx + nu
# nθ = nx * (nz)

# @variables fz[1:nx] fη[1:nx] η[1:nz] θ[1:nθ]

# function cost(fz, fη, η, θ)
#     M = reshape(θ, nx, nz)
#     r = fη - fz - M * η
#     return [transpose(r) * r]
#     return [transpose(θ) * θ]
# end

# c = cost(fz, fη, η, θ)
# cθ = Symbolics.gradient(c[1], θ)
# cθθ = Symbolics.hessian(c[1], θ)

# c_func = eval(Symbolics.build_function(c, fz, fη, η, θ)[2])
# cθ_func = eval(Symbolics.build_function(cθ, fz, fη, η, θ)[2])
# cθθ_func = eval(Symbolics.build_function(cθθ, fz, fη, η, θ)[2])

# fz0 = rand(nx)
# fη0 = rand(nx)
# η0 = rand(nz)
# θ0 = rand(nθ)
# θθ0 = rand(nθ, nθ)

# c0 = zeros(1)
# cθ0 = zeros(nθ)
# cθθ0 = zeros(nθ, nθ)

# c_func(c0, fz0, fη0, η0, θ0)
# @benchmark c_func($c0, $fz0, $fη0, $η0, $θ0)

# cθ_func(cθ0, fz0, fη0, η0, θ0)
# @benchmark cθ_func($cθ0, $fz0, $fη0, $η0, $θ0)

# cθθ_func(cθθ0, fz0, fη0, η0, θ0)
# @benchmark cθ_func($cθθ0, $fz0, $fη0, $η0, $θθ0)

# N = 2 * nz
# η = [zeros(nz) for i = 1:N]
# ϵ = 0.1
# for i = 1:nz
#     η[i][i] = ϵ
#     η[i+nz][i] = -ϵ
# end

# z0 = rand(nz)
# fz = f(z0)
# fη = [f(z0 + η[i]) for i = 1:N]

# ls = LeastSquares(N,
#     fz, fη, η,
#     c_func, cθ_func, cθθ_func,
#     [0.0], zeros(nθ), zeros(nθ, nθ),
#     [0.0], zeros(nθ), zeros(nθ, nθ),
#     zeros(nθ), zeros(nθ),
#     lu_solver(zeros(nθ, nθ)))

# # eval_cost!(ls, z0)
# # @benchmark eval_cost!($ls)

# # eval_grad!(ls, z0)
# # @benchmark eval_grad!($ls)

# # eval_hess!(ls, z0)
# # @benchmark eval_hess!($ls)

# update!(ls)
# reshape(ls.θ, nx, nz)
# # @benchmark update!($ls, z) setup=(z=copy(z0))


text:
# see issue #151
@testset "HTML escape" begin
    st = read(joinpath(D, "151.md"), String)
    @test isapproxstr(st |> conv,
            """

            add OhMyREPL#master

AAA

            \"\"\"
                bar(x[, y])

            BBB

            # Examples
            ```jldoctest
            D
            ```
            \"\"\"
            function bar(x, y)
                ...
            end

For complex functions with multiple arguments use a argument list, also if there are many keyword arguments use <keyword arguments>:

            \"\"\"
                matdiag(diag, nr, nc; &ltkeyword arguments&gt)

            Create Matrix with number `vdiag` on the super- or subdiagonals and `vndiag`
            in the rest.

            # Arguments
            - `diag::Number`: `Number` to write into created super- or subdiagonal

            # Examples
            ```jldoctest
            julia> matdiag(true, 5, 5, sr=2, ec=3)
            ```
            \"\"\"
            function matdiag(diag::Number, nr::Integer, nc::Integer;)
                ...
            end

""")
end

# see issue #182
@testset "Code blocks" begin
    st = read(joinpath(D, "182.md"), String)
    @test isapproxstr(st |> conv, """

Code block:

            The average temperature is 19.5°C.

The end.

            """)
end

@testset "Table" begin
    fs1()
    st = """
        @def fd_rpath = "pages/pg1.html"
        A

        ### Title

        No. | Graph | Vertices | Edges
        :---: | :---------: | :------------: | :-----------------:
        1 | Twitter Social Circles | 81,306 | 1,342,310
        2 | Astro-Physics Collaboration | 17,903 | 197,031
        3 | Facebook Social Circles | 4,039 | 88,234

        C
        """

    if VERSION >= v"1.4.0-"
        @test isapproxstr(st |> seval, raw"""
A

Title

No.	Graph	Vertices	Edges

1	Twitter Social Circles	81,306	1,342,310

2	Astro-Physics Collaboration	17,903	197,031

3	Facebook Social Circles	4,039	88,234

C
""")
    else
        @test isapproxstr(st |> seval, raw"""
A

Title

                No.	Graph	Vertices	Edges

                1	Twitter Social Circles	81,306	1,342,310

                2	Astro-Physics Collaboration	17,903	197,031

                3	Facebook Social Circles	4,039	88,234

C

            """)
    end
end

@testset "Table-fs2" begin
    fs2()
    st = """
        @def fd_rpath = "pages/pg1.html"
        A

        ### Title

        No. | Graph | Vertices | Edges
        :---: | :---------: | :------------: | :-----------------:
        1 | Twitter Social Circles | 81,306 | 1,342,310
        2 | Astro-Physics Collaboration | 17,903 | 197,031
        3 | Facebook Social Circles | 4,039 | 88,234

        C
        """

    if VERSION >= v"1.4.0-"
        @test isapproxstr(st |> seval, raw"""
A

Title

No.	Graph	Vertices	Edges

1	Twitter Social Circles	81,306	1,342,310

2	Astro-Physics Collaboration	17,903	197,031

3	Facebook Social Circles	4,039	88,234

C
""")
    else
        @test isapproxstr(st |> seval, raw"""
A

Title

                No.	Graph	Vertices	Edges

                1	Twitter Social Circles	81,306	1,342,310

                2	Astro-Physics Collaboration	17,903	197,031

                3	Facebook Social Circles	4,039	88,234

C

            """)
    end
end

@testset "Auto title" begin
    # if no title is set, then the first header is used
    s = raw"""
    # AAA
    etc
    ~~~{{fill title}}~~~
    """ |> fd2html_td
    @test isapproxstr(s, raw"""
AAA
  etc AAA""")
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((1, 3, 2, false, true, false, false, false),),
    ((1, 3, 2, false, false, true, false, false),),
    ((1, 3, 2, true, false, false, true, false),),
    ((1, 3, 2, true, false, true, false, false),),
    ((1, 3, 2, true, false, false, true, true),),
    ]
# insts["fast"] = [] # TODO
insts["various"] = [
    ((5, 20, 20, false, true, false, false, false),),
    ((5, 20, 20, false, false, true, false, false),),
    ((5, 20, 20, true, false, false, true, false),),
    ((5, 20, 20, true, false, true, false, false),),
    ((5, 20, 20, true, false, false, true, true),),
    ((10, 30, 20, false, true, false, false, false),),
    ((10, 30, 20, false, false, true, false, false),),
    ((10, 30, 20, true, false, false, true, false),),
    ((10, 30, 20, true, false, true, false, false),),
    ((10, 30, 20, true, false, false, true, true),),
    ((8, 50, 30, false, true, false, false, false),),
    ((8, 50, 30, false, false, true, false, false),),
    ((8, 50, 30, true, false, false, true, false),),
    ((8, 50, 30, true, false, true, false, false),),
    ((8, 50, 30, true, false, false, true, true),),
    ]
return (SparseLMIJuMP, insts)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Allocations" begin
    # ensure that there are 0 allocations for all the benchmarks
    for b in Benchmarks.all_benchmarks()
        @test Benchmarks.runbenchmark(b, samples = 1, evals = 1).allocs == 0
    end
end # testset Allocations


text:
export write_vcf_head, write_vcf_line, write_txt_head, write_txt_line

function write_vcf_head(f::IO, rules)
    vars = Dict(v.name => v for r in rules for v in r.var)
    info = join(info_head(vars[anno]) for r in rules for anno in r.anno)
    filters = join(filter_head(r) for r in rules if r.soft && !isempty(r.filt))

    f << """
    ##fileformat=VCFv4.2
    ##fileDate=$(Dates.format(now(), "yyyymmdd"))
    ##source=FalconV0.0.1
    """ << info << filters << """
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
    ##FORMAT=<ID=RD,Number=1,Type=Integer,Description="Depth of reference-supporting bases">
    ##FORMAT=<ID=AD,Number=1,Type=Integer,Description="Depth of variant-supporting bases">
    #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tCFDNA
    """
end

function write_vcf_line(f::IO, bam::Bam, chr::Int32, m::Mut, info::String, filters::Vector{String}, data::NTuple{3, Int})
    f << (chr == -1 ? '*' : car(bam.refs[chr+1])) << '\t'
    f << (isa(m, Deletion) ? m.pos-1 : m.pos) << '\t'
    f << '.'   << '\t' # ID, don't know what's that though
    λ"""switch(m
        (isa . SNP)       (write f m.ref   '\t' m.alt   '\t')
        (isa . Insertion) (write f '.'     '\t' m.bases '\t')
        (isa . Deletion)  (write f m.bases '\t' '.'     '\t'))
    """
    f << '.' << '\t' # Qual
    f << (isempty(filters) ? "PASS" : join(filters, ';')) << '\t'
    f << info << '\t'
    f << "DP:RD:AD" << '\t'
    f << data[1] << ':' << data[2] << ':' << data[3] << '\n'
end

function write_txt_head(f::IO, info::String, filters::String)

end

function write_txt_line(f::IO, bam::Bam, chr::Int32, m::Mut, info::String, filters::Vector{String}, data::NTuple{3, Int})

end

anno_info{T}(f, var::Variable{Vector{T}})         = :( $f << $(Meta.quot(var.name)) << '='; join($f, $(var.name), ','); $f << ';' )
anno_info{T<:Integer}(f, var::Variable{T})        = :( $f << $(Meta.quot(var.name)) << '=' << $(var.name) << ';' )
anno_info{T<:AbstractFloat}(f, var::Variable{T})  = :( $f << $(Meta.quot(var.name)) << '=' << $(var.name) << ';' )
anno_info{T<:AbstractString}(f, var::Variable{T}) = :( $f << $(Meta.quot(var.name)) << '=' << $(var.name) << ';' )
anno_info(f, var::Variable{Bool})                 = :( ($(var.name) && $f << $(Meta.quot(var.name))); $f << ';' )

_info_head(var, n, t) = "##INFO=<ID=$(var.name),Number=$n,Type=$t,Description=\"$(var.desc)\">\n"
info_head{T<:Integer}(var::Variable{Vector{T}})        = _info_head(var, 'A', "Integer")
info_head{T<:AbstractFloat}(var::Variable{Vector{T}})  = _info_head(var, 'A', "Float")
info_head{T<:AbstractString}(var::Variable{Vector{T}}) = _info_head(var, 'A', "String")
info_head{T<:Integer}(var::Variable{T})                = _info_head(var,   1, "Integer")
info_head{T<:AbstractFloat}(var::Variable{T})          = _info_head(var,   1, "Float")
info_head{T<:AbstractString}(var::Variable{T})         = _info_head(var,   1, "String")
info_head(var::Variable{Bool})                         = _info_head(var,   0, "Flag")

filter_head(rule) = "##Filter=<ID=$(rule.name),Description=\"$(rule.desc)\">\n"


text:
# Generate figures:
include(joinpath(@__DIR__, "src", "plot_profiles.jl"))

prof_lib = joinpath(@__DIR__, "src", "generated_profile_library.md")
open(prof_lib, "w") do io

    println(io, "```@meta")
    println(io, "# ****************************************************")
    println(io, "# ****************************************************")
    println(io, "#     THIS FILE IS AUTO-GENERATED, PLEASE DO NOT EDIT.")
    println(io, "# ****************************************************")
    println(io, "# ****************************************************")
    println(io, "```")

    println(io, "# Profile library")

    println(io, "")
    println(io, "```@meta")
    println(io, "CurrentModule = AtmosphericProfilesLibrary")
    println(io, "```")

    println(io, "")
    println(io, "```@example")
    println(io, "include(joinpath(@__DIR__, \"plot_profiles.jl\"))")
    println(io, "```")
    println(io, "")

    println(io, "## z-profiles")
    for profile in z_profiles
        fname = nameof(profile.func)
        println(io, "```@docs")
        println(io, "$fname")
        println(io, "```")
        println(io, "![](z_$fname.png)")
        println(io, "")
    end

    if !isempty(t_profiles)
        println(io, "## t-profiles")
        for profile in t_profiles
            fname = nameof(profile.func)
            println(io, "```@docs")
            println(io, "$fname")
            println(io, "```")
            println(io, "![](t_$fname.png)")
            println(io, "")
        end
    end

    println(io, "## tz-profiles")
    for profile in tz_profiles
        fname = nameof(profile.func)
        println(io, "```@docs")
        println(io, "$fname")
        println(io, "```")
        println(io, "![](tz_$fname.png)")
        println(io, "")
    end
end

@assert isfile(prof_lib)


text:
using Test
using ThreeBodyTB
using Suppressor

#include("../includes_laguerre.jl")
#include("../Ewald.jl")

TESTDIR=ThreeBodyTB.TESTDIR

function loaddata(dirs; scf=true)
    tbc_list  = []
    dft_list = []

    for t in dirs
        #                println(t*"/qe.save")
        tfull = "$TESTDIR/"*t
        dft = ThreeBodyTB.QE.loadXML(tfull*"/qe.save")
        tbc = []
        tbc_scf = []

        try
            if scf
                tbc = ThreeBodyTB.TB.read_tb_crys_kspace("projham_K.xml.gz", directory=tfull)
                tbc_scf = ThreeBodyTB.SCF.remove_scf_from_tbc(tbc)
            else
                tbc_scf = ThreeBodyTB.TB.read_tb_crys_kspace("projham_K.xml.gz", directory=tfull)
            end
        catch
            println("failed to load $t $tfull")
#            tbc,tbck,projwarn = ThreeBodyTB.AtomicProj.projwfc_workf(dft, directory=tfull, skip_og=true, skip_proj=true, freeze=true, localized_factor = 0.15,only_kspace=true )
#            if scf
#                tbc_scf = ThreeBodyTB.SCF.remove_scf_from_tbc(tbck)
#            else
#                tbc_scf = tbck
#            end

        end

        push!(dft_list, dft)
        push!(tbc_list, tbc_scf)
    end
    return     tbc_list, dft_list
end

function test_force()

    @testset "testing force dimer" begin

        if true
#        @suppress begin
            ft = open("$TESTDIR/data_forces/fil_MgS_dimer", "r");
            dirst = readlines(ft);
            close(ft);

            #        println(dirst)

#            for scf = [false true]
            f_cart = zeros(2,3)
            f_cartFD = 0.0

            for scf = [false, true]
                @suppress begin
                    tbc_list, dft_list = loaddata(dirst, scf=scf);

                    println("TBC list ", tbc_list)
                    println("DFT list ", dft_list)
                    database_rec = ThreeBodyTB.FitTB.do_fitting_recursive(tbc_list,dft_list = dft_list,  fit_threebody=false, fit_threebody_onsite=false);

                x = 4;
                smearing = 0.01;

                    #en, tbc_x, flag = scf_energy(tbc_list[x].crys, database = database_rec)
                    #en, f_cart,stress = ThreeBodyTB.Force_Stress.get_energy_force_stress(tbc_x, database_rec,   smearing = smearing);

                    en, f_cart,stress = ThreeBodyTB.Force_Stress.get_energy_force_stress(tbc_list[x].crys, database_rec,   smearing = smearing);

                enFD, f_cartFD = ThreeBodyTB.Force_Stress.finite_diff(tbc_list[x].crys, database_rec,1, 3,   smearing = smearing);
                end

#                println(scf, " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  ", f_cartFD, "  ", f_cart[1,3])

                @test abs(f_cartFD - f_cart[1,3]) < 1e-3
#                @test abs(f_cartFD - f_cart_fft[1,3]) < 1e-3
                #            println("SCF $scf TEST1 finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
                #            println("TEST1 dft ref ", dft_list[x].forces[1,3])

            end
            #        x = 3;
            #        smearing = 0.01;
            #        en, f_cart,stress = Force_Stress.get_energy_force_stress(tbc_list[x].crys, database_rec,   smearing = smearing);
            #        enFD, f_cartFD = Force_Stress.finite_diff(tbc_list[x].crys, database_rec,1, 3,   smearing = smearing);

            #        @test abs(f_cartFD - f_cart[1,3]) < 1e-3
            #        println("TEST3 finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST3 dft ref ", dft_list[x].forces[1,3])

        end
    end
end

function test_stress()

    @testset "testing force znse" begin
        @suppress begin

            ft = open("$TESTDIR/data_forces/fil_MgS_znse", "r");
            dirst = readlines(ft);
            close(ft);

            #        println(dirst)

            tbc_list, dft_list = loaddata(dirst, scf=false);

            println("DFT_list", dft_list)
            println("TBC_list", tbc_list)

            database = ThreeBodyTB.FitTB.do_fitting_recursive(tbc_list, dft_list = dft_list, fit_threebody=false, fit_threebody_onsite=false, do_plot=false, niters=5);

            #        database = FitTB.do_fitting_recursive(tbc_list,dft_list = dft_list,  fit_threebody=true, fit_threebody_onsite=false);

            x = 1;
            smearing = 0.01;
            en, f_cart, stress = ThreeBodyTB.Force_Stress.get_energy_force_stress(tbc_list[x].crys, database,   smearing = smearing);

            enFD, f_cartFD = ThreeBodyTB.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 3,   smearing = smearing);

            #        println("TEST force finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].forces[1,3])
            @test abs(f_cartFD - f_cart[1,3]) < 1e-4

            x=1
            enFD, stressFD = ThreeBodyTB.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 1, stress_mode=true,  smearing = smearing);

            #        println("TEST stress11 finite diff: ", stressFD , " autodiff:   ", stress[1,1], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].stress[1,1])
            @test abs(stressFD - stress[1,1]) < 1e-5

            x=1
            enFD, stressFD = ThreeBodyTB.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 2, stress_mode=true,  smearing = smearing);

            #        println("TEST stress12 finite diff: ", stressFD , " autodiff:   ", stress[1,2], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].stress[1,2])
            @test abs(stressFD - stress[1,2]) < 1e-5

            #        println("done")
        end
    end
end

test_force()

#println("sleep ")
#sleep(3)

test_stress()


text:
##

@info("-------- TEST FILTERING MECHANISM ---------")
using SHIPs, JuLIP, BenchmarkTools, LinearAlgebra, Test, Random, StaticArrays
using SHIPs: eval_basis!, eval_basis, PolyCutoff1s, PolyCutoff2s
using JuLIP.MLIPs: IPSuperBasis
using JuLIP.Testing: print_tf
using Printf

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N)
randR(N, syms) = randR(N)[1], rand( Int16.(atomic_number.(syms)), N )
function randiso()
   K = @SMatrix rand(3,3)
   K = K - K'
   Q = rand([-1,1]) * exp(K)
end
function randiso(Rs)
   Q = randiso()
   return [ Q * R for R in shuffle(Rs) ]
end

##

@info("Construct a 5B 🚢")
trans = PolyTransform(2, 1.0)
cutf = PolyCutoff2s(2, 0.5, 3.0)
ship = SHIPBasis(SparseSHIP(5,  12; wL = 1.5), trans, cutf,
                                    filter=false)

##
@info("Select a basis group and show it doesn't have full rank:")
maxgrp = maximum(SHIPs.alllen_bgrp(ship, 1))
@show maxgrp
igrp = findfirst(SHIPs.alllen_bgrp(ship, 1) .== maxgrp)
Igrp = SHIPs.I_bgrp(ship, igrp, 1)

G = zeros(length(Igrp), length(Igrp))
nsamples = 100 * length(Igrp)
Zs = zeros(SHIPs.IntS, 5)
for n = 1:nsamples
   Rs, Zs = randR(5)
   B = SHIPs.eval_basis(ship, Rs, Zs, 0)
   b = B[Igrp]
   global G += b * b' / nsamples
end

println(@test rank(G) < maxgrp)
@show rk = rank(G)

@show svd(G).S

##
@info("Now filter that basis and show that this basis group has reduced to length 1")

fship = SHIPs.filter_rpi_basis(ship, 1_000)
println(@test SHIPs.len_bgrp(fship, igrp, 1) == rk)
@show length(ship), length(fship)

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
export DeflatedNEP
export DeflatedNEPMM
export DeflatedGenericNEP
export DeflatedSPMF
export deflate_eigpair
export get_deflated_eigpairs

"""
    struct DeflatedNEPMM <: NEP

Represents a deflated NEP where the compute functions
are carried out via `compute_MM`. See further
documentation [`deflate_eigpair`](@ref).

"""
struct DeflatedNEPMM <: NEP
    orgnep::NEP
    S0
    V0
end

"""
    struct DeflatedSPMF <: AbstractSPMF

Represents a deflated NEP based on transforming
the deflated NEP to SPMF-form. See further
documentation [`deflate_eigpair`](@ref).
"""
struct DeflatedSPMF{T,NEP1,NEP2} <: AbstractSPMF{T}
    orgnep::AbstractSPMF{T}
    spmf::SPMFSumNEP{NEP1,NEP2}
    S0
    V0
end

"""
    struct DeflatedGenericNEP <: NEP

Represents a deflated NEP based on transforming
carrying out compute functions of derivatives by
binomial expansions. See further
documentation [`deflate_eigpair`](@ref).
"""
struct DeflatedGenericNEP <: NEP
    orgnep::NEP
    S0
    V0
end

DeflatedNEP=Union{DeflatedNEPMM,DeflatedGenericNEP,DeflatedSPMF{T} where T<:AbstractMatrix};

function size(nep::DeflatedNEP,dim=-1)
    n=size(nep.orgnep,1)+size(nep.V0,2);
    if (dim==-1)
        return (n,n)
    else
        return n
    end
end

## DeflatedGenericNEP optimized routines based on binomial expansion for derivative
function compute_Mlincomb(nep::DeflatedGenericNEP,λ::Number,
                 V::AbstractVecOrMat,a::Vector=ones(eltype(V),size(V,2)))

    X=nep.V0;
    S=nep.S0;
    T=promote_type(typeof(λ),eltype(X),eltype(S),eltype(a),eltype(V));
    # It probably only makes sense to do a factorization
    # if c*(p*p*p+k*k*p*p) < k*k*p*p*p
    F=factorize(λ*I-S);
    k=size(V,2);
    Xhat=X/F;
    n0=size(nep.orgnep,1);

    p=size(S,1);
    # precompute terms with inverses:
    Q=Vector{Matrix{T}}(undef,k)
    for i=0:k-1;
        QQ=zeros(T,p,k)
        QQ[:,i+1]=V[(n0+1):end,i+1];
        for j=(i-1):-1:0;
          QQ[:,j+1]=F\QQ[:,j+2];
          #QQ[:,j+1]=(λ*I-S)^(-(i-j))*V[(n0+1):end,i+1];
        end
        Q[i+1]= QQ
    end

    Z=zeros(T,n0,k);
    for j=0:k-1
        z=zeros(T,n0);
        for i=j:k-1
            #Vnew= (λ*I-S)^(-(i-j))*V[(n0+1):end,i+1]w
            Vnew = Q[i+1][:,j+1];
            factor=((-1)^(i-j))*(a[i+1]*factorial(i)/factorial(j));
            z+=factor*Xhat*Vnew;
        end
        Z[:,j+1]=z;
    end
    Vnew=V[1:n0,:]*Diagonal(a)+Z;
    z_top=compute_Mlincomb(nep.orgnep,λ,Vnew);

    z_bottom=X'*V[1:n0,1]*a[1];
    return [z_top ; z_bottom];
end

compute_Mder(nep::DeflatedGenericNEP,λ::Number)=compute_Mder(nep,λ,0)
function compute_Mder(nep::DeflatedGenericNEP,λ::Number,der::Integer)
    X=nep.V0;
    S=nep.S0;
    T=promote_type(typeof(λ),eltype(X),eltype(S))

    n0=size(nep.orgnep,1);
    p=size(S,1);

    Q = deflated_nep_compute_Q(nep,λ,der)

    M0=compute_Mder(nep.orgnep,λ,der);
    if (M0 isa SparseMatrixCSC)
        # return a sparse matrix
        (II,JJ,VV)=findnz(M0);
        (II2,JJ2,VV2)=findnz(sparse(Q))
        JJ2=JJ2 .+  n0;
        II_merged=[II;II2];
        JJ_merged=[JJ;JJ2];
        VV_merged=[VV;VV2];

        # If zero'th derivative add (2,1)-block
        if (der==0)
            (II3,JJ3,VV3)=findnz(sparse(Matrix(X')));
            II_merged=[II_merged;II3 .+ n0];
            JJ_merged=[JJ_merged;JJ3];
            VV_merged=[VV_merged;VV3];
        end
        return sparse(II_merged,JJ_merged,VV_merged,n0+p,n0+p);
    else
        # Return a dense matrix
        if (der==0)
            return [M0 Q ; X'  zeros(T,p,p)];
        else
            return [M0 Q ; zeros(T,p,n0)  zeros(T,p,p)];
        end
    end
end

function deflated_nep_compute_Q(nep::DeflatedNEP,λ::Number,der::Integer)
    X=nep.V0
    S=nep.S0
    T=promote_type(typeof(λ),eltype(X),eltype(S))

    n0=size(nep.orgnep,1)
    p=size(S,1)

    Q=zeros(T,n0,p)
    F=factorize(λ*I-S)
    Vnew = X
    for i=der:-1:0
        # Equivalent to: Vnew= X*(λ*I-S)^(-(der-i+1))
        Vnew = Vnew/F
        factor=((-1)^(der-i))*(factorial(der)/factorial(i))
        # Equivalent to: Mi=compute_Mder(nep.orgnep,λ,i); Q+=Mi*(Vnew*factor);
        for j = 1:p
            Q[:,j] = Q[:,j] + compute_Mlincomb(nep.orgnep, λ, Vnew[:,j], [factor], i)
        end
    end
    return Q
end

## The DeflatedSPMF just delegates to the nep.spmf
compute_Mlincomb(nep::DeflatedSPMF,λ::Number,V::AbstractVecOrMat,a::Vector)= compute_Mlincomb(nep.spmf,λ,V,a);
compute_Mlincomb(nep::DeflatedSPMF,λ::Number,V::AbstractVecOrMat)= compute_Mlincomb(nep.spmf,λ,V);
compute_Mder(nep::DeflatedSPMF,λ::Number)=compute_Mder(nep.spmf,λ,0)
compute_Mder(nep::DeflatedSPMF,λ::Number,der::Integer)=compute_Mder(nep.spmf,λ,der)
compute_MM(nep::DeflatedSPMF,par...)=compute_MM(nep.spmf,par...)
get_Av(nep::DeflatedSPMF)=get_Av(nep.spmf)
get_fv(nep::DeflatedSPMF)=get_fv(nep.spmf)

## The DeflatedNEPMM & DeflatedGenericNEP calls corresponding MM-functions
function compute_MM(nep::Union{DeflatedNEPMM,DeflatedGenericNEP},S,V)
    orgnep=nep.orgnep;
    n0=size(orgnep,1);
    S0=nep.S0; V0=nep.V0
    p0=size(S0,1); p=size(S,1);
    V1=view(V,1:n0,            1:size(V,2))
    V2=view(V,(n0+1):size(V,1),1:size(V,2))
    Stilde=hcat(vcat(S0,zeros(p,p0)),vcat(V2,S))
    Vtilde=hcat(V0,V1);
    R=compute_MM(orgnep, Stilde,Vtilde);
    return vcat(R[1:n0,(size(nep.S0,1)+1):end],V0'*V1);
end
# Use the MM to compute Mlincomb for DeflatedNEPMM
compute_Mlincomb(nep::DeflatedNEPMM,λ::Number,
                 V::AbstractVecOrMat,a::Vector=ones(eltype(V),size(V,2)))=
             compute_Mlincomb_from_MM(nep,λ,V,a)

function compute_Mder(nep::DeflatedNEPMM,λ::Number,i::Integer=0)
    return compute_Mder_from_MM(nep,λ,i);
end

# Creates and SPMF with deflated `S0` and `V0`.
function create_spmf_dnep(nep::AbstractSPMF,S0,V0)
    Av_org=get_Av(nep);
    fv_org=get_fv(nep);
    m=size(fv_org,1);
    p=size(V0,2);
    n0=size(nep,1);

    m1=m;     # size of "old" part
    m2=m*p+1; # size of "deflation" part

    # spmf1: Create the "old" part
    A1=Vector{eltype(Av_org)}(undef,m1);
    for k=1:m
        A0k=Av_org[k];
        if (eltype(A1) <:  SparseMatrixCSC)
            (II,JJ,VV)=findnz(A0k)
            A1[k]=sparse(II,JJ,VV,n0+p,n0+p);
        else
            A1[k]=zeros(eltype(A0k),n0+p,n0+p)
            A1[k][1:n0,1:n0]=A0k;
        end
    end
    spmf1=SPMF_NEP(A1,fv_org,check_consistency=false)
    # spmf2: Create the additional deflation terms:

    #    Promote rules for the eltype:
    #    We may need to increase the eltype type size, since S0,V0 can be complex
    T=promote_type(eltype(V0),eltype(S0),eltype(Av_org[1]));
    local T_LowRankFactor;
    if (eltype(Av_org) <: SparseMatrixCSC)
        T_LowRankFactor=SparseMatrixCSC{T,Int};
    else
        T_LowRankFactor=Matrix{T};
    end
    L2=Vector{T_LowRankFactor}(undef,m2);
    U2=Vector{T_LowRankFactor}(undef,m2);
    fv2=Vector{Function}(undef,m2);
    (λtmp,X)=eigen(S0);
    λ::Vector{T}=λtmp[:]; # Ensure type
    count=0;
    for i=1:p
        y=(V0*(X[:,i]));
        ei=zeros(p); ei[i]=1;
        x=(ei'/X);
        for r=1:m
            count=count+1;
            # This will automatically convert to sparse / full
            L2[count] = reshape([(Av_org[r]*y) ;zeros(p)],n0+p,1);
            U2[count] = reshape([zeros(n0);x'],n0+p,1);
            fv2[count]=S-> (S-λ[i]*one(S))\fv_org[r](S);
        end
    end
    # The constant term
    L2[m*p+1]=[zeros(n0,p);Matrix{T}(I,p,p)]
    U2[m*p+1]=[Matrix(V0);zeros(p,p)]
    fv2[m*p+1]= S->one(S);
    spmf2=LowRankFactorizedNEP(L2,U2,fv2);

    return SumNEP(spmf1,spmf2);
end

"""
    normalize_schur_pair!(S,V)

Makes the Schur pair ``(S,V)`` in the sense that `V'*V` is the
identity matrix. This will not work if `size(V,2)>size(V,1)`.

"""
function normalize_schur_pair!(S,V)
    if (size(V,2)>size(V,1))
        @warn "Cannot normalize short and skinny V-matrices."
    else
        (QQ,RR)=qr(V);
        V[:]=Matrix(QQ); # Use skinny QR-factorization
        S[:]=(RR*S)/RR;
    end
end

# Determine and verify that the deflate mode is correct.
function verify_deflate_mode(nep::NEP,mode)
    if (mode==:Auto)
        if (nep isa AbstractSPMF)
            mode=:SPMF
        else
            mode=:Generic
        end
    end
    if ((mode == :SPMF) || (mode == :SPMFPlain)) &&
        !(nep isa AbstractSPMF)
        error("SPMF-mode only possible for `AbstractSPMF`-NEPs")
    end
    return mode;
end
function verify_deflate_mode(nep::DeflatedNEP,mode)
    if (nep isa DeflatedSPMF && ((mode == :SPMF) || (mode == :Auto)) )
       return :SPMF
    elseif (nep isa DeflatedNEPMM && ((mode == :MM) || (mode == :Auto)) )
       return :MM
    elseif (nep isa DeflatedGenericNEP && ((mode == :Generic) || (mode == :Auto)) )
       return :Generic
    else
       error("Unknown mode / type");
    end
end

"""
    dnep=deflate_eigpair(orgnep::NEP,λ,v,[mode=:Auto])

This function creates a deflated NEP based on ``(λ,v)``, which
are assumed to an eigenpair of `nep`. Effectively,
the function will return `dnep::NEP` which has the same
solutions as orgnep, except those corresponding to ``(λ,v)``.
Deflation is typically used to avoid reconvergence.

If `orgnep` is a `DeflatedNEP`, the `orgnep` the
deflation in `orgnep` will be updated.

The `mode` kwarg can be `:Auto`, `:Generic`, `:SPMF`,
`:MM`. This specifies how the deflated NEP should be represented.
Which mode is the most efficient depends on
many problem properties.
If the original NEP is an `AbstractSPMF` with only
a few terms, `mode=:SPMF` may be efficient. The SPMF-mode
is based on a diagonalization of the deflated invariant
pair and is not necessarily robust when you deflate eigenvalues
near to each other. When `mode=:MM`
is used, all compute functions are implemented via calls
to the `compute_MM`. This can work well for small dense
problems. The `:Generic` is based on an explicit derivation
of the problem (via binomial expansion) which can be
efficient if low order derivates are needed. If
`:Auto` is selected, NEP-PACK tries to determine which
one is the most efficient based on the `orgnep`.

# Example:
```julia-repl
julia> nep=nep_gallery("dep0");
julia> (λ,v)=newton(nep,v=ones(size(nep,1)));
julia> dnep=deflate_eigpair(nep,λ,v)
julia> (λ2,v2)=augnewton(dnep,v=ones(size(dnep,1)));  # this converges to different eigval
julia> using LinearAlgebra;
julia> minimum(svdvals(compute_Mder(nep,λ2)))
2.5161012836775824e-17
```
The function [`get_deflated_eigpairs()`](@ref) extracts the eigenpairs that
have been deflated. The returned pairs are eigenpairs of
the original NEP:
```julia-repl
julia> dnep=deflate_eigpair(dnep,λ2,v2);
julia> (D,V)=get_deflated_eigpairs(dnep)
julia> norm(compute_Mlincomb(nep,D[1],V[:,1]))
2.3970746442479104e-16
julia> norm(compute_Mlincomb(nep,D[2],V[:,2]))
8.101585801848734e-16
```

# References
* C. Effenberger, Robust solution methods for nonlinear eigenvalue problems, PhD thesis, 2013, EPF Lausanne
"""
function deflate_eigpair(nep::NEP,λ,v;mode=:Auto)
    mode=verify_deflate_mode(nep,mode);
    n=size(nep,1);
    S0=reshape([λ],1,1);
    V0=reshape(v,n,1);
    normalize_schur_pair!(S0,V0);
    # Create it based on nep
    return internal_create_deflated_nep(nep,S0,V0,mode);
end

function deflate_eigpair(nep::DeflatedNEP,λ,v;mode=:Auto)
    mode=verify_deflate_mode(nep,mode)

    T=promote_type(typeof(λ),eltype(v),eltype(nep.V0),eltype(nep.S0));

    n=size(nep.orgnep,1);
    p0=size(nep.V0,2);
    # fetch pair + expand with λ v
    V1=zeros(T,n,p0+1);
    S1=zeros(T,p0+1,p0+1);
    V1[1:n,1:end-1]=nep.V0
    V1[1:n,end]=v[1:n];
    S1[1:end-1,1:end-1]=nep.S0;
    S1[1:end,end]=[v[n+1:end];λ];

    # normalize schur pair
    normalize_schur_pair!(S1,V1);
    # Create it based on nep.orgnep
    return internal_create_deflated_nep(nep.orgnep,S1,V1,mode);
end

# Internal helper function to create a nep based on mode
function internal_create_deflated_nep(nep,S1,V1,mode)
    # create new DeflatedNEP
    if (mode==:MM)
        newnep=DeflatedNEPMM(nep,S1,V1);
        return newnep;
    elseif (mode==:SPMF)
        spmf=create_spmf_dnep(nep,S1,V1);;
        newnep=DeflatedSPMF(nep,spmf,S1,V1);
        return newnep;
    elseif (mode==:Generic)
        newnep=DeflatedGenericNEP(nep,S1,V1);
        return newnep;
    end
end

"""
    (D,V)=get_deflated_eigpairs(S,V,n)
    (D,V)=get_deflated_eigpairs(dnep::DeflatedNEP [λ,v])

Returns a vector of eigenvalues `D` and a matrix with corresponding
eigenvectors `V` of the invariant pair `S,V`.
The eigenpairs correspond to the original problem, underlying the `DeflatedNEP`.
The optional parameters `λ,v` allows the inclusion
of an additional eigpair. Essentially, the optional parameters
are the expanding the deflation and the running `get_deflated_eigpairs`
 with kwarg, i.e.,
```julia
(D,V)=get_deflated_eigpairs(deflate_eigpair(dnep,λ,v))`
```
See example in [`deflate_eigpair`](@ref).
"""
function get_deflated_eigpairs(S,V,n)
    (D,X)=eigen(S);
    return D,V[1:n,:]*X;
end
function get_deflated_eigpairs(nep::DeflatedNEP)
    return get_deflated_eigpairs(nep.S0,nep.V0,size(nep.orgnep,1))
end
function get_deflated_eigpairs(nep::DeflatedNEP,λ,v)
   T=promote_type(typeof(λ),eltype(v),eltype(nep.V0),eltype(nep.S0));
   p0=size(nep.V0,2);
   n=size(nep.orgnep,1);
   # fetch pair + expand with λ v
   V1=zeros(T,n,p0+1);
   S1=zeros(T,p0+1,p0+1);
   V1[1:n,1:end-1]=nep.V0
   V1[1:n,end]=v[1:n];
   S1[1:end-1,1:end-1]=nep.S0;
   S1[1:end,end]=[v[n+1:end];λ];
   V=V1
   S=S1;
   (D,X)=eigen(S);
   return D,V[1:size(nep.orgnep,1),:]*X;
end


text:
using Catlab
using Catlab.Theories
using Catlab.Syntax

using Metatheory, Metatheory.EGraphs
using TermInterface

@metatheory_init ()

abstract type CatType end
struct ObType <: CatType
    ob
    mod
end
struct HomType <: CatType
    dom
    codom
    mod
end

# Custom type APIs for the GATExpr
using Metatheory.TermInterface
TermInterface.gethead(t::ObExpr) = :call
TermInterface.getargs(t::ObExpr) = [head(t), t.args...]
TermInterface.gethead(t::HomExpr) = :call
TermInterface.getargs(t::HomExpr) = [head(t), t.args...]

# Type information will be stored in the metadata
function TermInterface.metadata(t::HomExpr)
    return HomType(t.type_args[1], t.type_args[2], typeof(t).name.module)
end
TermInterface.metadata(t::ObExpr) = ObType(t, typeof(t).name.module)
TermInterface.isterm(t::GATExpr) = true
TermInterface.arity(t::GATExpr) = length(getargs(t))

struct CatlabAnalysis <: AbstractAnalysis end
function EGraphs.make(an::Type{CatlabAnalysis}, g::EGraph, n::ENode{T}) where T
    !(T <: GATExpr) && return T
    return metadata(n)
end
EGraphs.join(an::Type{CatlabAnalysis}, from, to) = from
EGraphs.islazy(x::Type{CatlabAnalysis}) = false

function infer(t::GATExpr)
    g = EGraph(t)
    analyze!(g, CatlabAnalysis)
    getdata(geteclass(g, g.root), CatlabAnalysis)
end

function EGraphs.extractnode(g::EGraph, n::ENode{T}, extractor::Function) where {T <: ObExpr}
    @assert n.head == :call
    return metadata(n).ob
end

function EGraphs.extractnode(g::EGraph, n::ENode{T}, extractor::Function) where {T <: HomExpr}
    @assert n.head == :call
    nargs = extractor.(n.args)
    nmeta = metadata(n)
    return nmeta.mod.Hom{nargs[1]}(nargs[2:end], GATExpr[nmeta.dom, nmeta.codom])
end

# ==============================================================================

using MatchCore

datasym(x::Symbol) = Symbol(String(x) * "_data")
extrsym(x::Symbol) = Symbol(String(x) * "_extr")

function build_rhs(x::Expr, pvars, mod)
    if Meta.isexpr(x, :call)
        if x.args[1] == :munit && length(x.args) == 1
            mod.munit(mod.Ob)
        else
            Expr(x.head, getfield(mod, x.args[1]), map(y -> build_rhs(y, pvars, mod), x.args[2:end])...)
        end
    else
        Expr(x.head, map(y -> build_rhs(y, pvars, mod), x.args)...)
    end
end
function build_rhs(x, pvars, mod)
    if x ∈ pvars
        extrsym(x)
    else
        x
    end
end
function gen_rule(axiom::Catlab.GAT.AxiomConstructor, mod; righttoleft=false)
    # left to right
    @assert axiom.name == :(==)

    ax_left = axiom.left
    ax_right = axiom.right
    if righttoleft
        ax_left = axiom.right
        ax_right = axiom.left
    end

    lhs = Pattern(ax_left, mod)
    l_pvars = patvars(lhs)

    rhs = build_rhs(ax_right, l_pvars, mod)
    # println(rhs)

    lines = []

    eq_ctx = Dict{Symbol, Vector{Any}}()
    for patvar in l_pvars
        # retrieve the catlab data
        data_var = datasym(patvar)
        data_expr = :($data_var = getdata($patvar, CatlabAnalysis))
        push!(lines, data_expr)
        # push!(lines, :(println($data_var)))

        extr_var = extrsym(patvar)
        extr_expr = :($extr_var = extract!(_egraph, astsize; root=($patvar).id))
        push!(lines, extr_expr)
        # push!(lines, :(println($extr_var)))

        ctxval = axiom.context[patvar]
        # TODO use GATTheory.types
        @smatch ctxval begin
            :Ob => begin
                check_type_line = :(!($data_var isa ObType) && (return _lhs_expr))
                aset = get!(()->[], eq_ctx, patvar)
                push!(lines, check_type_line)
                push!(aset, :($(data_var).ob))
            end
            :(Hom($(A::Symbol), $(B::Symbol))) => begin
                aset = get!(()->[], eq_ctx, A)
                bset = get!(()->[], eq_ctx, B)
                push!(aset, :($(data_var).dom))
                push!(bset, :($(data_var).codom))
                check_type_line = :(!($data_var isa HomType) && (return _lhs_expr))
                push!(lines, check_type_line)
            end
            _ => error("unrecognized GAT type context $patvar => $ctxval")
        end
    end

    for (ctxvar, eqset) in eq_ctx
        if ctxvar ∉ l_pvars
            push!(lines, :($ctxvar = $(eqset[1])))
        end
    end

    # conjunction of equalities needed
    conjunction = []

    for (ctxvar, eqset) in eq_ctx
        unique!(eqset)
        c = []
        if length(eqset) < 2
            continue
        end
        fst = first(eqset)
        for other in eqset[2:end]
            push!(c, :($fst == $other))
        end
        append!(conjunction, c)
    end

    if !isempty(conjunction)
        conj_expr = foldl((x,y) -> :($x && $y), conjunction)

        the_big_if = :(if $conj_expr
            # WORKAROUND FOR RuntimeGeneratedFunctions.jl `id` bug
            # return $(evalmod).eval($(Meta.quot(ax_right)))
            return $rhs
        else
            return _lhs_expr end) |> Metatheory.Util.rmlines
        push!(lines, the_big_if)
    else
        push!(lines, :(return $rhs))
    end

    block = Expr(:block, lines...)

    DynamicRule(lhs, block)
end

# test
tt = theory(SymmetricMonoidalCategory)
ax = tt.axioms[2]

gen_rule(tt.axioms[2], @__MODULE__)

# Generate a theory from a syntax system
function gen_theory(m::Module)
    gat_theory = theory(m.theory())
    mt_theory = Rule[]
    for axiom in gat_theory.axioms
        push!(mt_theory, gen_rule(axiom, m))
        push!(mt_theory, gen_rule(axiom, m, righttoleft=true))
    end
    mt_theory
end

# ====================================================
# TEST

# WE HAVE TO REDEFINE THE SYNTAX TO AVOID ASSOCIATIVITY AND N-ARY FUNCTIONS
import Catlab.Theories: id, compose, otimes, ⋅, braid, σ, ⊗, Ob, Hom
@syntax SMC{ObExpr,HomExpr} SymmetricMonoidalCategory begin
end

function simplify(ex, syntax)
    t = gen_theory(syntax)
    g = EGraph(ex)
    analyze!(g, CatlabAnalysis)
    params=SaturationParams(timeout=3)
    saturate!(g, t, params; mod=@__MODULE__)
    extract!(g, astsize)
end

Metatheory.options.printiter = true
Metatheory.options.verbose = true

A, B, C = Ob(SMC, :A, :B, :C)
f = Hom(:f, A, B)

ex = f ⋅ id(B)
simplify(ex, SMC) == f

ex = id(A) ⋅ id(A) ⋅ f ⋅ id(B)
simplify(ex, SMC) == f

ex = σ(A,B) ⋅ σ(B,A)
simplify(ex, SMC) == id(A ⊗ B)

# ======================================================
# another test

using Catlab.Graphics

l = (σ(C,B) ⊗ id(A)) ⋅ (id(B) ⊗ σ(C,A)) ⋅ (σ(B,A) ⊗ id(C))
r = (id(C) ⊗ σ(B,A)) ⋅ (σ(C,A) ⊗ id(B)) ⋅ (id(A) ⊗ σ(C,B))

to_graphviz(l)
to_graphviz(r)

g = EGraph()
analyze!(g, CatlabAnalysis)
l_ec, _ = addexpr!(g, l)
r_ec, _ = addexpr!(g, r)

in_same_class(g, l_ec, r_ec)

saturate!(g, gen_theory(SMC), SaturationParams(timeout=1, eclasslimit=6000); mod=@__MODULE__)

ll = extract!(g, astsize; root=l_ec.id)
rr = extract!(g, astsize; root=r_ec.id)

# ======================================================
# another test

# WE HAVE TO REDEFINE THE SYNTAX TO AVOID ASSOCIATIVITY AND N-ARY FUNCTIONS
import Catlab.Theories: id, compose, otimes, ⋅, braid, σ, ⊗, Ob, Hom, pair, proj1, proj2
@syntax BPC{ObExpr,HomExpr} BiproductCategory begin
end
A, B, C = Ob(BPC, :A, :B, :C)
f = Hom(:f, A, B)
k = Hom(:k, B, C)

l = Δ(A) ⋅ (delete(A) ⊗ id(A))
r = id(A)

g = EGraph(l)
analyze!(g, CatlabAnalysis)

saturate!(g, gen_theory(BPC), SaturationParams(timeout=1, eclasslimit=6000); mod=@__MODULE__)

extract!(g, astsize)

# ======================================================
# another test

l = σ(A, B ⊗ C)
# r = σ(B,A) ⊗ id(C)
r = (σ(A,B) ⊗ id(C)) ⋅ (id(B) ⊗ σ(A,C))
# r = σ(B ⊗ C, A)

to_composejl(l)
to_composejl(r)

g = EGraph(ex)
analyze!(g, CatlabAnalysis)
l_ec, _ = addexpr!(g, l)
r_ec, _ = addexpr!(g, r)

saturate!(g, gen_theory(SMC), SaturationParams(timeout=1, eclasslimit=6000); mod=@__MODULE__)

extract!(g, astsize; root=l_ec.id)

extract!(g, astsize; root=r_ec.id)

# ====================================================
# TEST
cc = gen_theory(FreeCartesianCategory)

A, B, C = Ob(FreeCartesianCategory, :A, :B, :C)
f = Hom(:f, A, B)

g = EGraph()
analyze!(g, CatlabAnalysis)
ex = id(A) ⊗ id(B)
to_composejl(ex)

l_ec, _ = addexpr!(g, ex)
saturate!(g, cc, SaturationParams(timeout=2); mod=@__MODULE__)
extract!(g, astsize; root=l_ec.id)

ex = pair(proj1(A, B), proj2(A, B))
to_composejl(ex)
r_ec, _ = addexpr!(g, ex)
saturate!(g, cc; mod=@__MODULE__)
extract!(g, astsize; root=r_ec.id)


text:
@testset "NLP with parameters" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawOptimizerAttribute("print_level"), 0)
    cached =
        () -> MOI.Bridges.full_bridge_optimizer(
            MOIU.CachingOptimizer(
                MOIU.UniversalFallback(MOIU.Model{Float64}()),
                ipopt,
            ),
            Float64,
        )
    POI_cached_optimizer() = ParametricOptInterface.Optimizer(cached())

    model = Model(() -> POI_cached_optimizer())

    @variable(model, x)
    @variable(model, y)
    @variable(model, z in ParametricOptInterface.Parameter(10))
    @constraint(model, x + y >= z)
    @NLobjective(model, Min, x^2 + y^2)

    optimize!(model)
    objective_value(model)
    @test value(x) ≈ 5
    MOI.get(model, ParametricOptInterface.ParameterDual(), z)
    MOI.set(model, ParametricOptInterface.ParameterValue(), z, 2.0)
    optimize!(model)
    @test objective_value(model) ≈ 2 atol = 1e-3
    @test value(x) ≈ 1

    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawOptimizerAttribute("print_level"), 0)
    model = Model(() -> ParametricOptInterface.Optimizer(ipopt))

    @variable(model, x)
    @variable(model, z in ParametricOptInterface.Parameter(10))
    MOI.set(
        model,
        ParametricOptInterface.ConstraintsInterpretation(),
        ParametricOptInterface.ONLY_BOUNDS,
    )
    @constraint(model, x >= z)
    @NLobjective(model, Min, x^2)

    optimize!(model)
    objective_value(model)
    @test value(x) ≈ 10
    MOI.get(model, ParametricOptInterface.ParameterDual(), z)
    MOI.set(model, ParametricOptInterface.ParameterValue(), z, 2.0)
    optimize!(model)
    @test objective_value(model) ≈ 4 atol = 1e-3
    @test value(x) ≈ 2
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
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.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting marginal electricity price for each model zone and time step. Marginal electricity price is equal to the dual variable of the load balance constraint. If GenX is configured as a mixed integer linear program, then this output is only generated if `WriteShadowPrices` flag is activated. If configured as a linear program (i.e. linearized unit commitment or economic dispatch) then output automatically available.
"""
function write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	## Extract dual variables of constraints
	# Electricity price: Dual variable of hourly power balance constraint = hourly price
	dfPrice = DataFrame(Zone = 1:Z) # The unit is $/MWh

	# Dividing dual variable for each hour with corresponding hourly weight to retrieve marginal cost of generation
	if setup["ParameterScale"] == 1
		dfPrice = hcat(dfPrice, DataFrame(transpose(dual.(EP[:cPowerBalance])./inputs["omega"]*ModelScalingFactor), :auto))
	else
		dfPrice = hcat(dfPrice, DataFrame(transpose(dual.(EP[:cPowerBalance])./inputs["omega"]), :auto))
	end

	auxNew_Names=[Symbol("Zone");[Symbol("t$t") for t in 1:T]]
	rename!(dfPrice,auxNew_Names)

	## Linear configuration final output
	CSV.write(joinpath(path, "prices.csv"), dftranspose(dfPrice, false), writeheader=false)
	return dfPrice
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
include("imports.jl")

@testset ExtendedTestSet "Learner callback utilities" begin
    @testset ExtendedTestSet "setcallbacks!" begin
        learner = Learner(nothing, nothing, nothing, nothing; usedefaultcallbacks = false)
        @test length(learner.callbacks.cbs) == 0
        FluxTraining.setcallbacks!(learner, [Metrics()])
        @test length(learner.callbacks.cbs) == 1
    end

    @testset ExtendedTestSet "addcallback!" begin
        learner = Learner(nothing, nothing, nothing, nothing; usedefaultcallbacks = false)
        @test length(learner.callbacks.cbs) == 0
        FluxTraining.addcallback!(learner, Metrics())
        @test length(learner.callbacks.cbs) == 1
    end

    @testset ExtendedTestSet "getcallback" begin
        learner = Learner(nothing, nothing, nothing, nothing; usedefaultcallbacks = false)
        @test isnothing(FluxTraining.getcallback(learner, Metrics))
        FluxTraining.addcallback!(learner, Metrics())
        @test !isnothing(FluxTraining.getcallback(learner, Metrics))
    end

    @testset ExtendedTestSet "replacecallback!" begin
        learner = Learner(nothing, nothing, nothing, nothing, Metrics(accuracy); usedefaultcallbacks = false)
        @test length(learner.callbacks.cbs[1].metrics) == 2
        FluxTraining.replacecallback!(learner, Metrics())
        @test length(learner.callbacks.cbs[1].metrics) == 1
    end
end