text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

"""
    data = SimpleIdealGasMediumData(;cp_const=nothing, R_gas=nothing, MM_const=nothing,
                                     eta_const=nothing, lambda_const=nothing, T_min=nothing,
                                     T_max=nothing, T0=nothing)

Generate a `SimpleIdealGasMediumData` object containing the data
for a SimpleIdealGas medium.
"""
mutable struct SimpleIdealGasMediumData
    cp_const::Float64       # Constant specific heat capacity at constant pressure
    cv_const::Float64       # Constant specific heat capacity at constant volume
    R_gas::Float64          # Medium specific gas constant
    MM_const::Float64       # Molar mass
    eta_const::Float64      # Constant dynamic viscosity
    lambda_const::Float64   # Constant thermal conductivity
    T_min::Float64          # Minimum temperature valid for medium model
    T_max::Float64          # Maximum temperature valid for medium model
    T0::Float64             # Zero enthalpy temperature

    function SimpleIdealGasMediumData(;cp_const=nothing, R_gas=nothing, MM_const=nothing,
                                       eta_const=nothing, lambda_const=nothing, T_min=nothing, T_max=nothing, T0=nothing)
         cv_const = cp_const - R_gas
         new(cp_const, cv_const, R_gas, MM_const, eta_const, lambda_const,
             T_min, T_max, T0)
    end
end

"""
    medium = SimpleIdealGasMedium(; mediumName     = nothing,
                                    reference_p    = 101325.0,
                                    reference_T    = 298.15,
                                    p_default      = 101325.0,
                                    T_default      = 293.15,
                                    fluidConstants = nothing,
                                    data           = nothing)

Generate a `SimpleIdealGasMedium <: PureSubstance` medium object.
"""
mutable struct SimpleIdealGasMedium <: PureSubstance
    infos::FluidInfos
    fluidConstants::SVector{1,BasicFluidConstants}
    fluidLimits::FluidLimits
    data::SimpleIdealGasMediumData
end

"""
    state = SimpleIdealGasMediumState(Medium, p, T)

Generate a `SimpleIdealGasMediumState <: ThermodynamicState` object containing
pressure `p` [Pa] and temperature `T` [K] as thermodynamic states.
"""
mutable struct SimpleIdealGasMediumState <: ThermodynamicState
    Medium::SimpleIdealGasMedium
    p::Float64
    T::Float64
end

function SimpleIdealGasMedium(; mediumName=nothing,
                                reference_p=101325.0,
                                reference_T=298.15,
                                p_default=101325.0,
                                T_default=293.15,
                                fluidConstants=nothing,
                                data=nothing)

    infos = FluidInfos(mediumName           = mediumName,
                       substanceNames       = [mediumName],
                       extraPropertiesNames = fill("",0),
                       ThermoStates         = IndependentVariables_pT,
                       singleState          = false,
                       reducedX             = true,
                       fixedX               = false,
                       reference_p          = reference_p,
                       reference_T          = reference_T,
                       reference_X          = fill(1.0,1),
                       p_default            = p_default,
                       T_default            = T_default)

    fluidLimits = FluidLimits(TMIN = data.T_min,
                              TMAX = data.T_max)

    Medium = SimpleIdealGasMedium(infos, fill(fluidConstants,1), fluidLimits, data)

    infos.h_default  = specificEnthalpy(data, SimpleIdealGasMediumState(Medium,p_default,T_default))
    fluidLimits.HMIN = specificEnthalpy(data, SimpleIdealGasMediumState(Medium,p_default,data.T_min))
    fluidLimits.HMAX = specificEnthalpy(data, SimpleIdealGasMediumState(Medium,p_default,data.T_max))

    return Medium
end

setState_pTX(m::SimpleIdealGasMedium,p,T,X) = SimpleIdealGasMediumState(m,p,T)
setState_phX(m::SimpleIdealGasMedium,p,h,X) = SimpleIdealGasMediumState(m,p,m.data.T0+h/m.data.cp_const)
setState_psX(m::SimpleIdealGasMedium,p,s,X) = SimpleIdealGasMediumState(m,p,exp(s/m.data.cp_const + log(m.infos.reference_T)+m.infos.R_gas*log(p/m.infos.reference_p)))
setState_dTX(m::SimpleIdealGasMedium,d,T,X) = SimpleIdealGasMediumState(m,d*m.data.R_gas*T,T)
isenthalpicState(m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState, dp::Float64) = SimpleIdealGasMediumState(m, state.p+dp, state.T)

setState_pTX!(state::SimpleIdealGasMediumState,p,T,X) = begin state.p=p; state.T=T; nothing end
setState_phX!(state::SimpleIdealGasMediumState,p,h,X) = begin state.p=p; state.T=state.Medium.data.T0+h/state.Medium.data.cp_const; nothing end
setState_psX!(state::SimpleIdealGasMediumState,p,s,X) = begin state.p=p; state.T=exp(s/state.Medium.data.cp_const + log(state.Medium.infos.reference_T)); nothing end
setState_dTX!(state::SimpleIdealGasMediumState,d,T,X) = begin state.p=d*state.Medium.data.R_gas*T; state.T=T; nothing end
isenthalpicState!(state_b::SimpleIdealGasMediumState, state_a::SimpleIdealGasMediumState, dp::Float64) = begin state_b.p = state_a.p+dp; state_b.T = state_a.T; nothing end

specificEnthalpy(data::SimpleIdealGasMediumData, state::SimpleIdealGasMediumState)::Float64 = data.cp_const*(state.T - data.T0)
pressure(               m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = state.p
temperature(            m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = state.T
density(                m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = state.p/(m.data.R_gas*state.T)
specificEnthalpy(       m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = specificEnthalpy(m.data,state)
specificInternalEnergy( m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cp_const*(state.T-m.data.T0)-m.data.R_gas*state.T
specificEntropy(        m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cp_const*log(state.T/m.data.T0)-m.data.R_gas*log(state.p/m.infos.reference_p)
specificGibbsEnergy(    m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cp_const*(state.T-m.data.T0)-state.T*specificEntropy(m,state)
specificHelmholtzEnergy(m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = specificInternalEnergy(state)-state.T*specificEntropy(state)
dynamicViscosity(       m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.eta_const
thermalConductivity(    m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.lambda_const
specificHeatCapacityCp( m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cp_const
specificHeatCapacityCv( m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cv_const
isentropicExponent(     m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = m.data.cp_const/m.data.cv_const
velocityOfSound(        m::SimpleIdealGasMedium, state::SimpleIdealGasMediumState)::Float64 = sqrt(m.data.cp_const/m.data.cv_const*m.data.R_gas*state.T)

density_pT(      m::SimpleIdealGasMedium,p,T)::Float64 = p/(m.data.R_gas*T)
density_pT_der_1(m::SimpleIdealGasMedium,p,T)::Float64 = 0.0
density_pT_der_2(m::SimpleIdealGasMedium,p,T)::Float64 = 1.0/(m.data.R_gas*T)
density_pT_der_3(m::SimpleIdealGasMedium,p,T)::Float64 = -p/(m.data.R_gas*T*T)

specificInternalEnergy_T(      m::SimpleIdealGasMedium,T)::Float64 = m.data.cp_const*(T-m.data.T0)-m.data.R_gas*T
specificInternalEnergy_T_der_1(m::SimpleIdealGasMedium,T)::Float64 = 0.0
specificInternalEnergy_T_der_2(m::SimpleIdealGasMedium,T)::Float64 = m.data.cp_const - m.data.R_gas

specificEnthalpy_T(m::SimpleIdealGasMedium, T)::Float64 = m.data.cp_const*(T - m.data.T0)

function standardCharacteristics(m::SimpleIdealGasMedium)::Dict{AbstractString,Any}
    p_ref = m.infos.reference_p
    T     = collect( range(m.fluidLimits.TMIN, stop=min(1600.0, m.fluidLimits.TMAX), length=501) )
    p     = [0.5e5, 1.0e5, 2.0e5]
    nT    = length(T)
    np    = length(p)
    h     = zeros(nT)
    u     = zeros(nT)
    cp    = zeros(nT)
    cv    = zeros(nT)
    d     = zeros(nT,np)

    for i in 1:nT
        state = setState_pT(m,p_ref,T[i])
        h[i]  = specificEnthalpy(      m, state)
        u[i]  = specificInternalEnergy(m, state)
        cp[i] = specificHeatCapacityCp(m, state)
        cv[i] = specificHeatCapacityCv(m, state)
    end

    for j in 1:np
        for i in 1:nT
            d[i,j] = to_DensityDisplayUnit( density(SimpleIdealGasMediumState(m,p[j],T[i]) ) )
        end
    end

    mediumDict = Dict{AbstractString,Any}()
    mediumDict["T"]  = uconvert.(u"°C", T*1u"K")
    mediumDict["h"]  = h*1u"J/kg"
    mediumDict["u"]  = u*1u"J/kg"
    mediumDict["cp"] = cp*1u"J/(kg*K)"
    mediumDict["cv"] = cv*1u"J/(kg*K)"
    mediumDict["d(p=0.5 bar)"] = d[:,1]*1u"g/cm^3"
    mediumDict["d(p=1.0 bar)"] = d[:,2]*1u"g/cm^3"
    mediumDict["d(p=2.0 bar)"] = d[:,3]*1u"g/cm^3"
    return mediumDict
end

function standardPlot(m::SimpleIdealGasMedium; figure=1)
    mediumDict = standardCharacteristics(m)
    ModiaMath.plot(mediumDict, [("h", "u"), ("cp","cv"), ("d(p=0.5 bar)" ,
                                                          "d(p=1.0 bar)" ,
                                                          "d(p=2.0 bar)")], xAxis="T", heading=m.infos.mediumName, figure=figure)
end


text:
struct MbplsrRosa
    T::Matrix{Float64}
    P::Matrix{Float64}
    R::Matrix{Float64}
    W::Matrix{Float64}
    C::Matrix{Float64}
    TT::Vector{Float64}
    xmeans::Vector
    ymeans::Vector{Float64}
    weights::Vector{Float64}
    bl::Vector
end

"""
    rosaplsr(X, Y, weights = ones(size(X, 1)); nlv)
Multi-block PLSR with the ROSA algorithm (Liland et al. 2016).
* `X` : List (vector) of blocks (matrices) of X-data.
    Each component of the list is a block.
* `Y` : Y-data.
* `weights` : Weights of the observations (rows).
* `nlv` : Nb. latent variables (LVs) to consider.

The function has the following differences with the original
algorithm of Liland et al. (2016):
* Scores T are not normed to 1.
* Multivariate `Y` is allowed. In such a case, , the squared residuals are summed
    over the columns for finding the winning blocks the entered (therefore Y-columns
    should have the same scale).

Vector `weights` (row-weighting) is internally normalized to sum to 1.
See the help of `plskern` for details.

`X` and `Y` are internally centered.

## References

Liland, K.H., Næs, T., Indahl, U.G., 2016. ROSA—a fast extension of partial least
squares regression for multiblock data analysis. Journal of Chemometrics 30,
651–662. https://doi.org/10.1002/cem.2824
"""
function mbplsr_rosa(X, Y, weights = ones(size(X[1], 1)); nlv)
    nbl = length(X)
    zX = list(nbl, Matrix{Float64})
    @inbounds for i = 1:nbl
        zX[i] = copy(ensure_mat(X[i]))
    end
    mbplsr_rosa!(zX, copy(Y), weights; nlv = nlv)
end

function mbplsr_rosa!(X, Y, weights = ones(size(X[1], 1)); nlv)
    Y = ensure_mat(Y)
    n = size(X[1], 1)
    q = size(Y, 2)
    nbl = length(X)
    weights = mweights(weights)
    D = Diagonal(weights)
    xmeans = list(nbl, Vector{Float64})
    p = fill(0, nbl)
    @inbounds for i = 1:nbl
        p[i] = size(X[i], 2)
        xmeans[i] = colmeans(X[i], weights)
        center!(X[i], xmeans[i])
    end
    ymeans = colmeans(Y, weights)
    center!(Y, ymeans)
    # Pre-allocation
    W = similar(X[1], sum(p), nlv)
    P = copy(W)
    T = similar(X[1], n, nlv)
    TT = similar(X[1], nlv)
    C = similar(X[1], q, nlv)
    DY = similar(X[1], n, q)
    t   = similar(X[1], n)
    dt  = similar(X[1], n)
    c   = similar(X[1], q)
    zp_bl = list(nbl, Vector{Float64})
    zp = similar(X[1], sum(p))
    ssr = similar(X[1], nbl)
    W_bl = list(nbl, Array{Float64})
    w_bl = list(nbl, Vector{Float64})  # List of the weights "w" by block for a given "a"
    zT = similar(X[1], n, nbl)         # Matrix gathering the nbl scores for a given "a"
    bl = fill(0, nlv)
    Res = zeros(n, q, nbl)
    ### Start
    @inbounds for a = 1:nlv
        DY .= D * Y  # apply the metric on covariance
        @inbounds for i = 1:nbl
            XtY = X[i]' * DY
            if q == 1
                w_bl[i] = vec(XtY)
                w_bl[i] ./= sqrt(dot(w_bl[i], w_bl[i]))
            else
                w_bl[i] = svd!(XtY).U[:, 1]
            end
            zT[:, i] .= X[i] * w_bl[i]
        end
        # GS Orthogonalization of the scores
        if a > 1
            z = vcol(T, 1:(a - 1))
            zT .= zT .- z * inv(z' * (D * z)) * z' * (D * zT)
        end
        # Selection of the winner block
        @inbounds for i = 1:nbl
            t = vcol(zT, i)
            dt .= weights .* t
            tt = dot(t, dt)
            Res[:, :, i] .= Y .- (t * t') * DY / tt
        end
        ssr = vec(sum(Res.^2, dims = (1, 2)))
        opt = findmin(ssr)[2][1]
        bl[a] = opt
        # Outputs for winner
        t .= zT[:, opt]
        dt .= weights .* t
        tt = dot(t, dt)
        mul!(c, Y', dt)
        c ./= tt
        T[:, a] .= t
        TT[a] = tt
        C[:, a] .= c
        Y .= Res[:, :, opt]
        for i = 1:nbl ; zp_bl[i] = X[i]' * dt ; end
        zp .= reduce(vcat, zp_bl)
        P[:, a] .= zp / tt
        # Orthogonalization of the weights "w" by block
        zw = w_bl[opt]
        if (a > 1) && isassigned(W_bl, opt)
            z = W_bl[opt]
            zw .= zw .- z * (z' * zw)
        end
        zw ./= sqrt(dot(zw, zw))
        if !isassigned(W_bl, opt)
            W_bl[opt] = reshape(zw, :, 1)
        else
            W_bl[opt] = hcat(W_bl[opt], zw)
        end
        # Build the weights over the overall matrix
        z = zeros(nbl) ; z[opt] = 1
        W[:, a] .= reduce(vcat, z .* w_bl)
    end
    R = W * inv(P' * W)
    MbplsrRosa(T, P, R, W, C, TT, xmeans, ymeans, weights, bl)
end

"""
    transform(object::MbplsrRosa, X; nlv = nothing)
Compute LVs ("scores" T) from a fitted model.
* `object` : The maximal fitted model.
* `X` : A list (vector) of blocks (matrices) of X-data for which LVs are computed.
* `nlv` : Nb. LVs to consider. If nothing, it is the maximum nb. LVs.
"""
function transform(object::MbplsrRosa, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    nbl = length(object.xmeans)
    zX = list(nbl, Matrix{Float64})
    @inbounds for i = 1:nbl
        zX[i] = center(X[i], object.xmeans[i])
    end
    reduce(hcat, zX) * vcol(object.R, 1:nlv)
end

"""
    coef(object::MbplsrRosa; nlv = nothing)
Compute the X b-coefficients of a model fitted with `nlv` LVs.
* `object` : The maximal fitted model.
* `nlv` : Nb. LVs to consider. If nothing, it is the maximum nb. LVs.
"""
function coef(object::MbplsrRosa; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    zxmeans = reduce(vcat, object.xmeans)
    beta = object.C[:, 1:nlv]'
    B = vcol(object.R, 1:nlv) * beta
    int = object.ymeans' .- zxmeans' * B
    (B = B, int = int)
end

"""
    predict(object::MbplsrRosa, X; nlv = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : A list (vector) of X-data for which predictions are computed.
* `nlv` : Nb. LVs, or collection of nb. LVs, to consider.
    If nothing, it is the maximum nb. LVs.
"""
function predict(object::MbplsrRosa, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    le_nlv = length(nlv)
    zX = reduce(hcat, X)
    pred = list(le_nlv, Matrix{Float64})
    @inbounds for i = 1:le_nlv
        z = coef(object; nlv = nlv[i])
        pred[i] = z.int .+ zX * z.B
    end
    le_nlv == 1 ? pred = pred[1] : nothing
    (pred = pred,)
end

ext: jl


text:
#-------------------------------------------------------------------------------
# This function creates an instance of SNEASY+MAGICC.
#-------------------------------------------------------------------------------

# Load required packages.
using Mimi
using MimiMAGICC
using MimiSNEASY
using DataFrames
using CSVFiles

function create_sneasy_magiccch4(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2300, etminan_ch4_forcing::Bool=true)

    # ---------------------------------------------
    # Load and clean up necessary data.
    # ---------------------------------------------

    # Load RCP emissions and concentration scenario values (RCP options = "RCP26" and "RCP85").
    rcp_emissions      = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_emissions.csv"), skiplines_begin=36))
    rcp_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_concentrations.csv"), skiplines_begin=37))
    rcp_forcing        = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_midyear_radforcings.csv"), skiplines_begin=58))

    # Load MAGICC data for historical fossil carbon dioxide emissions values used to calculate pre-2000 non-methane ozone forcing.
    foss_hist_for_O₃   = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", "FOSS_HIST_magicc.csv")))

    # Find start and end year indices to crop RCP scenario data to correct model time horizon.
    rcp_indices = findall((in)(collect(start_year:end_year)), rcp_emissions.YEARS)

    # Set pre-industrial atmospheric CO₂, CH₄, and N₂O concentrations to RCP values in 1765.
    CO₂_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CO2][1]
    CH₄_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CH4][1]
    N₂O_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :N2O][1]

    # Calculate carbon dioxide emissions as well as aerosol and exogenous RCP radiative forcing scenarios.
    rcp_co2_emissions     = rcp_emissions.FossilCO2 .+ rcp_emissions.OtherCO2
    rcp_aerosol_forcing   = rcp_forcing.TOTAER_DIR_RF .+ rcp_forcing.CLOUD_TOT_RF
    rcp_exogenous_forcing = rcp_forcing.TOTAL_INCLVOLCANIC_RF .- rcp_forcing.CO2_RF .- rcp_forcing.CH4_RF .- rcp_forcing.TROPOZ_RF .- rcp_forcing.CH4OXSTRATH2O_RF .- rcp_aerosol_forcing

    # ------------------------------------------------------------
    # Initialize Mimi-SNEASY and add new CH₄ components to model.
    # ------------------------------------------------------------

    # Get an instance of Mimi-SNEASY.
    m = MimiSNEASY.get_model(start_year=start_year, end_year=end_year)

    # Remove old radiative forcing components.
    delete!(m, :rfco2)
    delete!(m, :radiativeforcing)

    # Add in new components.
    add_comp!(m, rf_total, before = :doeclim)
    add_comp!(m, rf_co2_etminan, before = :rf_total)
    add_comp!(m, MimiMAGICC.rf_ch4h2o, before = :rf_co2_etminan)
    add_comp!(m, MimiMAGICC.rf_o3, before = :rf_ch4h2o)
    add_comp!(m, total_co2_emissions, before= :rf_o3)
    add_comp!(m, MimiMAGICC.ch4_cycle, before = :total_co2_emissions)

    # Add in user-specified CH₄ radiative forcing component.
    # Note: If not using Etminan et al. equations, use original forcing equations from parent CH₄ model.
    if etminan_ch4_forcing == true
        add_comp!(m, rf_ch4_etminan, before = :rf_ch4h2o)
    else
        add_comp!(m, MimiMAGICC.rf_ch4, before = :rf_ch4h2o)
    end

    # ---------------------------------------------
    # Set component parameters.
    # ---------------------------------------------

    # ---- Common parameters ----
    Mimi.set_external_param!(m, :CH₄_0, CH₄_0)
    Mimi.set_external_param!(m, :N₂O_0, N₂O_0)
    Mimi.set_external_param!(m, :N₂O, rcp_concentrations.N2O[rcp_indices], param_dims=[:time])
    Mimi.set_external_param!(m, :CO_emissions, rcp_emissions.CO[rcp_indices], param_dims=[:time])
    Mimi.set_external_param!(m, :NMVOC_emissions, rcp_emissions.NMVOC[rcp_indices], param_dims=[:time])
    Mimi.set_external_param!(m, :NOx_emissions, rcp_emissions.NOx[rcp_indices], param_dims=[:time])

    # ---- Carbon Cycle ---- #
    update_param!(m, :atmco20, CO₂_0)

    # ---- Methane Cycle ---- #
    set_param!(m, :ch4_cycle, :BBCH4, 2.78)
    set_param!(m, :ch4_cycle, :ANOX, 0.0042)
    set_param!(m, :ch4_cycle, :ACO, -0.000105)
    set_param!(m, :ch4_cycle, :AVOC, -0.000315)
    set_param!(m, :ch4_cycle, :TAUINIT, 7.73)
    set_param!(m, :ch4_cycle, :SCH4, -0.32)
    set_param!(m, :ch4_cycle, :TAUSOIL, 160.)
    set_param!(m, :ch4_cycle, :TAUSTRAT, 120.)
    set_param!(m, :ch4_cycle, :GAM, -1.0)
    set_param!(m, :ch4_cycle, :CH4_natural, 266.5)
    set_param!(m, :ch4_cycle, :fffrac, 0.18)
    connect_param!(m, :ch4_cycle, :CH₄_0, :CH₄_0)
    set_param!(m, :ch4_cycle, :CH4_emissions, rcp_emissions.CH4[rcp_indices])
    connect_param!(m, :ch4_cycle, :NOx_emissions, :NOx_emissions)
    connect_param!(m, :ch4_cycle, :CO_emissions, :CO_emissions)
    connect_param!(m, :ch4_cycle, :NMVOC_emissions, :NMVOC_emissions)

    # ---- Total Carbon Dioxide Emissions ---- #
    set_param!(m, :total_co2_emissions, :exogenous_CO₂_emissions, rcp_co2_emissions[rcp_indices])

    # ---- Methane Radiative Forcing ---- #
    if etminan_ch4_forcing == true
        connect_param!(m, :rf_ch4_etminan, :CH₄_0, :CH₄_0)
        connect_param!(m, :rf_ch4_etminan, :N₂O_0, :N₂O_0)
        set_param!(m, :rf_ch4_etminan, :scale_CH₄, 1.0)
        set_param!(m, :rf_ch4_etminan, :a₃, -1.3e-6)
        set_param!(m, :rf_ch4_etminan, :b₃, -8.2e-6)
        connect_param!(m, :rf_ch4_etminan, :N₂O, :N₂O)
    else
        connect_param!(m, :rf_ch4, :N₂O_0, :N₂O_0)
        connect_param!(m, :rf_ch4, :CH₄_0, :CH₄_0)
        set_param!(m, :rf_ch4, :scale_CH₄, 1.0)
    end

    # ---- Straospheric Water Vapor From Oxidized Methane Radiative Forcing ---- #
    connect_param!(m, :rf_ch4h2o, :CH₄_0, :CH₄_0)
    set_param!(m, :rf_ch4h2o, :STRATH2O, 0.15)

    # ---- Tropospheric Ozone Radiative Forcing ---- #
    connect_param!(m, :rf_o3, :CH₄_0, :CH₄_0)
    set_param!(m, :rf_o3, :OZ00CH4, 0.161)
    set_param!(m, :rf_o3, :OZNOX, 0.125)
    set_param!(m, :rf_o3, :OZCO, 0.0011)
    set_param!(m, :rf_o3, :OZVOC, 0.0033)
    connect_param!(m, :rf_o3, :NOx_emissions, :NOx_emissions)
    connect_param!(m, :rf_o3, :CO_emissions, :CO_emissions)
    connect_param!(m, :rf_o3, :NMVOC_emissions, :NMVOC_emissions)
    set_param!(m, :rf_o3, :FOSSHIST, foss_hist_for_O₃.EFOSS[rcp_indices])
    set_param!(m, :rf_o3, :TROZSENS, 0.042)
    set_param!(m, :rf_o3, :OZCH4, 5.0)

    # ---- Carbon Dioxide Radiative Forcing ---- #
    set_param!(m, :rf_co2_etminan, :a₁, -2.4e-7)
    set_param!(m, :rf_co2_etminan, :b₁, 7.2e-4)
    set_param!(m, :rf_co2_etminan, :c₁, -2.1e-4)
    set_param!(m, :rf_co2_etminan, :CO₂_0, CO₂_0)
    connect_param!(m, :rf_co2_etminan, :N₂O_0, :N₂O_0)
    connect_param!(m, :rf_co2_etminan, :N₂O, :N₂O)
    set_param!(m, :rf_co2_etminan, :rf_scale_CO₂, co2_rf_scale(3.7, CO₂_0, N₂O_0))

    # ---- Total Radiative Forcing ---- #
    set_param!(m, :rf_total, :α, 1.0)
    # TODO It would be nice if `rf_aerosol` would not be a model
    # parameter at this point
    connect_param!(m, :rf_total, :rf_aerosol, :rf_aerosol)
    update_param!(m, :rf_aerosol, rcp_aerosol_forcing[rcp_indices])
    set_param!(m, :rf_total, :rf_exogenous, rcp_exogenous_forcing[rcp_indices])

    # ----------------------------------------------------------
    # Create connections between Mimi SNEASY+Hector components.
    # ----------------------------------------------------------
    connect_param!(m, :doeclim,             :forcing,             :rf_total,            :total_forcing)
    connect_param!(m, :ch4_cycle,           :temperature,         :doeclim,             :temp)
    connect_param!(m, :total_co2_emissions, :oxidized_CH₄_to_CO₂, :ch4_cycle,           :emeth)
    connect_param!(m, :ccm,                 :CO2_emissions,       :total_co2_emissions, :total_CO₂_emissions)
    connect_param!(m, :ccm,                 :temp,                :doeclim,             :temp)
    connect_param!(m, :rf_co2_etminan,      :CO₂,                 :ccm,                 :atmco2)
    connect_param!(m, :rf_o3,               :CH₄,                 :ch4_cycle,           :CH₄)
    connect_param!(m, :rf_ch4h2o,           :CH₄,                 :ch4_cycle,           :CH₄)

    # Create different connections if using updated Etminan et al. CH₄ forcing equations.
    if etminan_ch4_forcing == true
        connect_param!(m, :rf_ch4_etminan, :CH₄,    :ch4_cycle,      :CH₄)
        connect_param!(m, :rf_total,       :rf_CH₄, :rf_ch4_etminan, :rf_CH₄)
    else
        connect_param!(m, :rf_ch4,         :CH₄,    :ch4_cycle,      :CH₄)
        connect_param!(m, :rf_total,       :rf_CH₄, :rf_ch4,         :QMeth)
    end

    connect_param!(m, :rf_total, :rf_O₃,      :rf_o3,          :rf_O₃)
    connect_param!(m, :rf_total, :rf_CH₄_H₂O, :rf_ch4h2o,      :QCH4H2O)
    connect_param!(m, :rf_total, :rf_CO₂,     :rf_co2_etminan, :rf_CO₂)

    # Return constructed model.
    return m
end


text:
words = ["Enjoy", "Rosetta", "Code"]

function sleepprint(s)
    sleep(rand())
    println(s)
end

@sync for word in words
    @async sleepprint(word)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using SpectralKit
using Test, DocStringExtensions, StaticArrays, BenchmarkTools, Sobol
using ForwardDiff: derivative, gradient

include("utilities.jl")

include("test_generic_api.jl")
include("test_chebyshev.jl")
include("test_univariate.jl")
include("test_smolyak.jl")


text:
using ITensors

# Text visualization, prints to screen
using ITensorUnicodePlots

# Visualization in an interactive window
using ITensorGLMakie

using LinearAlgebra

ITensorVisualizationBase.set_backend!("UnicodePlots")

println("
#######################################################
# Tutorial 1                                          #
#                                                     #
# 1-site state basics                                 #
#######################################################
")

i = Index(2)

#
# Making states
#

# Make an "up" state (+1 eigenvalue of `Z`),
# denoted as `|Z+⟩`.

# Define through setting elements
Zp = ITensor(i)
Zp[i => 1] = 1.0

# Alternative syntax, construct from a Vector
Zp = ITensor([1 0], i)

# Make an "down" state (-1 eigenvalue of `Z`),
# denoted as `|Z-⟩`.

# Define through setting elements
Zm = ITensor(i)
Zm[i => 2] = 1.0

# Alternative syntax, construct from a Vector
Zm = ITensor([0 1], i)

@visualize Zm

# Can do algebra, inner products, etc:
@show (Zp + Zm) / √2
@show dag(Zp) * Zm


text:
# Logger.jl

using Dates

const INFO     = 1
const WARNING  = 2
const ERROR    = 3

struct Logger
    filename   # log file name
    level      # minimum level acceptable to be logged
    handle     # file handle
end

Logger(filename, level) = Logger(filename, level, open(filename, "w"))

# info logging
function info!(logger::Logger, args...)
    if logger.level <= INFO
        let io = logger.handle
            print(io, trunc(now(), Dates.Second), " [INFO] ")
            for (idx, arg) in enumerate(args)
                idx > 0 && print(io, " ")
                print(io, arg)
            end
            println(io)
            flush(io)
        end
    end
end

# testing
info_logger = Logger("/tmp/info.log", INFO)
info!(info_logger, "hello", 123)

readlines("/tmp/info.log")
#=
julia> readlines("/tmp/info.log")
1-element Array{String,1}:
 "2019-09-07T22:22:46 [INFO]  hello 123"
=#

error_logger = Logger("/tmp/error.log", ERROR)
info!(error_logger, "hello", 123)

readlines("/tmp/error.log")
#=
julia> readlines("/tmp/error.log")
0-element Array{String,1}
=#

# now I'm tempted to just finish the project by duplicating the code
# and make slight changes.
# => Violation of DRY principle!
function warning!(logger::Logger, args...)
    if logger.level <= WARNING
        let io = logger.handle
            print(io, trunc(now(), Dates.Second), " [WARNING] ")
            for (idx, arg) in enumerate(args)
                idx > 0 && print(io, " ")
                print(io, arg)
            end
            println(io)
            flush(io)
        end
    end
end

# code gen
for level in (:info, :warning, :error)
    lower_level_str = String(level)
    upper_level_str = uppercase(lower_level_str)
    upper_level_sym = Symbol(upper_level_str)

    fn = Symbol(lower_level_str * "!")
    label = " [" * upper_level_str * "] "

    @eval function $fn(logger::Logger, args...)
        if logger.level <= $upper_level_sym
            let io = logger.handle
                print(io, trunc(now(), Dates.Second), $label)
                for (idx, arg) in enumerate(args)
                    idx > 0 && print(io, " ")
                    print(io, arg)
                end
                println(io)
                flush(io)
            end
        end
    end
end

info_logger = Logger("/tmp/info.log", INFO)
info!(info_logger,    "hello", 123)
warning!(info_logger, "hello", 123)
error!(info_logger,   "hello", 123)

warning_logger = Logger("/tmp/warning.log", WARNING)
info!(warning_logger,    "hello", 123)
warning!(warning_logger, "hello", 123)
error!(warning_logger,   "hello", 123)

error_logger = Logger("/tmp/error.log", ERROR)
info!(error_logger,    "hello", 123)
warning!(error_logger, "hello", 123)
error!(error_logger,   "hello", 123)

readlines("/tmp/info.log")
readlines("/tmp/warning.log")
readlines("/tmp/error.log")
#=
julia> readlines("/tmp/info.log")
3-element Array{String,1}:
 "2019-12-08T17:03:28 [INFO]  hello 123"
 "2019-12-08T17:03:28 [WARNING]  hello 123"
 "2019-12-08T17:03:28 [ERROR]  hello 123"

julia> readlines("/tmp/warning.log")
2-element Array{String,1}:
 "2019-12-08T17:03:29 [WARNING]  hello 123"
 "2019-12-08T17:03:29 [ERROR]  hello 123"

julia> readlines("/tmp/error.log")
1-element Array{String,1}:
 "2019-12-08T17:03:30 [ERROR]  hello 123"
=#

# debugging - https://github.com/timholy/CodeTracking.jl

methods(error!)
methods(error!) |> first

using Revise, CodeTracking
methods(error!) |> first |> definition
#=
julia> methods(error!)
# 1 method for generic function "error!":
[1] error!(logger::Logger, args...) in Main at REPL[70]:8

julia> methods(error!) |> first
error!(logger::Logger, args...) in Main at REPL[70]:8

julia> using Revise, CodeTracking

julia> methods(error!) |> first |> definition
:(function error!(logger::Logger, args...)
      #= REPL[24]:8 =#
      if logger.level <= ERROR
          #= REPL[24]:9 =#
          let io = logger.handle
              #= REPL[24]:10 =#
              print(io, trunc(now(), Dates.Second), " [ERROR] ")
              #= REPL[24]:11 =#
              for (idx, arg) = enumerate(args)
                  #= REPL[24]:12 =#
                  idx > 0 && print(io, " ")
                  #= REPL[24]:13 =#
                  print(io, arg)
              end
              #= REPL[24]:15 =#
              println(io)
              #= REPL[24]:16 =#
              flush(io)
          end
      end
  end)
=#

# MacroTools.prewalk

using MacroTools
MacroTools.prewalk(rmlines, definition(first(methods(error!))))
#=
julia> using MacroTools

julia> MacroTools.prewalk(rmlines, definition(first(methods(error!))))
:(function error!(logger::Logger, args...)
      if logger.level <= ERROR
          let io = logger.handle
              print(io, trunc(now(), Dates.Second), " [ERROR] ")
              for (idx, arg) = enumerate(args)
                  idx > 0 && print(io, " ")
                  print(io, arg)
              end
              println(io)
              flush(io)
          end
      end
  end)
=#

# using Lazy and chain with rmlines

using Lazy
@>> methods(error!) first definition MacroTools.prewalk(rmlines)
#=
using Lazy

julia> @>> methods(error!) first definition MacroTools.prewalk(rmlines)
:(function error!(logger::Logger, args...)
      if logger.level <= ERROR
          let io = logger.handle
              print(io, trunc(now(), Dates.Second), " [ERROR] ")
              for (idx, arg) = enumerate(args)
                  idx > 0 && print(io, " ")
                  print(io, arg)
              end
              println(io)
              flush(io)
          end
      end
  end)
=#

function logme!(level, label, logger::Logger, args...)
    if logger.level <= level
        let io = logger.handle
            print(io, trunc(now(), Dates.Second), label)
            for (idx, arg) in enumerate(args)
                idx > 0 && print(io, " ")
                print(io, arg)
            end
            println(io)
            flush(io)
        end
    end
end

info!(logger::Logger, msg...)    = logme!(INFO,    " [INFO] ",    logger, msg...)
warning!(logger::Logger, msg...) = logme!(WARNING, " [WARNING] ", logger, msg...)
error!(logger::Logger, msg...)   = logme!(ERROR,   " [ERROR] ",   logger, msg...)

#----- using closuure -----
function make_log_func(level, label)
    (logger::Logger, args...) -> begin
        if logger.level <= level
            let io = logger.handle
                print(io, trunc(now(), Dates.Second), " [", label, "] ")
                for (idx, arg) in enumerate(args)
                    idx > 0 && print(io, " ")
                    print(io, arg)
                end
                println(io)
                flush(io)
            end
        end
    end
end

info!    = make_log_func(INFO, "INFO")
warning! = make_log_func(WARNING, "WARNING")
error!   = make_log_func(ERROR, "ERROR")


text:
module Simulate_Pendulum

using  Modia3D
import Modia3D.ModiaMath

# visualization material
visuMaterial1 = Modia3D.Material(color="LightBlue", transparency=0.5)
visuMaterial2 = Modia3D.Material(color="Red")

@assembly Pendulum(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx) begin
   # world is the reference point and parent Object3D of all other Object3D
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.4*Lx))

   # Pendulum
   # the beam is made of Aluminium and is visualized with blue color
   beam     = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), "Aluminium", visuMaterial1))
   # the bearing is a helping Object3D, its parent is the beam
   bearing  = Modia3D.Object3D(beam; r=[-Lx/2, 0.0, 0.0])
   # for visualizing the revolute joint
   cylinder = Modia3D.Object3D(bearing, Modia3D.Cylinder(Ly/2,1.2*Ly; material=visuMaterial2))

   # the revolute joint connects the world with the bearing, which belongs to the beam
   revolute    = Modia3D.Revolute(world, bearing)
end
# a Pendulum object is created, where SceneOptions are set (in this case visualization properties of frames)
pendulum = Pendulum(Lx=1.6, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3))

# a simulation model is created, simulated and potted
model  = Modia3D.SimulationModel(pendulum)
result = ModiaMath.simulate!(model, stopTime=5.0, tolerance=1e-6,log=true)
ModiaMath.plot(result, ["revolute.phi", "revolute.w", "revolute.a"])

println("... success of Simulate_Pendulum.jl!")
end


text:
# Benchmarks for projection time of FiniteCone vs. PolyhedralCone

using BenchmarkLite
using Distributions, TailRiskScenGen

# FiniteCone Benchmark test definition

type ConeProject <: Proc
    dist::ContinuousMultivariateDistribution
    K::Cone
    function ConeProject(dist::ContinuousMultivariateDistribution, K::Cone)
        length(dist)==length(K) || throw(ArgumentError("Distribution and cone must be defined in same dimension"))
        new(dist, K)
    end
end

function Base.string(proc::ConeProject)
    dim = length(proc.K)
    ty = typeof(proc.K)
    "$(ty): Dim=$(dim)"
end
Base.length(proc,cfg) = cfg
Base.isvalid(proc::ConeProject, cfg) = true
function Base.start(proc::ConeProject, cfg)
    rand(proc.dist, cfg)
end

function Base.run(proc::ConeProject, cfg, s)
    for i in 1:cfg
        project(proc.K, s[:,i])
    end
end
Base.done(proc::ConeProject, cfg, s) = nothing

# Run test
dims = [2, 5, 10, 20]
procs = Proc[ConeProject(MvNormal(eye(40)), FiniteCone(eye(40))), ConeProject(MvNormal(eye(40)), PolyhedralCone(eye(40)))]
cfgs = [100, 1000, 2000]

rtable = run(procs, cfgs)
show(rtable)


text:
function basic_linalg_problem(;Ns=1000,Nc=100)

  A = Diagonal(2*ones(Ns))
  C = Diagonal(ones(Nc))
  B2 = zeros(Nc,Ns)
  for j in 1:min(Nc,Ns)
      B2[j,j] = 1.0
  end
  B1t = B2'
  rhs1v = ones(Ns)
  rhs2v = 2*ones(Nc)

  Abig = [A B1t;B2 C]
  rhsbig = [rhs1v;rhs2v]
  solex = similar(rhsbig)
  solex .= Abig\rhsbig;

  return A,B2,B1t,C,rhs1v,rhs2v, solex

end


text:
include("base.jl")

function analyze(widget::QStyle, ctx::PlotContext)
    return analyze(widget.data, ctx)
end

function renderWidget!(widget::QStyle, ctx::PlotContext)
    for (k, v) in widget.overloads
        ctx.styles[k] = v
    end
    return renderWidget!(widget.data, ctx)
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
module ExportMulti

using YAML
using ACE1
using ACE1: PIBasis, PIBasisFcn, PIPotential
using ACE1.OrthPolys: TransformedPolys
using ACE1: rand_radial, cutoff, numz, ZList
using JuLIP: energy, bulk, i2z, z2i, chemical_symbol

function export_ACE(fname, IP)
    #decomposing into V1, V2, V3 (One body, two body and ACE bases)
    V1 = IP.components[1]
    V2 = IP.components[2]
    V3 = IP.components[3]

    if hasproperty(V2, :basis)
        species = collect(string.(chemical_symbol.(V2.basis.zlist.list.data)))
    else hasproperty(V2, :Vout)
        species = collect(string.(chemical_symbol.(V2.Vout.basis.zlist.list.data)))
    end

    species_dict = Dict(zip(collect(0:length(species)-1), species))
    reversed_species_dict = Dict(zip(species, collect(0:length(species)-1)))

    data = Dict()

    data["deltaSplineBins"] = 0.001 #" none

    elements = Vector(undef, length(species))
    E0 = zeros(3)

    for (index, element) in species_dict
        E0[index+1] = V1(Symbol(element))
        elements[index+1] = element
    end

    # grabbing the elements key and E0 key from the onebody (V1)
    data["elements"] = elements
    data["E0"] = E0

    if hasproperty(V2, :basis)
        polypairpot = export_polypairpot(V2, reversed_species_dict)
    else hasproperty(V2, :Vin)
        polypairpot = export_polypairpot(V2.Vout, reversed_species_dict)
        reppot = export_reppot(V2, reversed_species_dict)
        data["reppot"] = reppot
    end

    data["polypairpot"] = polypairpot
    #creating "data" dict where we'll store everything

    embeddings, bonds = export_radial_basis(V3, species_dict)
    data["embeddings"] = embeddings
    data["bonds"] = bonds

    functions, lmax = export_ACE_functions(V3, species, reversed_species_dict)
    data["functions"] = functions
    data["lmax"] = lmax

    YAML.write_file(fname, data)
end

function export_reppot(Vrep, reversed_species_dict)
    reppot = Dict("coefficients" => Dict())

    zlist_dict = Dict(zip(1:length(Vrep.Vout.basis.zlist.list), [string(chemical_symbol(z)) for z in Vrep.Vout.basis.zlist.list]))

    for (index1, element1) in zlist_dict
        for (index2, element2) in zlist_dict
            pair = [reversed_species_dict[element1], reversed_species_dict[element2]]
            coefficients = Dict( "A" => Vrep.Vin[index1, index2].A,
                                "B" => Vrep.Vin[index1, index2].B,
                                "e0" => Vrep.Vin[index1, index2].e0,
                                "ri" => Vrep.Vin[index1, index2].ri)
            reppot["coefficients"][pair] = coefficients
        end
    end

    return reppot
end

function export_polypairpot(V2, reversed_species_dict)
    Pr = V2.basis.J

    p = Pr.trans.p
    r0 = Pr.trans.r0
    xr = Pr.J.tr
    xl = Pr.J.tl
    pr = Pr.J.pr
    pl = Pr.J.pl
    rcut = cutoff(Pr)
    maxn = length(Pr)

    if length(keys(reversed_species_dict)) == 1
        num_coeffs = length(V2.coeffs)
    else
        num_coeffs = vcat(V2.basis.bidx0...)[2]
    end

    zlist_dict = Dict(zip(1:length(V2.basis.zlist.list), [string(chemical_symbol(z)) for z in V2.basis.zlist.list]))

    polypairpot = Dict( "p" => p,
                        "r0" => r0,
                        "xr" => xr,
                        "xl" => xl,
                        "pr" => pr,
                        "pl" => pl,
                        "rcut" => rcut,
                        "maxn"=> maxn,
                        "recursion_coefficients" => Dict("A" => [Pr.J.A[i] for i in 1:maxn],
                                                         "B" => [Pr.J.B[i] for i in 1:maxn],
                                                         "C" => [Pr.J.C[i] for i in 1:maxn],),
                        "coefficients" => Dict())

    for (index1, element1) in zlist_dict
        for (index2, element2) in zlist_dict
            pair = [reversed_species_dict[element1], reversed_species_dict[element2]]
            ind = V2.basis.bidx0[index1, index2]
            polypairpot["coefficients"][pair] = V2.coeffs[ind+1:ind+num_coeffs]
        end
    end

    return polypairpot
end

function export_radial_basis(V3, species_dict)
    #grabbing the transform and basis
    transbasis = V3.pibasis.basis1p.J
    Pr = V3.pibasis.basis1p

    #grabbing all the required params
    p = transbasis.trans.p
    r0 = transbasis.trans.r0
    xr = transbasis.J.tr
    xl = transbasis.J.tl
    pr = transbasis.J.pr
    pl = transbasis.J.pl
    rcut = cutoff(Pr)
    maxn = length(V3.pibasis.basis1p.J.J.A)

    #guessing "radbasname" is that just "polypairpots"
    radbasename = "ACE1.jl.base"

    embeddings = Dict()

    for species_ind1 in sort(collect(keys(species_dict)))
        embeddings[species_ind1] = Dict("ndensity" => 1,
                    "FS_parameters" => [1.0, 1.0],
                    "npoti" => "FinnisSinclairShiftedScaled",
                    "drho_core_cutoff" => 1.000000000000000000,
                    "rho_core_cutoff" => 100000.000000000000000000)
    end

    bonds = Dict()
    #this does not respect the coefficient decompositions required per pair
    #need to figure out how to get the right coeffs per pair
    for species_ind1 in sort(collect(keys(species_dict)))
        for species_ind2 in sort(collect(keys(species_dict)))
            pair = [species_ind1, species_ind2]
            bonds[pair] = Dict("p" => p,
                "r0" => r0,
                "xl" => xl,
                "xr" => xr,
                "pr" => pr,
                "pl" => pl,
                "rcut" => rcut,
                "radbasename" => radbasename,
                "maxn" => maxn,
                "recursion_coefficients" => Dict("A" => [Pr.J.J.A[i] for i in 1:maxn],
                                                 "B" => [Pr.J.J.B[i] for i in 1:maxn],
                                                 "C" => [Pr.J.J.C[i] for i in 1:maxn],))
        end
    end

    return embeddings, bonds
end

function export_ACE_functions(V3, species, reversed_species_dict)
    functions = Dict()
    lmax = 0

    for i in 1:length(V3.pibasis.inner)
        sel_bgroups = []
        inner = V3.pibasis.inner[i]
        z0 = V3.pibasis.inner[i].z0
        coeffs = V3.coeffs[i]
        groups = _basis_groups(inner, coeffs)
        for group in groups
            for (m, c) in zip(group["M"], group["C"])
                c_ace = c / (4*π)^(group["ord"]/2)
                #@show length(c_ace)
                ndensity = 1
                push!(sel_bgroups, Dict("rank" => group["ord"],
                            "mu0" => reversed_species_dict[string(chemical_symbol(group["z0"]))],
                            "ndensity" => ndensity,
                            "ns" => group["n"],
                            "ls" => group["l"],
                            "mus" => [reversed_species_dict[i] for i in string.(chemical_symbol.(group["zs"]))],
                            "ctildes" => [c_ace],
                            "ms_combs" => m,
                            "num_ms_combs" => length([c_ace])))
                if maximum(group["l"]) > lmax
                    lmax = maximum(group["l"])
                end
            end
        end
        functions[reversed_species_dict[string(chemical_symbol(z0))]] = sel_bgroups
    end

    return functions, lmax
end

function _basis_groups(inner, coeffs)
    ## grouping the basis functions
    NLZZ = []
    M = []
    C = []
    for b in keys(inner.b2iAA)
       if coeffs[ inner.b2iAA[b] ] != 0
          push!(NLZZ, ( [b1.n for b1 in b.oneps], [b1.l for b1 in b.oneps], [b1.z for b1 in b.oneps], b.z0))
          push!(M, [b1.m for b1 in b.oneps])
          push!(C, coeffs[ inner.b2iAA[b] ])
       end
    end
    ords = length.(M)
    perm = sortperm(ords)
    NLZZ = NLZZ[perm]
    M = M[perm]
    C = C[perm]
    @assert issorted(length.(M))
    bgrps = []
    alldone = fill(false, length(NLZZ))
    for i = 1:length(NLZZ)
       if alldone[i]; continue; end
       nlzz = NLZZ[i]
       Inl = findall(NLZZ .== Ref(nlzz))
       alldone[Inl] .= true
       Mnl = M[Inl]
       Cnl = C[Inl]
       pnl = sortperm(Mnl)
       Mnl = Mnl[pnl]
       Cnl = Cnl[pnl]
       order = length(nlzz[1])
       push!(bgrps, Dict("n" => nlzz[1], "l" => nlzz[2], "z0" => nlzz[4], "zs" => nlzz[3],
                         "M" => Mnl, "C" => Cnl, "ord" => order)) #correct?
    end
    return bgrps
end

end


text:
export pairwise_identical, missing_data, geno_freqtable, allele_freqtable

#TODO update docs (API/dataexploration) and tests
"""
    missing_data(data::PopData; by::String = "sample")
Get missing genotype information in a `PopData`. Specify a mode of operation
to return a DataFrame corresponding with that missing information.

#### Modes
- "sample" - returns a count and list of missing loci per individual (default)
- "pop" - returns a count of missing genotypes per population
- "locus" - returns a count of missing genotypes per locus
- "full" - returns a count of missing genotypes per locus per population

### Example:
```
missing_data(@gulfsharks, by = "pop")
```
"""
@inline function missing_data(data::PopData; by::String = "sample")
    if by ∈ ["sample", "individual"]
        DataFrames.combine(
            DataFrames.groupby(data.loci, :name),
            :genotype => (i -> count(ismissing, i)) => :missing
        )
    elseif by ∈ ["pop", "population"]
                DataFrames.combine(
            DataFrames.groupby(data.loci, :population),
            :genotype => (i -> count(ismissing, i)) => :missing
        )
    elseif by ∈ ["locus", "loci"]
        DataFrames.combine(
            DataFrames.groupby(data.loci, :locus),
            :genotype => (i -> count(ismissing, i)) => :missing
        )

    elseif by ∈ ["detailed", "full"]
        DataFrames.combine(
            DataFrames.groupby(data.loci, [:locus, :population]),
            :genotype => (i -> count(ismissing, i)) => :missing
        )
    else
        @error "Mode \"$by\" not recognized. Please specify one of: sample, pop, locus, or full"
        missing_data(data)
    end
end

"""
    pairwise_identical(data::PopData)
Return a table of the percent of identical genotypes at each locus between all pairs of individuals.

### Example:
```
julia> cats = @nancycats ;

julia> pairwise_identical(cats)
27966×4 DataFrame
   Row │ sample_1  sample_2  identical  n
       │ String    String    Float64    Int64
───────┼──────────────────────────────────────
     1 │ N215      N216           0.5       8
     2 │ N215      N217           0.25      8
     3 │ N215      N218           0.38      8
     4 │ N215      N219           0.38      8
   ⋮   │    ⋮         ⋮          ⋮        ⋮
 27963 │ N297      N290           0.29      7
 27964 │ N281      N289           0.25      8
 27965 │ N281      N290           0.43      7
 27966 │ N289      N290           0.14      7
                            27958 rows omitted
```

"""
function pairwise_identical(data::PopData)
    sample_names = collect(samples(data))
    pairwise_identical(data, sample_names)
end

"""
    pairwise_identical(data::PopData, sample_names::Vector{String})
Return a table of the percent of identical genotypes at each locus
between all pairs of provided `sample_names`.

### Example:
```
julia> cats = @nancycats ;

julia> interesting_cats = samples(cats)[1:5]
5-element Array{String,1}:
 "N215"
 "N216"
 "N217"
 "N218"
 "N219"

julia> pairwise_identical(cats, interesting_cats)
10×4 DataFrame
 Row │ sample_1  sample_2  identical  n
     │ String    String    Float64    Int64
─────┼──────────────────────────────────────
   1 │ N215      N216           0.5       8
   2 │ N215      N217           0.25      8
   3 │ N215      N218           0.38      8
   4 │ N215      N219           0.38      8
   5 │ N216      N217           0.12      8
   6 │ N216      N218           0.25      8
   7 │ N216      N219           0.38      8
   8 │ N217      N218           0.0       9
   9 │ N217      N219           0.11      9
  10 │ N218      N219           0.33      9
```
"""
function pairwise_identical(data::PopData, sample_names::Vector{String})
    errs = ""
    all_samples = samples(data)
    if sample_names != all_samples
        [errs *= "\n  $i" for i in sample_names if i ∉ all_samples]
        errs != "" && error("Samples not found in the PopData: " * errs)
    end
    sample_pairs = pairwise_pairs(sample_names)
    n = length(sample_pairs)
    perc_ident_vec = Vector{Float64}(undef, n)
    n_vec = Vector{Int}(undef, n)
    popdata_idx = groupby(data.loci, :name)
    idx = 0
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    @inbounds @sync for i in 1:length(sample_pairs)
        Base.Threads.@spawn begin
            @inbounds geno_1 = popdata_idx[(sample_pairs[i][1],)].genotype
            @inbounds geno_2 = popdata_idx[(sample_pairs[i][2],)].genotype
            len_1 = nonmissing(geno_1)
            len_2 = nonmissing(geno_2)
            shared_geno = minimum([len_1, len_2])
                shared_geno = minimum([len_1, len_2])
                shared = sum(skipmissing(geno_1 .== geno_2))
                @inbounds perc_ident_vec[i] = round(shared/shared_geno, digits = 2)
                @inbounds n_vec[i] = shared_geno
            next!(p)
        end
    end
    DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :identical => perc_ident_vec, :n => n_vec)
end

"""
    geno_freqtable(data::PopData; by::String = "global")
Return a table of the observed `global` (default) or `population` genotype frequencies in a PopData object.

### Example:
```
julia> cats = @nancycats ;

julia> geno_freqtable(cats)

341×4 DataFrame
 Row │ locus   genotype    count  frequency
     │ String  Tuple…      Int64  Float64
─────┼───────────────────────────────────────
   1 │ fca8    (135, 143)     16  0.0737327
   2 │ fca8    (133, 135)      9  0.0414747
   3 │ fca8    (135, 135)     23  0.105991
   4 │ fca8    (137, 143)      8  0.0368664
  ⋮  │   ⋮         ⋮         ⋮        ⋮
 338 │ fca37   (206, 220)      1  0.00421941
 339 │ fca37   (208, 218)      1  0.00421941
 340 │ fca37   (184, 184)      3  0.0126582
 341 │ fca37   (208, 210)      3  0.0126582
                             333 rows omitted

julia> geno_freqtable(cats, by = "population")
1094×5 DataFrame
  Row │ locus   population  genotype    count  frequency
      │ String  String      Tuple…      Int64  Float64
──────┼──────────────────────────────────────────────────
    1 │ fca8    1           (135, 143)      3  0.375
    2 │ fca8    1           (133, 135)      2  0.25
    3 │ fca8    1           (135, 135)      2  0.25
    4 │ fca8    1           (137, 143)      1  0.125
  ⋮   │   ⋮         ⋮           ⋮         ⋮        ⋮
 1091 │ fca37   17          (208, 208)     10  0.769231
 1092 │ fca37   17          (182, 182)      1  0.0769231
 1093 │ fca37   17          (182, 208)      1  0.0769231
 1094 │ fca37   17          (208, 220)      1  0.0769231
                                        1086 rows omitted
```
"""
function geno_freqtable(data::PopData; by::String = "global")
    if by == "global"
        grp = groupby(dropmissing(data.loci, :genotype), [:locus, :genotype])
        counts = DataFrames.combine(grp, nrow => :count)
        counts = DataFrames.combine(
            groupby(counts, :locus), :genotype, :count,
            :count => (x -> x / sum(x)) => :frequency
        )
    elseif by in ["local", "population"]
        grp = groupby(dropmissing(data.loci, :genotype), [:locus, :population, :genotype])
        counts = DataFrames.combine(grp, nrow => :count)
        counts = DataFrames.combine(
            groupby(counts, [:locus, :population]), :genotype, :count,
            :count => (x -> x / sum(x)) => :frequency
        )
    else
        throw(ArgumentError("Please use by = \"global\" (default) or \"population\""))
    end
end

"""
    allele_freqtable(data::PopData; by::String = "global")
Return a table of the observed `global` (default) or `population` allele frequencies in a PopData object.

### Example:
```
julia> cats = @nancycats ;

julia> allele_freqtable(cats)
108×4 DataFrame
 Row │ locus   allele  count  frequency
     │ String  Int16?  Int64  Float64
─────┼───────────────────────────────────
   1 │ fca8       135    105  0.241935
   2 │ fca8       143     44  0.101382
   3 │ fca8       133     33  0.0760369
   4 │ fca8       137     83  0.191244
  ⋮  │   ⋮       ⋮       ⋮        ⋮
 105 │ fca37      226      2  0.00421941
 106 │ fca37      216      7  0.0147679
 107 │ fca37      224      2  0.00421941
 108 │ fca37      204      6  0.0126582
                         100 rows omitted

julia> allele_freqtable(cats, by = "population")
839×5 DataFrame
 Row │ locus   population  allele  count  frequency
     │ String  String      Int16?  Int64  Float64
─────┼──────────────────────────────────────────────
   1 │ fca8    1              135      9  0.5625
   2 │ fca8    1              143      4  0.25
   3 │ fca8    1              133      2  0.125
   4 │ fca8    1              137      1  0.0625
  ⋮  │   ⋮         ⋮         ⋮       ⋮        ⋮
 836 │ fca37   16             210      5  0.208333
 837 │ fca37   17             208     22  0.846154
 838 │ fca37   17             182      3  0.115385
 839 │ fca37   17             220      1  0.0384615
                                    831 rows omitted
```
"""
function allele_freqtable(data::PopData; by::String = "global")
    if by == "global"
        grp = groupby(dropmissing(data.loci, :genotype), :locus)
        _alleles = DataFrames.combine(grp, :genotype => alleles => :allele, ungroup = false)
        counts = DataFrames.combine(
            _alleles,
            :allele,
            nrow => :n,
            ungroup = false
        )
        counts = DataFrames.combine(
            counts,
            :locus, :allele, :n
        )
        counts = DataFrames.combine(
            groupby(counts, [:locus, :allele]),
            nrow => :count,
            [:n, :allele] => ((n,al) -> length(al)/first(n)) => :frequency
        )
    elseif by in ["local", "population"]
        grp = groupby(dropmissing(data.loci, :genotype), [:locus, :population])
        _alleles = DataFrames.combine(grp, :genotype => alleles => :allele, ungroup = false)
        counts = DataFrames.combine(
            _alleles,
            :allele,
            nrow => :n,
            ungroup = false
        )
        counts = DataFrames.combine(
            counts,
            :locus, :population, :allele, :n
        )
        counts = DataFrames.combine(
            groupby(counts, [:locus, :population, :allele]),
            nrow => :count,
            [:n, :allele] => ((n,al) -> length(al)/first(n)) => :frequency
        )
        #return counts
    else
        throw(ArgumentError("Please use by = \"global\" (default) or \"population\""))
    end
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    abstract type AbstractBridge <: MOI.Bridges.AbstractBridge end

Subtype of [`MathOptInterface.Bridges.AbstractBridge`](@ref) for variable
bridges.
"""
abstract type AbstractBridge <: MOI.Bridges.AbstractBridge end

"""
    supports_constrained_variable(
        BT::Type{<:AbstractBridge},
        S::Type{<:MOI.AbstractSet},
    )::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging
constrained variables in `S`. That is, it returns `true` if the bridge of type
`BT` converts constrained variables of type `S` into a form supported by the
solver.

## Implementation notes

 * This method depends only on the type of the bridge and set, not the runtime
   values.
 * There is a default fallback, so you need only implement this method for sets
   that the bridge implements.

## Example

```jldoctest; setup=(using MathOptInterface; const MOI = MathOptInterface)
julia> MOI.Bridges.Variable.supports_constrained_variable(
           MOI.Bridges.Variable.NonposToNonnegBridge{Float64},
           MOI.Nonpositives,
       )
true

julia> MOI.Bridges.Variable.supports_constrained_variable(
           MOI.Bridges.Variable.NonposToNonnegBridge{Float64},
           MOI.Nonnegatives,
       )
false
```
"""
function supports_constrained_variable(
    ::Type{<:AbstractBridge},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

"""
    concrete_bridge_type(
        BT::Type{<:AbstractBridge},
        S::Type{<:MOI.AbstractSet},
    )::Type

Return the concrete type of the bridge supporting variables in `S` constraints.

This function can only be called if `MOI.supports_constrained_variable(BT, S)`
is `true`.

## Examples

As a variable in [`MathOptInterface.GreaterThan`](@ref) is bridged into
variables in [`MathOptInterface.Nonnegatives`](@ref) by the
[`VectorizeBridge`](@ref):

```jldoctest; setup=:(using MathOptInterface; const MOI = MathOptInterface)
julia> MOI.Bridges.Variable.concrete_bridge_type(
           MOI.Bridges.Variable.VectorizeBridge{Float64},
           MOI.GreaterThan{Float64},
       )
MathOptInterface.Bridges.Variable.VectorizeBridge{Float64, MathOptInterface.Nonnegatives}
```
"""
concrete_bridge_type(::Type{BT}, ::Type{<:MOI.AbstractSet}) where {BT} = BT

function concrete_bridge_type(
    b::MOI.Bridges.AbstractBridgeOptimizer,
    S::Type{<:MOI.AbstractSet},
)
    return concrete_bridge_type(MOI.Bridges.bridge_type(b, S), S)
end

"""
    bridge_constrained_variable(
        BT::Type{<:AbstractBridge},
        model::MOI.ModelLike,
        set::MOI.AbstractSet,
    )::BT

Bridge the constrained variable in `set` using bridge `BT` to `model` and
returns a bridge object of type `BT`.

## Implementation notes

 * The bridge type `BT` must be a concrete type, that is, all the type
   parameters of the bridge must be set.
"""
function bridge_constrained_variable end

function MOI.get(
    ::MOI.ModelLike,
    attr::MOI.AbstractVariableAttribute,
    bridge::AbstractBridge,
)
    return throw(
        ArgumentError(
            "Variable bridge of type `$(typeof(bridge))` does not support " *
            "accessing the attribute `$attr`.",
        ),
    )
end

function MOI.get(
    ::MOI.ModelLike,
    attr::MOI.AbstractVariableAttribute,
    bridge::AbstractBridge,
    ::MOI.Bridges.IndexInVector,
)
    return throw(
        ArgumentError(
            "Variable bridge of type `$(typeof(bridge))` does not support " *
            "accessing the attribute `$attr`.",
        ),
    )
end

function MOI.supports(
    ::MOI.ModelLike,
    ::MOI.AbstractVariableAttribute,
    ::Type{<:AbstractBridge},
)
    return false
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.AbstractVariableAttribute,
    bridge::AbstractBridge,
    value,
    ::MOI.Bridges.IndexInVector...,
)
    if MOI.is_copyable(attr) && !MOI.supports(model, attr, typeof(bridge))
        throw(MOI.UnsupportedAttribute(attr))
    else
        throw(MOI.SetAttributeNotAllowed(attr))
    end
end

"""
    unbridged_map(
       bridge::MOI.Bridges.Variable.AbstractBridge,
        vi::MOI.VariableIndex,
    )

For a bridged variable in a scalar set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the
bridged variable `vi`.

```julia
unbridged_map(
    bridge::MOI.Bridges.Variable.AbstractBridge,
    vis::Vector{MOI.VariableIndex},
)
```

For a bridged variable in a vector set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the bridged
variable `vis`. If this method is not implemented, it falls back to calling
the following method for every variable of `vis`.

```julia
unbridged_map(
    bridge::MOI.Bridges.Variable.AbstractBridge,
    vi::MOI.VariableIndex,
    i::MOI.Bridges.IndexInVector,
)
```

For a bridged variable in a vector set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the bridged
variable `vi` corresponding to the `i`th variable of the vector.

If there is no way to recover the expression in terms of the bridged variable(s)
`vi(s)`, return `nothing`. See [`ZerosBridge`](@ref) for an example of bridge
returning `nothing`.
"""
function unbridged_map end

function unbridged_map(bridge::AbstractBridge, vis::Vector{MOI.VariableIndex})
    mappings = Pair{MOI.VariableIndex,MOI.AbstractScalarFunction}[]
    for (i, vi) in enumerate(vis)
        vi_mappings = unbridged_map(bridge, vi, MOI.Bridges.IndexInVector(i))
        if vi_mappings === nothing
            return
        end
        for mapping in vi_mappings
            push!(mappings, mapping)
        end
    end
    return mappings
end


text:
function record_positions!(
    msd::MeanSquaredDisplacement,
    positions::AbstractArray,
    ttime::AbstractFloat,
)
    # If we record positions, we update the counter here
    msd.naverage += 1
    # We keep record of the time the positions are collected
    msd.wt[msd.naverage, 1] = ttime * msd.interval * (msd.naverage - 1.0)

    # For each average value, we record the positions of all the particles
    # in that specific point in time
    for i in axes(msd.displacement, 2)
        msd.displacement[msd.naverage, i, :] = @view positions[i, :]
    end
end

function difusion!(msd::MeanSquaredDisplacement, sft::SelfScatteringFunction)
    num_particle = size(msd.displacement, 2)

    for i = 1:msd.naverage
        total_msd = 0.0
        total_scattering = 0.0

        for j = 1:msd.naverage-i
            # Obtain the distance between particles along all spatial dimensions
            ddiff = msd.displacement[(j+i), :, :] - msd.displacement[j, :, :]
            # We sum only the dimensions, not the particles
            ddiff = sum(ddiff .^ 2; dims = 2)

            # For the MeanSquareDisplacement, we actually need the sum of all the particles
            total_msd += sum(ddiff)

            # For the SelfScatteringFunction, we need the actual Euclidean distance
            # per particle, so we need to loop and extract the distance for each particle
            @inbounds @fastmath for d in ddiff
                # We already have the sum of square, just need the sqrt to complete
                # the Euclidean distance computation
                distance = √d
                # Given that the system is isotropic, we compute it as such
                total_scattering += sin(distance * sft.κ) / (distance * sft.κ)
            end
        end

        # Normalize the MeanSquaredDisplacement
        norm_val = num_particle * (msd.naverage - i)
        msd.wt[i, 2] = total_msd / norm_val

        # Normalize the SelfScatteringFunction
        sft.dft[i, 2] = total_scattering / norm_val
    end
end

function difusion!(msd::MeanSquaredDisplacement)
    num_particle = size(msd.displacement, 2)

    for i = 1:msd.naverage
        total_msd = 0.0

        @inbounds @fastmath for j = 1:msd.naverage-i
            # Obtain the distance between particles along all spatial dimensions
            ddiff = msd.displacement[(j+i), :, :] - msd.displacement[j, :, :]
            # And we can sum both the dimensions and the particles
            total_msd += sum(ddiff .^ 2)
        end

        # Normalize the MeanSquaredDisplacement
        norm_val = num_particle * (msd.naverage - i)
        msd.wt[i, 2] = total_msd / norm_val
    end
end


text:
module FilesWrapper

using JavaCall

include("CcsdsWrapper/CcsdsWrapper.jl")
include("GeneralWrapper/GeneralWrapper.jl")
include("Sp3Wrapper/Sp3Wrapper.jl")

end


text:
## Setup the Julia optimization software
#=
 The following packages are for Julia optimization software
 http://www.juliaopt.org/
=#

Pkg.add("Optim")     # This is the closest to the software we'll build in this class

Pkg.add("JuMP")         # high level modeling language like CVX, AMPL
Pkg.add("Convex")       # Convex models like CVX
Pkg.add("MathProgBase") # probably included in above, but why not?

# These all involve installing software and are likely to be less reliable.
# But don't let that discourage you! Give them try as the authors are all
# committeed to wide compatiblity and will likely be amenable to help solving
# problems.

Pkg.add("NLopt")
Pkg.add("Ipopt")
Pkg.add("Clp")
Pkg.add("SCS")
Pkg.add("GLPK")
Pkg.add("CSDP")


text:
#!/usr/bin/env julia

# Create image data VTK file.
# This corresponds to a rectilinear grid with uniform spacing in each direction.

using WriteVTK
import Compat.UTF8String

function main()
    Ni, Nj, Nk = 20, 30, 40
    outfiles = UTF8String[]

    # Initialise random number generator with deterministic seed (for
    # reproducibility).
    rng = MersenneTwister(42)
    data = randn(rng, Ni, Nj, Nk)

    # Write 2D and 3D arrays.
    append!(outfiles, vtk_write_array("arraydata_2D", data[:, :, 1], "array2D"))
    append!(outfiles, vtk_write_array("arraydata_3D", data, "array3D"))

    println("Saved:  ", join(outfiles, "  "))

    return outfiles
end

ProjDir = dirname(@__FILE__)
cd(ProjDir*"/output") do

  main()

end


text:
function autotransition_init(K::Integer, order::Integer, partition::String="")
        a = rand(Dirichlet(ones(K)/K)) #uninformative prior on initial state probabilities
        A = strong_autotrans_matrix(K)
        no_emission_symbols = Int(4^(order+1)) #alphabet size for the order
        emission_dists = [generate_emission_dist(no_emission_symbols) for i in 1:K]
        #generate the BHMM with the appropriate transition matrix and emissions distributions
        hmm = BHMM(a, A, emission_dists, partition)
end
            #function to construct BHMM transition matrix with strong priors on auto-transition
            function strong_autotrans_matrix(states::Integer, prior_dope::AbstractFloat=(states*250.0), prior_background::AbstractFloat=.1)
                transition_matrix=zeros(states,states)
                for k in 1:states
                    dirichlet_params = fill(prior_background, states)
                    dirichlet_params[k] = prior_dope
                    transition_matrix[k,:] = rand(Dirichlet(dirichlet_params))
                end
                return transition_matrix
            end

            #function to construct BHMM state emission distribution from uninformative dirichlet over the alphabet size
            function generate_emission_dist(no_emission_symbols, prior=Dirichlet(ones(no_emission_symbols)/no_emission_symbols))
                return Categorical(rand(prior))
            end


text:
# Accelerated component-wise product

C = zeros(m,p)

for j=1:m
    for l=1:p
        for k=1:n
            @inbounds C[j,l] += A[j,k]*B[k,l]
        end
    end
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Numerical Linear Algebra/ComponentProduct.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
#####
##### Prognostic Variables
#####

prognostic_vars(water::PrescribedWaterModel) = ()
prognostic_vars(water::SoilWaterModel) =
    (VolumetricLiquidFraction(), VolumetricIceFraction())
prognostic_vars(heat::PrescribedTemperatureModel) = ()
prognostic_vars(heat::SoilHeatModel) = (VolumetricInternalEnergy(),)

prognostic_vars(land::LandModel) =
    (prognostic_vars(land.soil.water)..., prognostic_vars(land.soil.heat)...)

get_prog_state(state, ::VolumetricLiquidFraction) = (state.soil.water, :ϑ_l)
get_prog_state(state, ::VolumetricIceFraction) = (state.soil.water, :θ_i)
get_prog_state(state, ::VolumetricInternalEnergy) = (state.soil.heat, :ρe_int)


text:
using ITensors

import ITensors: optimal_contraction_sequence

# This is for daggering and Empty ITensor
using ITensors.NDTensors
Base.conj(s::Empty; kwargs...) = s

# TODO: use ITensors.contraction_cost
contraction_cost(A::Vector{<: ITensor}, sequence) = _contraction_cost(A, sequence, 0)
function _contraction_cost(A::Vector{<: ITensor}, sequence, previouscost)
  inds1, cost1 = _contraction_cost(A, sequence[1], previouscost)
  inds2, cost2 = _contraction_cost(A, sequence[2], previouscost)
  return _contraction_cost(inds1, inds2, previouscost + cost1 + cost2)
end
_contraction_cost(As::Vector{<: ITensor}, sequence::Int, previouscost) = inds(As[sequence]), previouscost
function _contraction_cost(A1::IndexSet, A2::IndexSet, previouscost)
  return noncommoninds(A1, A2), previouscost + dim(unioninds(A1, A2))
end

# "Contract" a set of labels according to a sequence
function contract_labels(network_labels::Vector{Symbol}, sequence)
  branch1 = contract_labels(network_labels, sequence[1])
  branch2 = contract_labels(network_labels, sequence[2])
  contract_labels(branch1, branch2)
end
contract_labels(network_labels::Vector{Symbol}, n::Int) = network_labels[n]
contract_labels(label1::Symbol, label2::Symbol) = (label1, label2)
contract_labels(label1::Symbol, label2::Tuple) = (label1, label2)
contract_labels(label1::Tuple, label2::Tuple) = (label1, label2)

# Examples from Section III of https://arxiv.org/abs/1304.6112

english_alphabet = Char.(97:97+25)
# Include Π (Char(928))
# Include ϕ (Char(981))
# Include ϵ (Char(1013))
# Don't include \delta Char(948) since it conflicts with ITensors.δ
# Don't include \pi Char(960) since it conflicts with MathConstants.π
greek_alphabet = push!(vcat(Char.(945:947), Char.(949:959), Char.(961:969)), 928, 981, 1013)
greek_alphabet_uppercase = Char.(913:913+24)
indlabels = vcat(english_alphabet, greek_alphabet, greek_alphabet_uppercase)

# Add labels with overbars as well
# \overbar Char(773)
indlabels = vcat(indlabels, indlabels .* Char(773))

function define_inds(dχ, indlabels)
  for i in indlabels
    @eval begin
      $(Symbol(i)) = Index($dχ, string($i))
    end
  end
end

# Fig. 6(a)
# Optimal sequence: Any[6, Any[Any[1, 2], Any[5, Any[3, 4]]]]
# Optimal cost: 10χ^3+ 16χ^2
function tebd(dχ, ds = 2; indlabels = indlabels)
  define_inds(dχ, indlabels)

  cₛ = Index(ds, "c")
  fₛ = Index(ds, "f")
  iₛ = Index(ds, "i")
  jₛ = Index(ds, "j")

  λ¹ = randomITensor(a, b)
  Γ¹ = randomITensor(b, d, cₛ)
  λ² = randomITensor(d, e)
  Γ² = randomITensor(e, g, fₛ)
  λ³ = randomITensor(g, h)
  U = randomITensor(cₛ, fₛ, iₛ, jₛ)
  network = [λ¹, Γ¹, λ², Γ², λ³, U]
  sequence = @time optimal_contraction_sequence(network)

  # Analyze the sequence
  network_labels = [:λ¹, :Γ¹, :λ², :Γ², :λ³, :U]
  sequence_labels = @show contract_labels(network_labels, sequence)
  sequence_inds, sequence_cost = contraction_cost(network, sequence)
  @show hassameinds(sequence_inds, noncommoninds(network...))
  @show sequence_cost
  @show 10dχ^3+ 16dχ^2
  return (sequence = sequence, sequence_labels = sequence_labels, sequence_cost = sequence_cost)
end

# Fig. 6(b)
# Optimal sequence:  Any[3, Any[5, Any[2, Any[1, 4]]]]
# Optimal cost: 4χ^6
function ttn_1d_3_to_1(dχ; indlabels = indlabels)
  define_inds(dχ, indlabels)

  W² = emptyITensor(c, d, e, f)
  H = emptyITensor(a, b, k, c)
  W¹ᴴ = emptyITensor(i, j, a, h)
  W²ᴴ = replaceinds(dag(W²), (c, f) => (b, g))
  Ρ = emptyITensor(l, f, h, g)
  network = [W², H, W¹ᴴ, W²ᴴ, Ρ]
  sequence = @time optimal_contraction_sequence(network)

  # Analyze the sequence
  network_labels = [:W², :H, :W¹ᴴ, :W²ᴴ, :Ρ]
  sequence_labels = @show contract_labels(network_labels, sequence)
  sequence_inds, sequence_cost = contraction_cost(network, sequence)
  @show hassameinds(sequence_inds, noncommoninds(network...))
  @show sequence_cost
  @show 4dχ^6
  return (sequence = sequence, sequence_labels = sequence_labels, sequence_cost = sequence_cost)
end

# Fig. 7(a)
# Optimal sequence:  Any[5, Any[Any[4, Any[1, 6]], Any[Any[2, 7], Any[9, Any[3, 8]]]]]
# Optimal cost: 4χ^12+ 4χ^10
# WARNING: this takes a while to run
function ttn_2d_9_to_1(dχ; indlabels = indlabels)
  define_inds(dχ, indlabels)
  W² = emptyITensor(ϵ, ζ, η, θ, ι, κ, j, λ, μ, c)
  W³ = emptyITensor(o, p, h, q, r, s, t, u, v, a)
  W⁴ = emptyITensor(i, w, x, y, z, α, β, γ, Δ, b)
  H = emptyITensor(m, n, k, l, χ, j, h, i)
  W¹ᴴ = dag(emptyITensor(ν, ξ, Π, ρ, σ, τ, υ, ϕ, m, f))
  W²ᴴ = dag(replaceinds(W², (j, c) => (n, g)))
  W³ᴴ = dag(replaceinds(W³, (h, a) => (k, d)))
  W⁴ᴴ = dag(replaceinds(W⁴, (i, b) => (l, e)))
  Ρ = emptyITensor(ψ, c, a, b, f, g, d, e)
  network = [W², W³, W⁴, H, W¹ᴴ, W²ᴴ, W³ᴴ, W⁴ᴴ, Ρ]
  sequence = @time optimal_contraction_sequence(network)

  # Analyze the sequence
  network_labels = [:W², :W³, :W⁴, :H, :W¹ᴴ, :W²ᴴ, :W³ᴴ, :W⁴ᴴ, :Ρ]
  sequence_labels = @show contract_labels(network_labels, sequence)
  sequence_inds, sequence_cost = contraction_cost(network, sequence)
  @show hassameinds(sequence_inds, noncommoninds(network...))
  @show sequence_cost
  @show 4dχ^12+ 4dχ^10
  return (sequence = sequence, sequence_labels = sequence_labels, sequence_cost = sequence_cost)
end

# Fig. 7(b)
#
# TODO
#W(2) αδεζηθ
#W(3)  βικλμν
#W(4)  γξπρστ
#V(1)   ̄πιωabc
#V(2)   ̄ξευφχψ
#V(3)  θξdefg
#V(4)  μπhijk
#U     ψbehlmno
#h      ̄νφωlpqrs
#U∗    smnotuvw
#V(1)∗ raucxy
#V(2)∗ υqχtz ̄α
#V(3)∗ dvfg ̄β ̄γ
#V(4)∗ wijk ̄δ ̄ε
#W(1)∗  ̄λ ̄μpzx ̄ζ
#W(2)∗ δ ̄αζη ̄β ̄η
#W(3)∗ yκλ ̄δν ̄θ
#W(4)∗  ̄γ ̄ερστ ̄ι
#ρ      ̄ζ ̄η ̄θ ̄ι ̄καβγ

# Fig. 7(c)
# 4χ^26 + 2χ^25 + 2χ^23+ 3χ^22 + 3χ^20 + χ^16 + χ^14 + χ^13 + χ^12 + 4χ^8 + 4χ^7
function mera_2d_4_to_1(dχ; indlabels = indlabels)
  define_inds(dχ, indlabels)
  W²  =     emptyITensor(σ̅, d̅, e̅, ξ, σ)
  W³  =     emptyITensor(τ̅, f̅, g̅, τ, j̅)
  W⁴  =     emptyITensor(υ̅, k̅, Π, m̅, ω)
  W⁵  =     emptyITensor(φ̅, ρ, υ, a, f)
  W⁶  =     emptyITensor(χ̅, φ, l̅, g, n̅)
  W⁷  =     emptyITensor(ψ̅, o̅, b, q̅, r̅)
  W⁸  =     emptyITensor(ω̅, c, h, s̅, t̅)
  W⁹  =     emptyITensor(a̅, i, p̅, u̅, v̅)
  U¹  =     emptyITensor(ν, ξ, Π, ρ, κ̅, α, γ, λ̅)
  U²  =     emptyITensor(σ, τ, υ, φ, β, χ, μ̅, ψ)
  U³  =     emptyITensor(ω, a, b, c, ν̅, ξ̅, d, e)
  U⁴  =     emptyITensor(f, g, h, i, Π̅, j, k, l)
  H   =     emptyITensor(κ̅, α, β, γ, λ̅, μ̅, ν̅, ξ̅, Π̅, Δ, ε, ζ, η, θ, ι, κ, λ, μ)
  U¹ᴴ = dag(emptyITensor(Δ, ε, η, θ, m, n, o, p))
  U²ᴴ = dag(emptyITensor(ζ, χ, ι, ψ, s, t, u, v))
  U³ᴴ = dag(emptyITensor(κ, λ, d, e, y, z, α̅, β̅))
  U⁴ᴴ = dag(emptyITensor(μ, j, k, l, ζ̅, η̅, θ̅, ι̅))
  W¹ᴴ = dag(emptyITensor(b̅, c̅, h̅, m, w̅))
  W²ᴴ = dag(emptyITensor(d̅, e̅, n, s, x̅))
  W³ᴴ = dag(emptyITensor(f̅, g̅, t, j̅, q))
  W⁴ᴴ = dag(emptyITensor(k̅, o, m̅, y, r))
  W⁵ᴴ = dag(emptyITensor(p, u, z, ζ̅, w))
  W⁶ᴴ = dag(emptyITensor(v, l̅, η̅, n̅, x))
  W⁷ᴴ = dag(emptyITensor(o̅, α̅, q̅, r̅, γ̅))
  W⁸ᴴ = dag(emptyITensor(β̅, θ̅, s̅, t̅, Δ̅))
  W⁹ᴴ = dag(emptyITensor(ι̅, p̅, u̅, v̅, ε̅))
  Ρ   =     emptyITensor(w̅, x̅, q, r, w, x, γ̅, Δ̅, ε̅, ρ̅, σ̅, τ̅, υ̅, φ̅, χ̅, ψ̅, ω̅, a̅)
  network = [W², W³, W⁴, W⁵, W⁶, W⁷, W⁸, W⁹, U¹, U², U³, U⁴, H, U¹ᴴ, U²ᴴ, U³ᴴ, U⁴ᴴ, W¹ᴴ, W²ᴴ, W³ᴴ, W⁴ᴴ, W⁵ᴴ, W⁶ᴴ, W⁷ᴴ, W⁸ᴴ, W⁹ᴴ, Ρ]
  sequence = @time optimal_contraction_sequence(network)

  # Analyze the sequence
  network_labels = [:W², :W³, :W⁴, :W⁵, :W⁶, :W⁷, :W⁸, :W⁹, :U¹, :U², :U³, :U⁴, :H, :U¹ᴴ, :U²ᴴ, :U³ᴴ, :U⁴ᴴ, :W¹ᴴ, :W²ᴴ, :W³ᴴ, :W⁴ᴴ, :W⁵ᴴ, :W⁶ᴴ, :W⁷ᴴ, :W⁸ᴴ, :W⁹ᴴ, :Ρ]
  sequence_labels = @show contract_labels(network_labels, sequence)
  sequence_inds, sequence_cost = contraction_cost(network, sequence)
  @show hassameinds(sequence_inds, noncommoninds(network...))
  @show sequence_cost
  @show 4dχ^26 + 2dχ^25 + 2dχ^23 + 3dχ^22 + 3dχ^20 + dχ^16 + dχ^14 + dχ^13 + dχ^12 + 4dχ^8+ 4dχ^7
  return (sequence = sequence, sequence_labels = sequence_labels, sequence_cost = sequence_cost)
end


text:
"""
ClusterManager for a Slurm allocation

Represents the resources available within a slurm allocation created by salloc/sbatch.\\
The environment variables `SLURM_JOB_ID` or `SLURM_JOBID` and `SLURM_NTASKS` must be defined to construct this object.
"""
mutable struct SlurmManager <: ClusterManager
  jobid::Int
  ntasks::Int
  verbose::Bool
  launch_timeout::Float64
  srun_post_exit_sleep::Float64
  srun_proc

  function SlurmManager(;verbose=false, launch_timeout=60.0, srun_post_exit_sleep=0.01)

    jobid =
    if "SLURM_JOB_ID" in keys(ENV)
        ENV["SLURM_JOB_ID"]
    elseif "SLURM_JOBID" in keys(ENV)
        ENV["SLURM_JOBID"]
    else
        error("""
              SlurmManager must be constructed inside a slurm allocation environemnt.
              SLURM_JOB_ID or SLURM_JOBID must be defined.
              """)
    end

    ntasks =
    if "SLURM_NTASKS" in keys(ENV)
      ENV["SLURM_NTASKS"]
    else
      error("""
            SlurmManager must be constructed inside a slurm environment with a specified number of tasks.
            SLURM_NTASKS must be defined.
            """)
    end

    jobid = parse(Int, jobid)
    ntasks = parse(Int, ntasks)

    new(jobid, ntasks, verbose, launch_timeout, srun_post_exit_sleep, nothing)
  end
end

function launch(manager::SlurmManager, params::Dict, instances_arr::Array, c::Condition)
    try
        exehome = params[:dir]
        exename = params[:exename]
        exeflags = params[:exeflags]

        srun_cmd = `srun -D $exehome $exename $exeflags --worker=$(cluster_cookie())`
        manager.srun_proc = open(srun_cmd)

        t = @async for i in 1:manager.ntasks
          manager.verbose && println("connecting to worker $i out of $(manager.ntasks)")

          line = readline(manager.srun_proc)
          m = match(r".*:(\d*)#(.*)", line)
          m === nothing && error("could not parse $line")

          config = WorkerConfig()
          config.port = parse(Int, m[1])
          config.host = strip(m[2])

          push!(instances_arr, config)
          notify(c)
        end

        # workers must be launched before timeout otherwise interrupt
        status = timedwait(() -> istaskdone(t), manager.launch_timeout)
        if status !== :ok
          @async Base.throwto(t, ErrorException("launch_timeout exceeded"))
        end
        wait(t)

        # redirect output
        @async while !eof(manager.srun_proc)
          line = readline(manager.srun_proc)
          println(line)
        end

        # wait to make sure that srun_proc exits before main program to avoid slurm complaining
        # avoids "Job step aborted: Waiting up to 32 seconds for job step to finish" message
        finalizer(manager) do manager
          wait(manager.srun_proc)
          # need to sleep briefly here to make sure that srun exit is recorded by slurm daemons
          # TODO find a way to wait on the condition directly instead of just sleeping
          sleep(manager.srun_post_exit_sleep)
        end

    catch e
        println("Error launching Slurm job:")
        rethrow(e)
    end
end

function manage(manager::SlurmManager, id::Integer, config::WorkerConfig, op::Symbol)
    # This function needs to exist, but so far we don't do anything
end


text:
abstract StreamStat{
    F <: Distributions.VariateForm,
    S <: Distributions.ValueSupport
}

typealias DiscreteUnivariateStreamStat StreamStat{
    Distributions.Univariate,
    Distributions.Discrete
}

typealias ContinuousUnivariateStreamStat StreamStat{
    Distributions.Univariate,
    Distributions.Continuous
}

typealias DiscreteMultivariateStreamStat StreamStat{
    Distributions.Multivariate,
    Distributions.Discrete
}

typealias ContinuousMultivariateStreamStat StreamStat{
    Distributions.Multivariate,
    Distributions.Continuous
}


text:
using Distributions

function surrogate_splits(y_obs_split::Vector, X::DataFrame, col_indcs::Vector{Int}, max_surrogates::Int, weights::Vector)
    p = ncol(X)
    n_surr = p - 1 < max_surrogates ? p - 1: max_surrogates
    surr = Array{Tuple}(n_surr)

    for i = 1:n_surr
        surr[i] = _split_classifcation_error_loss(y_obs_split, X, col_indcs, weights)
        col_indcs = setdiff(col_indcs, surr[i][1])          # considerable mem-alloc here
    end
    return surr
end

# This method is dispatched when weights are omitted. This
# allows us to compute the loss function 5x faster
function surrogate_splits(y_obs_split::Vector, X::DataFrame, row_indcs::Vector{Int}, col_indcs::Vector{Int}, max_surrogates::Int)
    p = ncol(X)
    n_surr = p - 1 < max_surrogates ? p - 1: max_surrogates
    surr = Array{Tuple}(n_surr)

    for i = 1:n_surr
        surr[i] = _split_classifcation_error_loss(y_obs_split, X, row_indcs, col_indcs)
        col_indcs = setdiff(col_indcs, surr[i][1])
    end
    return surr
end

# For now, this function doesn't solve the edge case in which
# a given record has missing data on all surrogates.
function apply_surrogates(split_with_na::Vector, X::DataFrame, surr::Array{Tuple, 1})
    n = length(split_with_na)
    col_indcs = [x[1] for x in surr]
    col_thresh = [x[2] for x in surr]

    split = falses(n)

    for i = 1:n
        if isna(split_with_na[i])
            for (idx, j) in enumerate(col_indcs)

                if j != 0
                    if !isna(X[i, j])
                        split[i] = X[i, j] .< col_thresh[idx]
                        break
                    end
                # when NO_BEST was result from surrogate_splits() we
                # assign the split value at random
                elseif j == 0
                    split[i] = bitrand(1)[1]
                end
            end
        else
            split[i] = split_with_na[i]
        end
    end
    return split
end

count_surrogates(node::Node) = isdefined(node.surrogates) ? length(node.surrogates) : 0


text:
function LogNormalDistribution()
    return LogNormalDistribution(())
end

function LogNormalDistribution(arg0::jdouble, arg1::jdouble)
    return LogNormalDistribution((jdouble, jdouble), arg0, arg1)
end

function LogNormalDistribution(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return LogNormalDistribution((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function cumulative_probability(obj::LogNormalDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::LogNormalDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_location(obj::LogNormalDistribution)
    return jcall(obj, "getLocation", jdouble, ())
end

function get_numerical_mean(obj::LogNormalDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::LogNormalDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_shape(obj::LogNormalDistribution)
    return jcall(obj, "getShape", jdouble, ())
end

function get_support_lower_bound(obj::LogNormalDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::LogNormalDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function is_support_connected(obj::LogNormalDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::LogNormalDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end

function probability(obj::LogNormalDistribution, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "probability", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
function rotate_result(
    model,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintDual},
    ci::MOI.ConstraintIndex,
)
    x = MOI.get(model, attr, ci)
    T = eltype(x)
    s2 = √T(2)
    return [x[1] / s2 + x[2] / s2; x[1] / s2 - x[2] / s2; x[3:end]]
end

function rotate_result(
    model,
    attr::MOI.VariablePrimal,
    variables,
    i::IndexInVector,
)
    if i.value == 1 || i.value == 2
        t, u = MOI.get(model, attr, variables[1:2])
        T = typeof(t)
        s2 = √T(2)
        if i.value == 1
            return t / s2 + u / s2
        else
            return t / s2 - u / s2
        end
    else
        return MOI.get(model, attr, variables[i.value])
    end
end

function rotate_bridged_function(T::Type, variables, i::IndexInVector)
    s2 = √T(2)
    if i.value == 1 || i.value == 2
        t = MOIU.operate(/, T, MOI.SingleVariable(variables[1]), s2)
        u = MOIU.operate(/, T, MOI.SingleVariable(variables[2]), s2)
        if i.value == 1
            return MOIU.operate!(+, T, t, u)
        else
            return MOIU.operate!(-, T, t, u)
        end
    else
        return convert(
            MOI.ScalarAffineFunction{T},
            MOI.SingleVariable(variables[i.value]),
        )
    end
end

function rotate_unbridged_map(
    T::Type,
    inner_variables::Vector{MOI.VariableIndex},
    bridged_variables::Vector{MOI.VariableIndex},
)
    F = MOI.ScalarAffineFunction{T}
    umap = Pair{MOI.VariableIndex,F}[]
    s2 = √T(2)
    t = MOIU.operate(/, T, MOI.SingleVariable(bridged_variables[1]), s2)
    u = MOIU.operate(/, T, MOI.SingleVariable(bridged_variables[2]), s2)
    push!(umap, inner_variables[1] => MOIU.operate(+, T, t, u))
    push!(umap, inner_variables[2] => MOIU.operate(-, T, t, u))
    for i in 3:length(bridged_variables)
        func = convert(F, MOI.SingleVariable(bridged_variables[i]))
        push!(umap, inner_variables[i] => func)
    end
    return umap
end

"""
    RSOCtoSOCBridge{T} <: Bridges.Variable.AbstractBridge

Same transformation as [`MOI.Bridges.Constraint.RSOCBridge`](@ref).
"""
struct RSOCtoSOCBridge{T} <: AbstractBridge
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.SecondOrderCone}
end
function bridge_constrained_variable(
    ::Type{RSOCtoSOCBridge{T}},
    model::MOI.ModelLike,
    set::MOI.RotatedSecondOrderCone,
) where {T}
    variables, constraint = MOI.add_constrained_variables(
        model,
        MOI.SecondOrderCone(MOI.dimension(set)),
    )
    return RSOCtoSOCBridge{T}(variables, constraint)
end

function supports_constrained_variable(
    ::Type{<:RSOCtoSOCBridge},
    ::Type{MOI.RotatedSecondOrderCone},
)
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:RSOCtoSOCBridge})
    return [(MOI.SecondOrderCone,)]
end
function MOIB.added_constraint_types(::Type{<:RSOCtoSOCBridge})
    return Tuple{DataType,DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::RSOCtoSOCBridge, ::MOI.NumberOfVariables)
    return length(bridge.variables)
end
function MOI.get(bridge::RSOCtoSOCBridge, ::MOI.ListOfVariableIndices)
    return bridge.variables
end
function MOI.get(
    bridge::RSOCtoSOCBridge,
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone},
)
    return 1
end
function MOI.get(
    bridge::RSOCtoSOCBridge,
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.SecondOrderCone},
)
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::RSOCtoSOCBridge)
    return MOI.delete(model, bridge.variables)
end

# Attributes, Bridge acting as a constraint

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::RSOCtoSOCBridge{T},
) where {T}
    return MOI.RotatedSecondOrderCone(length(bridge.variables))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintDual},
    bridge::RSOCtoSOCBridge,
)
    return rotate_result(model, attr, bridge.constraint)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimal,
    bridge::RSOCtoSOCBridge,
    i::IndexInVector,
)
    return rotate_result(model, attr, bridge.variables, i)
end

function MOIB.bridged_function(
    bridge::RSOCtoSOCBridge{T},
    i::IndexInVector,
) where {T}
    return rotate_bridged_function(T, bridge.variables, i)
end
function unbridged_map(
    bridge::RSOCtoSOCBridge{T},
    vis::Vector{MOI.VariableIndex},
) where {T}
    return rotate_unbridged_map(T, bridge.variables, vis)
end


text:
#=

Macro interface is inspired by InfiniteOpt.jl. Many of the code structures are
reused. Transitioning this over to a fully InfiniteOpt.jl in the future is desirable
but most of our use cases involve nonlinear functions so that InfiniteOpt.jl
nonlinear support is gating.

The general idea here is pretty simple. When we add a variable to the model
the user specifies whether or not it is a decision variable or an uncertain
variable we then classify constraints. See https://github.com/PSORLab/EAGO.jl/issues/61
for a discussion of the roadmap here.
=#

#=
@enum(SIPCons, DECISION, UNCERTAIN, SEMIINFINITE, SIPNOTSET)
const MOIU_AUTO_CACHE = MOIU.CachingOptimizerMode.AUTOMATIC

Base.@kwdef mutable struct SIPData
    p::Dict{JuMP.VariableRef,Bool}                  = Dict{JuMP.VariableRef,Bool}()
    constr_type::Dict{JuMP.ConstraintRef, SIPCons}  = Dict{JuMP.ConstraintRef, SIPCons}()
    nl_constraint_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
    nl_expression_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
    model_decision::Model                           = Model(caching_mode = MOIU_AUTO_CACHE)
    model_uncertain::Model                          = Model((aching_mode = MOIU_AUTO_CACHE)
    model_semiinfinite::Model                       = Model(caching_mode = MOIU_AUTO_CACHE)
end
uncertain_variable_num(d::SIPData) = any(x -> x, d.p)

function _copy_nlp_objective!(out::JuMP.Model, m::JuMP.Model, nlp_data::JuMP._NLPData)
    if nlp_data.nlobj !== nothing
        # TODO CHECK THAT THE FORMULATION IS PEACHY
        obj_expr = MOI.objective_expr(nlp_new.evaluator)
        # obj_expr_jref = substitute moi variables to jump refs...
        @NLobjective(m.ext[:sip].model_decision, Min, obj_expr_jref)
    end
    return nothing
end

function _copy_nlp_constraints!(m::JuMP.Model, nlp_data::JuMP._NLPData)
    for i = 1:length(nlp_data.nlconstr)
        con_expr = MOI.constraint_expr(nlp_data.evaluator, i)
    end
    return nothing
end

"""

Populates variables and constraints into one of three models based on classification:
1) store decision variables and constraints containing only decision variables in
`m.model_decision`, 2) store uncertain variables and constraints containing only
uncertain variables in `m.model_uncertain`, 3) store all variables and all
constraints containing both decision variables to `m.model_semiinfinite`.
"""
function initialize_pure_models!(m::JuMP.Model)
    if haskey(m.ext, :sip)

        # extract sip data and storage models
        sip_data = _get_sip_data(m)::SIPData
        m_dec = sip_data.model_decision
        m_unc = sip_data.model_uncertain
        m_sip = sip_data.model_semiinfinite

        # get index map
        indx_map_dec = MOI.copy_to(backend(m_dec), backend(m), copy_names = true)
        indx_map_unc = MOI.copy_to(backend(m_unc), backend(m), copy_names = true)
        indx_map_sip = MOI.copy_to(backend(m_sip), backend(m), copy_names = true)

        # copy extension data
        for (key, data) in m.ext
            if key != :sip
                m_dec.ext[key] = copy_extension_data(data, m_dec, m)
                m_unc.ext[key] = copy_extension_data(data, m_unc, m)
                m_sip.ext[key] = copy_extension_data(data, m_sip, m)
            end
        end

        # copy objects (TODO need to check on SIP, DEC, UNC)
        reference_map = ReferenceMap(new_model, index_map)
        for (name, value) in object_dictionary(model)
            new_model[name] = getindex.(reference_map, value)
        end

        if m.nlp_data !== nothing
            nlp_data = deepcopy(m.nlp_data)
            nlp_data.evaluator = NLPEvaluator(m)
            MOI.initialize(nlp_data.evaluator, Symbol[:ExprGraph])
            copy_nlp_objective!(m, nlp_data)
            copy_nlp_expressions!(m, nlp_data)
            copy_nlp_constraints!(m, nlp_data)
        end

        for (vi, is_uncertain) in m.ext.p
            model_semiinfinite
            model_decision_
    end
    return nothing
end

function initialize_sip_data(m::JuMP.Model)
    m.ext[:sip] = SIPData()
end

function sip_optimizehook(m::JuMP.Model; kwargs...)
    data = _getsipdata(m)::SIPData
    if uncertain_variable_num(data) == 0.0
        ret = JuMP.optimize!(m::JuMP.Model, ignore_optimize_hook = true, kwargs...)
    else
        initialize_pure_models!(m)
        ret = sip_solve(data)
    end
    return ret
end

function enable_semiinfinite(m::JuMP.Model)
    haskey(m.ext, :sip) && error("Model already has semiinfinite programs enabled")
    initialize_sip_data(m)
    JuMP.set_optimize_hook(m, sip_optimizehook)
    return nothing
end

function ModelWithSIP(args...; kwargs...)
    m = JuMP.Model(args...; kwargs...)
    set_optimizer(SIPOptimizer)
    enable_semiinfinite(m)
    return m
end

"""
    @uncertain_variable

Add an *anonymous* variable to `m::SIPModel` described by the keyword
arguments `kw_args` and returns the variable.

    @uncertain_variable(m, expr, args..., kw_args...)

This macro simply denotes the variable as belonging to the set of uncertain
variables and then performs `@variable(m, expr, args..., kw_args...)` to
add the variable to the basic JuMP model for storage. All JuMP syntax is
supported.
"""
macro uncertain_variable(args...)
    esc(quote
        vi = @variable($(args...))
        $(args[1]).ext[:sip].p[vi] = true
        vi
        end)
end

"""
    @decision_variable

Add an *anonymous* variable to `m::SIPModel` described by the keyword
arguments `kw_args` and returns the variable.

    @decision_variable(m, expr, args..., kw_args...)

This macro simply denotes the variable as belonging to the set of decision
variables and then performs `@variable(m, expr, args..., kw_args...)` to
add the variable to the basic JuMP model for storage. All JuMP syntax is
supported.
"""
macro decision_variable(args...)
    esc(quote
        vi = @variable($(args...))
        $(args[1]).ext[:sip].p[vi] = false
        vi
        end)
end
=#