text:
# Quadratic decision function #
# ========================== #

# Helper structs to differentiate between
# QuadraticDecisionFunctions that are ultimately
# linear or
abstract type LinearQuadraticPart{T} end
struct LinearPart{T} <: LinearQuadraticPart{T}
    variable_part::MOI.ScalarAffineFunction{T}
    decision_part::MOI.ScalarAffineFunction{T}
end
struct QuadraticPart{T} <: LinearQuadraticPart{T}
    variable_part::MOI.ScalarQuadraticFunction{T}
    decision_part::MOI.ScalarQuadraticFunction{T}
    cross_terms::MOI.ScalarQuadraticFunction{T}
end

struct QuadraticDecisionFunction{T, LQ <: LinearQuadraticPart{T}} <: MOI.AbstractScalarFunction
    linear_quadratic_terms::LQ
    known_part::MOI.ScalarQuadraticFunction{T}
    known_variable_terms::MOI.ScalarQuadraticFunction{T}
    known_decision_terms::MOI.ScalarQuadraticFunction{T}
end

# Base overrides #
# ========================== #
function Base.copy(lq::L) where L <: LinearPart
    return L(copy(lq.variable_part),
             copy(lq.decision_part))
end
function Base.copy(lq::Q) where Q <: QuadraticPart
    return Q(copy(lq.variable_part),
             copy(lq.decision_part),
             copy(lq.cross_terms))
end
function Base.copy(f::F) where F <: QuadraticDecisionFunction
    return F(copy(f.linear_quadratic_terms),
             copy(f.known_part),
             copy(f.known_variable_terms),
             copy(f.known_decision_terms))
end

function Base.convert(::Type{LinearPart{T}}, α::T) where T
    return LinearPart{T}(
        convert(MOI.ScalarAffineFunction{T}, α),
        convert(MOI.ScalarAffineFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticPart{T}}, α::T) where T
    return QuadraticPart{T}(
        convert(MOI.ScalarQuadraticFunction{T}, α),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}}, α::T) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T,LQ}(
        convert(LQ, α),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{LinearPart{T}},
                      f::MOI.SingleVariable) where T
    return LinearPart{T}(
        convert(MOI.ScalarAffineFunction{T}, f),
        convert(MOI.ScalarAffineFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticPart{T}},
                      f::MOI.SingleVariable) where T
    return QuadraticPart{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}},
                      f::MOI.SingleVariable) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T}(
        convert(LQ, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{LinearPart{T}},
                      f::SingleDecision) where T
    return LinearPart{T}(
        convert(MOI.ScalarAffineFunction{T}, zero(T)),
        convert(MOI.ScalarAffineFunction{T}, f))
end
function Base.convert(::Type{QuadraticPart{T}},
                      f::SingleDecision) where T
    return QuadraticPart{T}(
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}},
                      f::SingleDecision) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T}(
        convert(LQ, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}},
                      f::SingleKnown) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T}(
        convert(LQ, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{LinearPart{T}},
                      f::MOI.ScalarAffineFunction{T}) where T
    return LinearPart{T}(
        convert(MOI.ScalarAffineFunction{T}, f),
        convert(MOI.ScalarAffineFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticPart{T}},
                      f::MOI.ScalarAffineFunction{T}) where T
    return QuadraticPart{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}},
                      f::MOI.ScalarAffineFunction{T}) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T}(
        convert(LQ, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{QuadraticPart{T}},
                      f::MOI.ScalarQuadraticFunction{T}) where T
    return QuadraticPart{T}(
        f,
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,QuadraticPart{T}}},
                      f::MOI.ScalarQuadraticFunction{T}) where T
    return QuadraticDecisionFunction{T}(
        convert(QuadraticPart{T}, f),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{LinearPart{T}},
                      f::AffineDecisionFunction{T}) where T
    return LinearPart{T}(
        f.variable_part,
        f.decision_part)
end
function Base.convert(::Type{QuadraticPart{T}},
                      f::AffineDecisionFunction{T}) where T
    return QuadraticPart{T}(
        convert(MOI.ScalarQuadraticFunction{T}, f.variable_part),
        convert(MOI.ScalarQuadraticFunction{T}, f.decision_part),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end
function Base.convert(::Type{QuadraticDecisionFunction{T,LQ}},
                      f::AffineDecisionFunction{T}) where {T, LQ <: LinearQuadraticPart{T}}
    return QuadraticDecisionFunction{T}(
        convert(LQ, f),
        convert(MOI.ScalarQuadraticFunction{T}, f.known_part),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
end

function Base.convert(::Type{QuadraticDecisionFunction{T,QuadraticPart{T}}},
                      f::QuadraticDecisionFunction{T,LinearPart{T}}) where T
    lq = f.linear_quadratic_terms
    return QuadraticDecisionFunction{T}(
        QuadraticPart(
            convert(MOI.ScalarQuadraticFunction{T}, lq.variable_part),
            convert(MOI.ScalarQuadraticFunction{T}, lq.decision_part),
            convert(MOI.ScalarQuadraticFunction{T}, zero(T))),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
end

function Base.convert(::Type{MOI.SingleVariable},
                      f::QuadraticDecisionFunction{T}) where T
    return convert(MOI.SingleVariable, convert(AffineDecisionFunction{T}, f))
end

function Base.convert(::Type{SingleDecision},
                      f::QuadraticDecisionFunction{T}) where T
    return convert(SingleDecision, convert(AffineDecisionFunction{T}, f.variable_part))
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::QuadraticDecisionFunction{T}) where T
    lq = f.linear_quadratic_terms
    if !Base.isempty(lq.cross_terms.quadratic_terms)
        throw(InexactError(:convert, AffineDecisionFunction{T}, f))
    end
    if !Base.isempty(f.known_cross_terms.quadratic_terms)
        throw(InexactError(:convert, AffineDecisionFunction{T}, f))
    end
    return AffineDecisionFunction{T}(
        convert(MOI.ScalarAffineFunction{T}, lq.variable_part),
        convert(MOI.ScalarAffineFunction{T}, lq.decision_part),
        convert(MOI.ScalarAffineFunction{T}, f.known_part))
end

function Base.convert(::Type{QuadraticDecisionFunction{T,LinearPart{T}}},
                      f::QuadraticDecisionFunction{T,QuadraticPart{T}}) where T
    lq = f.linear_quadratic_terms
    if !Base.isempty(lq.cross_terms.quadratic_terms)
        throw(InexactError(:convert, QuadraticDecisionFunction{T,LinearPart{T}}, f))
    end
    return QuadraticDecisionFunction(
        LinearPart(
            convert(MOI.ScalarAffineFunction{T}, lq.variable_part),
            convert(MOI.ScalarAffineFunction{T}, lq.decision_part)),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
end

function Base.zero(F::Type{QuadraticDecisionFunction{T,LQ}}) where {T, LQ <: LinearQuadraticPart{T}}
    return convert(F, zero(T))
end

function Base.iszero(f::QuadraticDecisionFunction)
    return iszero(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.isone(f::QuadraticDecisionFunction)
    return isone(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.one(F::Type{QuadraticDecisionFunction{T}}) where T
    return convert(F, one(T))
end

# JuMP overrides #
# ========================== #
function QuadraticDecisionFunction(quad::DQE)
    JuMP._assert_isfinite(quad)
    # Variable part
    variable_part = JuMP.moi_function(quad.variables)
    # Decision part
    decision_affine_terms =  Vector{MOI.ScalarAffineTerm{Float64}}()
    for (coef, dvar) in linear_terms(quad.decisions.aff)
        # Any fixed decision value is set in the decision bridge
        push!(decision_affine_terms, MOI.ScalarAffineTerm(coef, index(dvar)))
    end
    decision_quad_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for t in quad_terms(quad.decisions)
        # Any fixed decision value is set in the decision bridge
        push!(decision_quad_terms, JuMP._moi_quadratic_term(t))
    end
    # Cross terms
    cross_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for (cross_term, coeff) in quad.cross_terms
        # Any fixed decision value is set in the decision bridge
        push!(cross_terms,
              MOI.ScalarQuadraticTerm(coeff, index(cross_term.decision), index(cross_term.variable)))
    end
    # Known part
    known_affine_terms = Vector{MOI.ScalarAffineTerm{Float64}}()
    known_value = 0.0
    for (coef, kvar) in linear_terms(quad.knowns.aff)
        # Save the known terms to handle modifications in bridges
        push!(known_affine_terms, MOI.ScalarAffineTerm(coef, index(kvar)))
        # Then give the known values directly
        known_value += coef * value(kvar)
    end
    known_quad_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for t in quad_terms(quad.knowns)
        # Save the known terms to handle modifications in bridges
        push!(known_quad_terms, JuMP._moi_quadratic_term(t))
        # Then give the known values directly
        (coeff, var_1, var_2) = t
        known_value += coeff * value(var_1) * value(var_2)
    end
    # Known/variable cross terms
    known_variable_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for (cross_term, coeff) in quad.known_variable_terms
        # Any known decision is set in the known bridge
        push!(known_variable_terms,
              MOI.ScalarQuadraticTerm(coeff, index(cross_term.known), index(cross_term.variable)))
    end
    # Known/decision cross terms
    known_decision_terms = Vector{MOI.ScalarQuadraticTerm{Float64}}()
    for (cross_term, coeff) in quad.known_decision_terms
        # Any known decision is set in the known bridge
        push!(known_decision_terms,
              MOI.ScalarQuadraticTerm(coeff, index(cross_term.known), index(cross_term.variable)))
    end
    # Check if function has quadratic terms
    if isempty(variable_part.quadratic_terms) &&
        isempty(decision_quad_terms) &&
        isempty(cross_terms)
        # Create QuadraticDecisionFunction with LinearPart
        return QuadraticDecisionFunction(
            LinearPart(convert(MOI.ScalarAffineFunction{Float64}, variable_part),
                       MOI.ScalarAffineFunction{Float64}(decision_affine_terms, 0.0)),
            MOI.ScalarQuadraticFunction(known_affine_terms,
                                        known_quad_terms,
                                        known_value),
            MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                        known_variable_terms,
                                        0.0),
            MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                        known_decision_terms,
                                        0.0)
        )
    end
    # Return QuadraticDecisionFunction with QuadraticPart
    return QuadraticDecisionFunction(
        QuadraticPart(
            variable_part,
            MOI.ScalarQuadraticFunction(decision_affine_terms,
                                        decision_quad_terms,
                                        0.0),
            MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                        cross_terms,
                                        0.0)),
        MOI.ScalarQuadraticFunction(known_affine_terms,
                                    known_quad_terms,
                                    known_value),
        MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                    known_variable_terms,
                                    0.0),
        MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                    known_decision_terms,
                                    0.0))
end
JuMP.moi_function(quad::DQE) = QuadraticDecisionFunction(quad)
function JuMP.moi_function_type(::Type{DecisionQuadExpr{T}}) where T
    return QuadraticDecisionFunction{T,QuadraticPart{T}}
end

JuMP.moi_function(quad::BDE) = QuadraticDecisionFunction(quad)
function JuMP.moi_function_type(::Type{BilinearDecisionExpr{T}}) where T
    return QuadraticDecisionFunction{T,LinearPart{T}}
end

QuadraticDecisionFunction(quad::_DQE) = QuadraticDecisionFunction(convert(DQE, quad))
QuadraticDecisionFunction(quad::_KQE) = QuadraticDecisionFunction(convert(DQE, quad))

function _DecisionQuadExpr(m::Model, f::MOI.ScalarQuadraticFunction)
    quad = _DQE(_DecisionAffExpr(m, MOI.ScalarAffineFunction(f.affine_terms, 0.0)))
    for t in f.quadratic_terms
        v1 = t.variable_index_1
        v2 = t.variable_index_2
        coef = t.coefficient
        if v1 == v2
            coef /= 2
        end
        JuMP.add_to_expression!(quad, coef, DecisionRef(m, v1), DecisionRef(m, v2))
    end
    # There should not be any constants in the decision terms
    return quad
end

function _KnownQuadExpr(m::Model, f::MOI.ScalarQuadraticFunction)
    quad = _KQE(_KnownAffExpr(m, MOI.ScalarAffineFunction(f.affine_terms, 0.0)))
    for t in f.quadratic_terms
        v1 = t.variable_index_1
        v2 = t.variable_index_2
        coef = t.coefficient
        if v1 == v2
            coef /= 2
        end
        JuMP.add_to_expression!(quad, coef, KnownRef(m, v1), KnownRef(m, v2))
    end
    # There should not be any constants in the known terms
    return quad
end

function DQE(model::Model, f::QuadraticDecisionFunction{T}) where T
    lq = f.linear_quadratic_terms
    # Decision part
    cross_terms = OrderedDict{DecisionCrossTerm, Float64}()
    if lq isa QuadraticPart
        for term in lq.cross_terms.quadratic_terms
            cross_term = DecisionCrossTerm(
                DecisionRef(model, term.variable_index_1),
                VariableRef(model, term.variable_index_2))
            JuMP._add_or_set!(cross_terms, cross_term, term.coefficient)
        end
    end
    # Known/variable cross terms
    known_variable_terms = OrderedDict{KnownVariableCrossTerm, Float64}()
    for term in f.known_variable_terms.quadratic_terms
        cross_term = KnownVariableCrossTerm(
            KnownRef(model, term.variable_index_1),
            VariableRef(model, term.variable_index_2))
        JuMP._add_or_set!(known_variable_terms, cross_term, term.coefficient)
    end
    # Known/decision cross terms
    known_decision_terms = OrderedDict{KnownDecisionCrossTerm, Float64}()
    for term in f.known_decision_terms.quadratic_terms
        cross_term = KnownDecisionCrossTerm(
            KnownRef(model, term.variable_index_1),
            DecisionRef(model, term.variable_index_2))
        JuMP._add_or_set!(known_decision_terms, cross_term, term.coefficient)
    end
    return DQE(QuadExpr(model,
                        convert(MOI.ScalarQuadraticFunction{T}, lq.variable_part)),
               _DecisionQuadExpr(model,
                                convert(MOI.ScalarQuadraticFunction{T}, lq.decision_part)),
               _KnownQuadExpr(model, f.known_part),
               cross_terms,
               known_variable_terms,
               known_decision_terms)
end

function JuMP.jump_function_type(::Model,
                                 ::Type{QuadraticDecisionFunction{T,LinearPart{T}}}) where T
    return BilinearDecisionExpr{T}
end
function JuMP.jump_function_type(::Model,
                                 ::Type{QuadraticDecisionFunction{T,QuadraticPart{T}}}) where T
    return DecisionQuadExpr{T}
end
function JuMP.jump_function(model::Model, f::QuadraticDecisionFunction{T}) where T
    return DecisionQuadExpr{T}(model, f)
end

# MOI Function interface #
# ========================== #
MOI.constant(f::QuadraticDecisionFunction) =
    MOI.constant(f.linear_quadratic_terms.variable_part)
MOI.constant(f::QuadraticDecisionFunction, T::DataType) = MOI.constant(f)

function MOIU.eval_variables(varval::Function, f::LinearPart)
    var_value = MOIU.eval_variables(varval, f.variable_part)
    dvar_value = MOIU.eval_variables(varval, f.decision_part)
    return var_value + dvar_value
end
function MOIU.eval_variables(varval::Function, f::QuadraticPart)
    var_value = MOIU.eval_variables(varval, f.variable_part)
    dvar_value = MOIU.eval_variables(varval, f.decision_part)
    cross_value = MOIU.eval_variables(varval, f.cross_terms)
    return var_value + dvar_value + cross_value
end
function MOIU.eval_variables(varval::Function, f::QuadraticDecisionFunction)
    lq_value = MOIU.eval_variables(varval, f.linear_quadratic_terms)
    known_value = f.known_part.constant
    known_variable_value = MOIU.eval_variables(varval, f.known_variable_terms)
    known_decision_value = MOIU.eval_variables(varval, f.known_decision_terms)
    return lq_value + known_value + known_variable_value + known_decision_value
end

function MOIU.map_indices(index_map::Function, f::LinearPart{T}) where T
    return typeof(f)(MOIU.map_indices(index_map, f.variable_part),
                     MOIU.map_indices(index_map, f.decision_part))
end
function MOIU.map_indices(index_map::Function, f::QuadraticPart{T}) where T
    return typeof(f)(MOIU.map_indices(index_map, f.variable_part),
                     MOIU.map_indices(index_map, f.decision_part),
                     MOIU.map_indices(index_map, f.cross_terms))
end
function MOIU.map_indices(index_map::Function, f::QuadraticDecisionFunction{T}) where T
    return typeof(f)(MOIU.map_indices(index_map, f.linear_quadratic_terms),
                     MOIU.map_indices(index_map, f.known_part),
                     MOIU.map_indices(index_map, f.known_variable_terms),
                     MOIU.map_indices(index_map, f.known_decision_terms))
end

function MOIU.substitute_variables(variable_map::Function, f::QuadraticDecisionFunction{T, LinearPart{T}}) where T
    lq = f.linear_quadratic_terms
    g = QuadraticDecisionFunction(
        LinearPart(convert(MOI.ScalarAffineFunction{T}, MOI.constant(lq.variable_part)),
                   convert(MOI.ScalarAffineFunction{T}, zero(T))),
        MOI.ScalarQuadraticFunction{T}(
            f.known_part.affine_terms,
            f.known_part.quadratic_terms,
            zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
    g_lq = g.linear_quadratic_terms
    # Substitute variables
    for term in lq.variable_part.terms
        new_term = MOIU.substitute_variables(variable_map, term)
        MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarAffineFunction{T}
    end
    # Substitute decisions
    for term in lq.decision_part.terms
        new_term = MOIU.substitute_variables(variable_map, term)
        if iszero(new_term.constant) && term != new_term.terms[1]
            # If there no fixed value, the decision has been mapped to a variable
            MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarAffineFunction{T}
            # Otherwise, decision has been fixed and is added to decision part constant below
        end
        # Always keep the full term in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g_lq.decision_part, new_term)::MOI.ScalarAffineFunction{T}
    end
    # Substitute knowns
    _substitute_knowns(variable_map, f, g)
    return g
end

function MOIU.substitute_variables(variable_map::Function, f::QuadraticDecisionFunction{T, QuadraticPart{T}}) where T
    lq = f.linear_quadratic_terms
    g = QuadraticDecisionFunction(
        QuadraticPart(convert(MOI.ScalarQuadraticFunction{T}, MOI.constant(lq.variable_part)),
                      convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
                      convert(MOI.ScalarQuadraticFunction{T}, zero(T))),
        MOI.ScalarQuadraticFunction{T}(
            f.known_part.affine_terms,
            f.known_part.quadratic_terms,
            zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
    g_lq = g.linear_quadratic_terms
    # Substitute variables
    for term in lq.variable_part.affine_terms
        new_term = MOIU.substitute_variables(variable_map, term)
        MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    for term in lq.variable_part.quadratic_terms
        new_term = MOIU.substitute_variables(variable_map, term)
        MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    # Substitute decisions
    for term in lq.decision_part.affine_terms
        new_term = MOIU.substitute_variables(variable_map, term)
        if iszero(new_term.constant) && term != new_term.terms[1]
            # If there no fixed value, the decision has been mapped to a variable
            MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarQuadraticFunction{T}
            # Otherwise, decision has been fixed and is added to decision part constant below
        end
        # Always keep the full term in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g_lq.decision_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    for term in lq.decision_part.quadratic_terms
        new_term = MOIU.substitute_variables(variable_map, term)
        if iszero(new_term.constant) && term != new_term.quadratic_terms[1]
            # If there no fixed value, the decision has been mapped to a variable
            MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarQuadraticFunction{T}
            # Otherwise, decision has been fixed and is added to decision part constant below
        end
        # Remove any linear cross terms
        empty!(new_term.affine_terms)
        # Always keep the full term in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g_lq.decision_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    # Substitute decisions in cross terms
    for term in lq.cross_terms.quadratic_terms
        new_term = MOIU.substitute_variables(variable_map, term)
        if isempty(new_term.affine_terms)
            # If there is no fixed value, the decision in the cross term has been mapped to a variable
            MOIU.operate!(+, T, g_lq.variable_part, new_term)::MOI.ScalarQuadraticFunction{T}
        else
            # Decision has been fixed to a value,
            # giving an affine term, which can be
            # repeated in g already
            for term in new_term.affine_terms
                add_term!(g_lq.variable_part.affine_terms, term)
            end
        end
        # Always keep the cross terms in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g_lq.cross_terms, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    # Substitute knowns
    _substitute_knowns(variable_map, f, g)
    return g
end

function _substitute_knowns(variable_map::Function, f::QuadraticDecisionFunction{T}, g::QuadraticDecisionFunction{T}) where T
    lq = g.linear_quadratic_terms
    # Substitute known decisions
    for term in f.known_part.affine_terms
        # Known value acquired during substitution
        new_term = MOIU.substitute_variables(variable_map, term)
        # Add the known value to the constant of the known part
        MOIU.operate!(+, T, g.known_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    for term in f.known_part.quadratic_terms
        # Known value acquired during substitution
        new_term = MOIU.substitute_variables(variable_map, term)
        # Add the known value to the constant of the known part
        MOIU.operate!(+, T, g.known_part, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    # Substitute knowns in known/variable cross terms
    for term in f.known_variable_terms.quadratic_terms
        # First map only the variable part of the cross term
        known = MOI.SingleVariable(term.variable_index_1)
        mapped_variable = variable_map(term.variable_index_2)
        new_term = MOIU.operate(*, T, known, mapped_variable)::MOI.ScalarQuadraticFunction{T}
        MOIU.operate!(*, T, new_term, term.coefficient)
        # Just keep the mapped cross term to properly unbridge and handle modifications
        MOIU.operate!(+, T, g.known_variable_terms, new_term)::MOI.ScalarQuadraticFunction{T}
        # Now, perform the actual substitution
        new_term = MOIU.substitute_variables(variable_map, term)
        for term in new_term.affine_terms
            if lq isa LinearPart
                add_term!(lq.variable_part.terms, term)
            elseif lq isa QuadraticPart
                add_term!(lq.variable_part.affine_terms, term)
            end
        end
        # Always keep the full term in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g.known_variable_terms, new_term)::MOI.ScalarQuadraticFunction{T}
    end
    # Substitute knowns/decisions in known/decision cross terms
    for term in f.known_decision_terms.quadratic_terms
        # First map only the decision part of the cross term
        known = MOI.SingleVariable(term.variable_index_1)
        mapped_variable = convert(MOI.SingleVariable, variable_map(term.variable_index_2))
        new_term = MOIU.operate(*, T, known, mapped_variable)::MOI.ScalarQuadraticFunction{T}
        MOIU.operate!(*, T, new_term, term.coefficient)
        # Just keep the mapped cross term to properly unbridge and handle modifications
        MOIU.operate!(+, T, g.known_decision_terms, new_term)::MOI.ScalarQuadraticFunction{T}
        # Now, perform the actual substitution
        new_term = MOIU.substitute_variables(variable_map, term)::MOI.ScalarQuadraticFunction{T}
        if iszero(new_term.constant)
            for term in new_term.affine_terms
                # Decision has been mapped to a variable
                # giving an affine term, which can be
                # included in g already
                if lq isa LinearPart
                    add_term!(lq.variable_part.terms, term)
                    add_term!(lq.decision_part.terms, term)
                elseif lq isa QuadraticPart
                    add_term!(lq.variable_part.affine_terms, term)
                    add_term!(lq.decision_part.affine_terms, term)
                end
                # Otherwise, decision has been fixed and is added to the known/decision constant below
            end
        end
        # Always keep the full term in order to properly unbridge and handle modifications
        MOIU.operate!(+, T, g.known_decision_terms, new_term)::MOI.ScalarQuadraticFunction{T}
    end
end

MOIU.is_canonical(f::LinearPart) =
    MOIU.is_canonical(f.variable_part) &&
    MOIU.is_canonical(f.decision_part)
MOIU.is_canonical(f::QuadraticPart) =
    MOIU.is_canonical(f.variable_part) &&
    MOIU.is_canonical(f.decision_part) &&
    MOIU.is_canonical(f.cross_terms)
MOIU.is_canonical(f::QuadraticDecisionFunction) =
    MOIU.is_canonical(f.linear_quadratic_terms) &&
    MOIU.is_canonical(f.known_part) &&
    MOIU.is_canonical(f.known_variable_terms) &&
    MOIU.is_canonical(f.known_decision_terms)

MOIU.canonical(f::QuadraticDecisionFunction) = MOIU.canonicalize!(copy(f))

function MOIU.canonicalize!(f::LinearPart)
    MOIU.canonicalize!(f.variable_part)
    MOIU.canonicalize!(f.decision_part)
    return f
end
function MOIU.canonicalize!(f::QuadraticPart)
    MOIU.canonicalize!(f.variable_part)
    MOIU.canonicalize!(f.decision_part)
    MOIU.canonicalize!(f.cross_terms)
    return f
end
function MOIU.canonicalize!(f::QuadraticDecisionFunction)
    MOIU.canonicalize!(f.linear_quadratic_terms)
    MOIU.canonicalize!(f.known_part)
    MOIU.canonicalize!(f.known_variable_terms)
    MOIU.canonicalize!(f.known_decision_terms)
    return f
end

MOIU._is_constant(f::LinearPart) =
    MOIU._is_constant(f.variable_part) &&
    MOIU._is_constant(f.decision_part)
MOIU._is_constant(f::QuadraticPart) =
    MOIU._is_constant(f.variable_part) &&
    MOIU._is_constant(f.decision_part) &&
    MOIU._is_constant(f.cross_terms)
MOIU._is_constant(f::QuadraticDecisionFunction) =
    MOIU._is_constant(f.linear_quadratic_terms) &&
    MOIU._is_constant(f.known_part) &&
    MOIU._is_constant(f.known_variable_terms) &&
    MOIU._is_constant(f.known_decision_terms)

MOIU.all_coefficients(p::Function, f::LinearPart) =
    MOIU.all_coefficients(p, f.variable_part) &&
    MOIU.all_coefficients(p, f.decision_part)
MOIU.all_coefficients(p::Function, f::QuadraticPart) =
    MOIU.all_coefficients(p, f.variable_part) &&
    MOIU.all_coefficients(p, f.decision_part) &&
    MOIU.all_coefficients(p, f.cross_terms)
MOIU.all_coefficients(p::Function, f::QuadraticDecisionFunction) =
    MOIU.all_coefficients(p, f.linear_quadratic_terms) &&
    MOIU.all_coefficients(p, f.known_part) &&
    MOIU.all_coefficients(p, f.known_variable_terms) &&
    MOIU.all_coefficients(p, f.known_decision_terms)

MOIU.isapprox_zero(f::LinearPart, tol) =
    MOIU.isapprox_zero(f.variable_part, tol) &&
    MOIU.isapprox_zero(f.decision_part, tol)
MOIU.isapprox_zero(f::QuadraticPart, tol) =
    MOIU.isapprox_zero(f.variable_part, tol) &&
    MOIU.isapprox_zero(f.decision_part, tol) &&
    MOIU.isapprox_zero(f.cross_terms, tol)
MOIU.isapprox_zero(f::QuadraticDecisionFunction, tol) =
    MOIU.isapprox_zero(f.linear_quadratic_terms, tol) &&
    MOIU.isapprox_zero(f.known_part, tol) &&
    MOIU.isapprox_zero(f.known_variable_terms, tol) &&
    MOIU.isapprox_zero(f.known_decision_terms, tol)

function MOIU.filter_variables(keep::Function, f::LinearPart)
    return typeof(f)(MOIU.filter_variables(keep, f.variable_part),
                     MOIU.filter_variables(keep, f.decision_part))
end
function MOIU.filter_variables(keep::Function, f::QuadraticPart)
    return typeof(f)(MOIU.filter_variables(keep, f.variable_part),
                     MOIU.filter_variables(keep, f.decision_part),
                     MOIU.filter_variables(keep, f.cross_terms))
end
function MOIU.filter_variables(keep::Function, f::QuadraticDecisionFunction)
    # Only filter variable part, decision part and any cross terms
    return typeof(f)(MOIU.filter_variables(keep, f.linear_quadratic_terms),
                     copy(f.known_part),
                     MOIU.filter_variables(keep, f.known_variable_terms),
                     MOIU.filter_variables(keep, f.known_decision_terms))
end

function add_term!(terms::Vector{MOI.ScalarQuadraticTerm{T}},
                   term::MOI.ScalarQuadraticTerm{T}) where T
    index_1 = term.variable_index_1
    index_2 = term.variable_index_2
    coefficient = term.coefficient
    i = something(findfirst(t -> t.variable_index_1 == index_1 &&
                            t.variable_index_2 == index_2,
                            terms), 0)
    if iszero(i)
        if !iszero(coefficient)
            # The variable was not already included in the terms
            push!(terms, MOI.ScalarQuadraticTerm(coefficient, index_1, index_2))
        end
    else
        # Variable already included, increment the coefficient
        new_coeff = terms[i].coefficient + coefficient
        if iszero(new_coeff)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarQuadraticTerm(new_coeff, index_1, index_2)
        end
    end
    return nothing
end


text:
struct Prop2DAcoVTIDenQ_DEO2_FDTD
    p::Ptr{Cvoid}
end

function Prop2DAcoVTIDenQ_DEO2_FDTD(;
        nz=0,
        nx=0,
        nsponge=60,
        nbz=512,
        nbx=8,
        dz=5.0,
        dx=5.0,
        dt=0.001,
        nthreads=Sys.CPU_THREADS,
        freesurface=true,
        freqQ=5.0,
        qMin=0.1,
        qInterior=100.0)
    nz,nx,nsponge,nbz,nbx,nthreads = map(x->round(Int,x), (nz,nx,nsponge,nbz,nbx,nthreads))
    dz,dx,dt = map(x->Float32(x), (dz,dx,dt))

    @assert nx > 0
    @assert nz > 0
    @assert nsponge > 0
    @assert nthreads > 0
    @assert nbx > 0
    @assert nbz > 0

    fs = freesurface ? 1 : 0

    p = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_alloc, libprop2DAcoVTIDenQ_DEO2_FDTD),
        Ptr{Cvoid},
        (Clong, Clong,    Clong, Clong, Clong,   Cfloat, Cfloat, Cfloat, Clong, Clong),
         fs,    nthreads, nx,    nz,    nsponge, dx,     dz,     dt,     nbx,   nbz)

    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_SetupDtOmegaInvQ, libprop2DAcoVTIDenQ_DEO2_FDTD),
        Cvoid,
        (Ptr{Cvoid}, Cfloat, Cfloat, Cfloat),
         p,          freqQ,  qMin,   qInterior)

    Prop2DAcoVTIDenQ_DEO2_FDTD(p)
end

free(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_free, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

function size(prop::Prop2DAcoVTIDenQ_DEO2_FDTD)
    nz = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_getNz, libprop2DAcoVTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nx = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_getNx, libprop2DAcoVTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    (nz,nx)
end

for _f in (:V, :Eps, :Eta, :B, :F, :PSpace, :MSpace, :PCur, :POld, :MCur, :MOld, :DtOmegaInvQ)
    symf = "Prop2DAcoVTIDenQ_DEO2_FDTD_get" * string(_f)
    @eval $(_f)(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) = unsafe_wrap(Array, ccall(($symf, libprop2DAcoVTIDenQ_DEO2_FDTD), Ptr{Float32}, (Ptr{Cvoid},), prop.p), size(prop), own=false)
end

propagateforward!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_TimeStep, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)
propagateforward_linear!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_TimeStepLinear, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)
propagateadjoint!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) = propagateforward_linear!(prop) # self-adjoint

scale_spatial_derivatives!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD) =
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_ScaleSpatialDerivatives, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

abstract type Prop2DAcoVTIDenQ_DEO2_FDTD_Model end
struct Prop2DAcoVTIDenQ_DEO2_FDTD_Model_VEA <: Prop2DAcoVTIDenQ_DEO2_FDTD_Model end
struct Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V <: Prop2DAcoVTIDenQ_DEO2_FDTD_Model end

# v,ϵ and η in model-space
function forwardBornInjection!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_VEA, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_ForwardBornInjection_VEA, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat},        Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["ϵ"], dmodel["η"], wavefields["pold"], wavefields["mold"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_VEA,
        imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_AdjointBornAccumulation_VEA, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat},        Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["ϵ"], dmodel["η"], wavefields["pold"], wavefields["mold"], wavefields["pspace"], wavefields["mspace"])
end

# v in model-space
function forwardBornInjection!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_ForwardBornInjection_V, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_AdjointBornAccumulation_V, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionWaveFieldSeparationFWI, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep_V, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat},          Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"], 1)
end

function adjointBornAccumulation!(prop::Prop2DAcoVTIDenQ_DEO2_FDTD, modeltype::Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionWaveFieldSeparationRTM, dmodel, wavefields)
    ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep_V, libprop2DAcoVTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat},          Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"], 0)
end

function show(io::IO, prop::Prop2DAcoVTIDenQ_DEO2_FDTD)
    nx = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_getNx, libprop2DAcoVTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nz = ccall((:Prop2DAcoVTIDenQ_DEO2_FDTD_getNz, libprop2DAcoVTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    write(io, "Prop2DAcoVTIDenQ_DEO2_FDTD -- nx,nz; $nx,$nz")
end

ext: jl

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
include("type.jl")

function create_opf_model!(data::PowersenseData;
    formulation = PBRAPVmodel,
    initialize = true,
    FACTS_bShunt = true,
    FACTS_bSeries = true,
    box_constraints = false,
    obj_type = "linear")    #obj_type = "linear" or "quadratic"

    "JuMP Model for ACOPF formulations"
    data.model = JuMP.Model();

    #Variable initialization
    if initialize
        #Adding variables: modeling voltages for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel]
            JuMP.@variable(data.model, θ[i=1:data.nbus], start=data.θ[i]);
            JuMP.@variable(data.model, V[i=1:data.nbus], start=data.V[i]);
        elseif formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel, PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, vi[i=1:data.nbus], start=data.vi[i]);
            JuMP.@variable(data.model, vr[i=1:data.nbus], start=data.vr[i]);
        end
        if formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, Wd[i=1:data.nbus], start=data.Wd[i]);
            JuMP.@variable(data.model, Wr[i=1:data.nbr], start=data.Wr[i]);
            JuMP.@variable(data.model, Wi[i=1:data.nbr], start=data.Wi[i]);
        end

        #Adding variables: modeling generator injections for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel,
                            PBRARVmodel, CBRARVmodel, PNRARVmodel,
                            PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, pg[i=1:data.ngen], start=data.Pg[i]);
            JuMP.@variable(data.model, qg[i=1:data.ngen], start=data.Qg[i]);
        end

        #Adding variables: modeling branch flows for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PBRARVmodel, PBRAWVmodel]
            JuMP.@variable(data.model, Pij[i=1:data.nbr], start=data.Pij[i]);
            JuMP.@variable(data.model, Qij[i=1:data.nbr], start=data.Qij[i]);
            JuMP.@variable(data.model, Pji[i=1:data.nbr], start=data.Pji[i]);
            JuMP.@variable(data.model, Qji[i=1:data.nbr], start=data.Qji[i]);
        end
        if formulation ∈ [CBRARVmodel, CBRAWVmodel]
            JuMP.@variable(data.model, Irij[i=1:data.nbr], start=data.Irij[i]);
            JuMP.@variable(data.model, Iiij[i=1:data.nbr], start=data.Iiij[i]);
            JuMP.@variable(data.model, Irji[i=1:data.nbr], start=data.Irji[i]);
            JuMP.@variable(data.model, Iiji[i=1:data.nbr], start=data.Iiji[i]);
        end

        #Adding variables: modeling linearization using piecewise linear interpolation for different ACOPF formulations
        if obj_type == "linear"
            JuMP.@variable(data.model, cg[i=1:data.ngen], start=data.cg[i]);
        end

        #Adding variables: modeling shunt FACTS devices as continious switchable shunt succeptance
        if FACTS_bShunt
            JuMP.@variable(data.model, bss[i=1:data.nbss], start=data.Bss[i]);
        end

    else
        #Adding variables: modeling voltages for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel]
            JuMP.@variable(data.model, θ[1:data.nbus]);
            JuMP.@variable(data.model, V[1:data.nbus]);
        elseif formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel, PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, vi[1:data.nbus]);
            JuMP.@variable(data.model, vr[1:data.nbus]);
        end
        if formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, Wd[i=1:data.nbus]);
            JuMP.@variable(data.model, Wr[i=1:data.nbr]);
            JuMP.@variable(data.model, Wi[i=1:data.nbr]);
        end

        #Adding variables: modeling generator injections for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel,
                            PBRARVmodel, CBRARVmodel, PNRARVmodel,
                            PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@variable(data.model, pg[i=1:data.ngen]);
            JuMP.@variable(data.model, qg[i=1:data.ngen]);
        end

        #Adding variables: modeling branch flows for different ACOPF formulations
        if formulation ∈ [PBRAPVmodel, PBRARVmodel, PBRAWVmodel]
            JuMP.@variable(data.model, Pij[1:data.nbr]);
            JuMP.@variable(data.model, Qij[1:data.nbr]);
            JuMP.@variable(data.model, Pji[1:data.nbr]);
            JuMP.@variable(data.model, Qji[1:data.nbr]);
        end
        if formulation ∈ [CBRARVmodel, CBRAWVmodel]
            JuMP.@variable(data.model, Irij[1:data.nbr]);
            JuMP.@variable(data.model, Iiij[1:data.nbr]);
            JuMP.@variable(data.model, Irji[1:data.nbr]);
            JuMP.@variable(data.model, Iiji[1:data.nbr]);
        end

        #Adding variables: modeling linearization using piecewise linear interpolation for different ACOPF formulations
        if obj_type == "linear"
            JuMP.@variable(data.model, cg[i=1:data.ngen]);
        end

        #Adding variables: modeling shunt FACTS devices as continious switchable shunt succeptance
        if FACTS_bShunt
            JuMP.@variable(data.model, bss[i=1:data.nbss]);
        end
    end

    #Defining objective function
    if obj_type == "linear"
        JuMP.@variable(data.model, tgh[i = 1:data.ngen, j = 1:data.lps[i]] >= 0);
        JuMP.@objective(data.model, Min, sum(cg));
        for i=1:data.ngen
            JuMP.@constraint(data.model, cg[i] == sum(data.cgh[i][a] * tgh[i,a] for a = 1:data.lps[i]));
            JuMP.@constraint(data.model, pg[i] == sum(data.pgh[i][a] * tgh[i,a] for a = 1:data.lps[i]));
            JuMP.@constraint(data.model, sum(tgh[i,a] for a=1:data.lps[i]) == 1);
        end
    elseif obj_type == "quadratic"
        c0 = sum(data.c0);
    	if data.cost_order == 2
    		JuMP.@objective(data.model, Min, data.c2' * (pg .* pg) + (data.c1)' * pg + c0);
    	elseif data.cost_order == 3
    		JuMP.@objective(data.model, Min, (data.c1)' * pg + c0);
        else
            JuMP.@objective(data.model, Min, c0);
        end
    end

    #Adding constraints: defining generator injections bounds for different ACOPF formulations
    if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel,
                        PBRARVmodel, CBRARVmodel, PNRARVmodel,
                        PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
        JuMP.@constraint(data.model, data.Pmin .<= pg .<= data.Pmax);
        JuMP.@constraint(data.model, data.Qmin .<= qg .<= data.Qmax);
    end

    #Adding constraint: defining bounds for shunt FACTS devices modeled as continious switchable shunt succeptance
    if FACTS_bShunt
        JuMP.@constraint(data.model, data.bmin .<= bss .<= data.bmax);
    end

    #Adding constraints: defining voltages bounds for different ACOPF formulations
    if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel]
        JuMP.@constraint(data.model, data.Vmin .<= V .<= data.Vmax);
    elseif formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel, PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
        JuMP.@constraint(data.model, data.Vmin .* data.Vmin .<= vr .* vr + vi .* vi .<= data.Vmax .* data.Vmax);
    end

    if formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
        JuMP.@constraint(data.model, data.Vmin .* data.Vmin .<= Wd .<= data.Vmax .* data.Vmax);
        JuMP.@constraint(data.model, Wd .== vr .* vr + vi .* vi);
    end

    for i=1:data.nbus
        PNLexpr = JuMP.@NLexpression(data.model, 0.0);
        QNLexpr = JuMP.@NLexpression(data.model, 0.0);
        nbal = data.node[string(i)];

        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel,
            PBRARVmodel, CBRARVmodel, PNRARVmodel,
            PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            for g in nbal["gen"]
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr + pg[g]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + qg[g]);
            end
        end

        if formulation ∈ [PBRAPVmodel, PBRARVmodel, PBRAWVmodel]
            for k in nbal["Sij"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr + Pij[k]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + Qij[k]);
            end
            for k in nbal["Sji"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr + Pji[k]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + Qji[k]);
            end
        elseif formulation ∈ [CBRARVmodel, CBRAWVmodel]
            for k in nbal["Sij"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr + vr[f] * Irij[k] + vi[f] * Iiij[k]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + vi[f] * Irij[k] - vr[f] * Iiij[k]);
            end
            for k in nbal["Sji"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr + vr[t] * Irji[k] + vi[t] * Iiji[k]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + vi[t] * Irji[k] - vr[t] * Iiji[k]);
            end
        elseif formulation == PNPAPVmodel
            for k in nbal["Sij"]        f = data.br[k][1];                 t = data.br[k][2];
                Y = abs(data.Gf[k] + (data.Bf[k])im);         θY = angle(data.Gf[k] + (data.Bf[k])im);
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gf[k] * V[f]^2 - Y * cos(θ[f] - θ[t] - θY) * V[f] * V[t]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bf[k] * V[f]^2 - Y * sin(θ[f] - θ[t] - θY) * V[f] * V[t]);
            end
            for k in nbal["Sji"]        f = data.br[k][1];                 t = data.br[k][2];
                Y = abs(data.Gt[k] + (data.Bt[k])im);         θY = angle(data.Gt[k] + (data.Bt[k])im);
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gt[k] * V[t]^2 - Y * cos(θ[t] - θ[f] - θY) * V[f] * V[t]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bt[k] * V[t]^2 - Y * sin(θ[t] - θ[f] - θY) * V[f] * V[t]);
            end
        elseif formulation == PNRAPVmodel
            for k in nbal["Sij"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gf[k] * V[f]^2 - (data.Gf[k] * cos(θ[f] - θ[t]) + data.Bf[k] * sin(θ[f] - θ[t])) * V[f] * V[t]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bf[k] * V[f]^2 + (data.Bf[k] * cos(θ[f] - θ[t]) - data.Gf[k] * sin(θ[f] - θ[t])) * V[f] * V[t]);
            end
            for k in nbal["Sji"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gt[k] * V[t]^2 - (data.Gt[k] * cos(θ[t] - θ[f]) + data.Bt[k] * sin(θ[t] - θ[f])) * V[f] * V[t]);
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bt[k] * V[t]^2 + (data.Bt[k] * cos(θ[t] - θ[f]) - data.Gt[k] * sin(θ[t] - θ[f])) * V[f] * V[t]);
            end
        elseif formulation == PNRARVmodel
            for k in nbal["Sij"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gf[k] * (vr[f]^2+vi[f]^2) - data.Gf[k] * (vr[f]*vr[t]+vi[f]*vi[t]) + data.Bf[k] * (vr[f]*vi[t]-vi[f]*vr[t]));
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bf[k] * (vr[f]^2+vi[f]^2) + data.Bf[k] * (vr[f]*vr[t]+vi[f]*vi[t]) + data.Gf[k] * (vr[f]*vi[t]-vi[f]*vr[t]));
            end
            for k in nbal["Sji"]        f = data.br[k][1];                 t = data.br[k][2];
                PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.gt[k] * (vr[t]^2+vi[t]^2) - data.Gt[k] * (vr[t]*vr[f]+vi[t]*vi[f]) + data.Bt[k] * (vr[t]*vi[f]-vi[t]*vr[f]));
                QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.bt[k] * (vr[t]^2+vi[t]^2) + data.Bt[k] * (vr[t]*vr[f]+vi[t]*vi[f]) + data.Gt[k] * (vr[t]*vi[f]-vi[t]*vr[f]));
            end
        end

        if formulation ∈ [PBRAPVmodel, PNPAPVmodel, PNRAPVmodel]
            PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.Gl[i] * V[i]^2);
            QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.Bl[i] * V[i]^2);
            if FACTS_bShunt
                for sh in nbal["shunt"]
                    QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + bss[sh] * V[i]^2);
                end
            end
        elseif formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel]
            PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.Gl[i] * (vr[i]^2 + vi[i]^2));
            QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.Bl[i] * (vr[i]^2 + vi[i]^2));
            if FACTS_bShunt
                for sh in nbal["shunt"]
                    QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + bss[sh] * (vr[i]^2 + vi[i]^2));
                end
            end
        elseif formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            PNLexpr = JuMP.@NLexpression(data.model, PNLexpr - data.Gl[i] * Wd[i]);
            QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + data.Bl[i] * Wd[i]);
            if FACTS_bShunt
                for sh in nbal["shunt"]
                    QNLexpr = JuMP.@NLexpression(data.model, QNLexpr + bss[sh] * Wd[i]);
                end
            end
        end

        JuMP.@NLconstraint(data.model, PNLexpr == data.Pd[i]);
        JuMP.@NLconstraint(data.model, QNLexpr == data.Qd[i]);
    end

    for k=1:data.nbr                   f = data.br[k][1];                         t = data.br[k][2];
        if formulation == PBRAPVmodel
            JuMP.@NLconstraint(data.model, Pij[k] == - data.gf[k] * V[f]^2 - (data.Gf[k] * cos(θ[f] - θ[t]) + data.Bf[k] * sin(θ[f] - θ[t])) * V[f] * V[t]);
            JuMP.@NLconstraint(data.model, Pji[k] == - data.gt[k] * V[t]^2 - (data.Gt[k] * cos(θ[t] - θ[f]) + data.Bt[k] * sin(θ[t] - θ[f])) * V[f] * V[t]);
            JuMP.@NLconstraint(data.model, Qij[k] ==   data.bf[k] * V[f]^2 + (data.Bf[k] * cos(θ[f] - θ[t]) - data.Gf[k] * sin(θ[f] - θ[t])) * V[f] * V[t]);
            JuMP.@NLconstraint(data.model, Qji[k] ==   data.bt[k] * V[t]^2 + (data.Bt[k] * cos(θ[t] - θ[f]) - data.Gt[k] * sin(θ[t] - θ[f])) * V[f] * V[t]);
            JuMP.@NLconstraint(data.model, Pij[k]^2 + Qij[k]^2 <= data.Imax[k]^2 * V[f]^2);
            JuMP.@NLconstraint(data.model, Pji[k]^2 + Qji[k]^2 <= data.Imax[k]^2 * V[t]^2);
        elseif formulation == PBRARVmodel
            JuMP.@NLconstraint(data.model, Pij[k] == - data.gf[k] * (vr[f]^2+vi[f]^2) - data.Gf[k] * (vr[f]*vr[t]+vi[f]*vi[t]) + data.Bf[k] * (vr[f]*vi[t]-vi[f]*vr[t]));
            JuMP.@NLconstraint(data.model, Pji[k] == - data.gt[k] * (vr[t]^2+vi[t]^2) - data.Gt[k] * (vr[t]*vr[f]+vi[t]*vi[f]) + data.Bt[k] * (vr[t]*vi[f]-vi[t]*vr[f]));
            JuMP.@NLconstraint(data.model, Qij[k] == + data.bf[k] * (vr[f]^2+vi[f]^2) + data.Bf[k] * (vr[f]*vr[t]+vi[f]*vi[t]) + data.Gf[k] * (vr[f]*vi[t]-vi[f]*vr[t]));
            JuMP.@NLconstraint(data.model, Qji[k] == + data.bt[k] * (vr[t]^2+vi[t]^2) + data.Bt[k] * (vr[t]*vr[f]+vi[t]*vi[f]) + data.Gt[k] * (vr[t]*vi[f]-vi[t]*vr[f]));
            JuMP.@NLconstraint(data.model, Pij[k]^2 + Qij[k]^2 <= data.Imax[k]^2 * (vr[f]^2+vi[f]^2));
            JuMP.@NLconstraint(data.model, Pji[k]^2 + Qji[k]^2 <= data.Imax[k]^2 * (vr[t]^2+vi[t]^2));
        elseif formulation == CBRARVmodel
            JuMP.@constraint(data.model, Irij[k] == - data.gf[k] * vr[f] + data.bf[k] * vi[f] - data.Gf[k] * vr[t] + data.Bf[k] * vi[t]);
            JuMP.@constraint(data.model, Irji[k] == - data.gt[k] * vr[t] + data.bt[k] * vi[t] - data.Gt[k] * vr[f] + data.Bt[k] * vi[f]);
            JuMP.@constraint(data.model, Iiij[k] == - data.gf[k] * vi[f] - data.bf[k] * vr[f] - data.Gf[k] * vi[t] - data.Bf[k] * vr[t]);
            JuMP.@constraint(data.model, Iiji[k] == - data.gt[k] * vi[t] - data.bt[k] * vr[t] - data.Gt[k] * vi[f] - data.Bt[k] * vr[f]);
            JuMP.@NLconstraint(data.model, Irij[k]^2 + Iiij[k]^2 <= data.Imax[k]^2);
            JuMP.@NLconstraint(data.model, Irji[k]^2 + Iiji[k]^2 <= data.Imax[k]^2);
        elseif formulation == CBRAWVmodel
            JuMP.@constraint(data.model, Irij[k] == - data.gf[k] * vr[f] + data.bf[k] * vi[f] - data.Gf[k] * vr[t] + data.Bf[k] * vi[t]);
            JuMP.@constraint(data.model, Irji[k] == - data.gt[k] * vr[t] + data.bt[k] * vi[t] - data.Gt[k] * vr[f] + data.Bt[k] * vi[f]);
            JuMP.@constraint(data.model, Iiij[k] == - data.gf[k] * vi[f] - data.bf[k] * vr[f] - data.Gf[k] * vi[t] - data.Bf[k] * vr[t]);
            JuMP.@constraint(data.model, Iiji[k] == - data.gt[k] * vi[t] - data.bt[k] * vr[t] - data.Gt[k] * vi[f] - data.Bt[k] * vr[f]);
            JuMP.@NLconstraint(data.model, Irij[k]^2 + Iiij[k]^2 <= data.Imax[k]^2);
            JuMP.@NLconstraint(data.model, Irji[k]^2 + Iiji[k]^2 <= data.Imax[k]^2);
        elseif formulation ∈ [PNPAPVmodel, PNRAPVmodel]
            y¹ = abs2(data.gf[k] + (data.bf[k])im);       y² = abs2(data.Gf[k] + (data.Bf[k])im);       y³ = 2 * abs(data.gf[k] + (data.bf[k])im) * abs(data.Gf[k] + (data.Bf[k])im);
            θ¹ = angle(data.gf[k] + (data.bf[k])im);      θ² = angle(data.Gf[k] + (data.Bf[k])im);
            JuMP.@NLconstraint(data.model, y¹ * V[f]^2 + y² * V[t]^2 + y³ * V[f] * V[t] * cos(θ[f] - θ[t] + θ¹ - θ²) <= data.Imax[k]^2);
            y¹ = abs2(data.gt[k] + (data.bt[k])im);       y² = abs2(data.Gt[k] + (data.Bt[k])im);       y³ = 2 * abs(data.gt[k] + (data.bt[k])im) * abs(data.Gt[k] + (data.Bt[k])im);
            θ¹ = angle(data.gt[k] + (data.bt[k])im);      θ² = angle(data.Gt[k] + (data.Bt[k])im);
            JuMP.@NLconstraint(data.model, y¹ * V[t]^2 + y² * V[f]^2 + y³ * V[t] * V[f] * cos(θ[t] - θ[f] + θ¹ - θ²) <= data.Imax[k]^2);
        elseif formulation == PNRARVmodel
            y = abs2(data.gf[k] + (data.bf[k])im);                    Y = abs2(data.Gf[k] + (data.Bf[k])im);
            y¹ = 2 * (data.gf[k] * data.Gf[k] + data.bf[k] * data.Bf[k]);   y² = 2 * (data.bf[k] * data.Gf[k] - data.gf[k] * data.Bf[k]);
            JuMP.@NLconstraint(data.model, y¹ * (vr[f]*vr[t]+vi[f]*vi[t]) + y² * (vr[f]*vi[t]-vi[f]*vr[t]) + y * (vr[f]^2+vi[f]^2) + Y * (vr[t]^2+vi[t]^2) <= data.Imax[k]^2);
            y = abs2(data.gt[k] + (data.bt[k])im);                    Y = abs2(data.Gt[k] + (data.Bt[k])im);
            y¹ = 2 * (data.gt[k] * data.Gt[k] + data.bt[k] * data.Bt[k]);   y² = 2 * (data.bt[k] * data.Gt[k] - data.gt[k] * data.Bt[k]);
            JuMP.@NLconstraint(data.model, y¹ * (vr[t]*vr[f]+vi[t]*vi[f]) + y² * (vr[t]*vi[f]-vi[t]*vr[f]) + y * (vr[t]^2+vi[t]^2) + Y * (vr[f]^2+vi[f]^2) <= data.Imax[k]^2);
        end
    end

    if box_constraints
        # add_box_constraints(data, formulation)
        if formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel, PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            JuMP.@constraint(data.model, - data.Vmax .<= vi .<= data.Vmax);
            JuMP.@constraint(data.model, - data.Vmax .<= vr .<= data.Vmax);
        end
        if formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel]
            for k=1:data.nbr                   f = data.br[k][1];                         t = data.br[k][2];
                JuMP.@constraint(data.model, - data.Vmax[f] * data.Vmax[t] <= Wr[k] <= data.Vmax[f] * data.Vmax[t]);
                JuMP.@constraint(data.model, - data.Vmax[f] * data.Vmax[t] <= Wi[k] <= data.Vmax[f] * data.Vmax[t]);
            end
        end
        if formulation ∈ [PBRAPVmodel, PBRARVmodel, PBRAWVmodel]
            for k=1:data.nbr                   f = data.br[k][1];                         t = data.br[k][2];
                JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[f] <= Pij[k] <= data.Imax[k] * data.Vmax[f]);
                JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[t] <= Pji[k] <= data.Imax[k] * data.Vmax[t]);
                JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[f] <= Qij[k] <= data.Imax[k] * data.Vmax[f]);
                JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[t] <= Qji[k] <= data.Imax[k] * data.Vmax[t]);
            end
        elseif formulation ∈ [CBRARVmodel, CBRAWVmodel]
            JuMP.@constraint(data.model, - data.Imax .<= Irij .<= data.Imax);
            JuMP.@constraint(data.model, - data.Imax .<= Iiij .<= data.Imax);
            JuMP.@constraint(data.model, - data.Imax .<= Irji .<= data.Imax);
            JuMP.@constraint(data.model, - data.Imax .<= Iiji .<= data.Imax);
        end
    end
end

# function add_box_constraints(data::PowersenseData, formulation::ACOPF_fromulation; voltage = true, flow = true)
#     if formulation ∈ [PBRARVmodel, CBRARVmodel, PNRARVmodel, PBRAWVmodel, CBRAWVmodel, PNRAWVmodel] && voltage
#         JuMP.@constraint(data.model, data.Vmin .<= vi .<= data.Vmax);
#         JuMP.@constraint(data.model, data.Vmin .<= vr .<= data.Vmax);
#     end
#     if formulation ∈ [PBRAWVmodel, CBRAWVmodel, PNRAWVmodel] && voltage
#         for k=1:data.nbr                   f = data.br[k][1];                         t = data.br[k][2];
#             JuMP.@constraint(data.model, - data.Vmax[f] * data.Vmax[t] <= Wr[k] <= data.Vmax[f] * data.Vmax[t]);
#             JuMP.@constraint(data.model, - data.Vmax[f] * data.Vmax[t] <= Wi[k] <= data.Vmax[f] * data.Vmax[t]);
#         end
#     end
#     if formulation ∈ [PBRAPVmodel, PBRARVmodel, PBRAWVmodel] && flow
#         for k=1:data.nbr                   f = data.br[k][1];                         t = data.br[k][2];
#             JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[f] <= Pij[k] <= data.Imax[k] * data.Vmax[f]);
#             JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[t] <= Pji[k] <= data.Imax[k] * data.Vmax[t]);
#             JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[f] <= Qij[k] <= data.Imax[k] * data.Vmax[f]);
#             JuMP.@constraint(data.model, - data.Imax[k] * data.Vmax[t] <= Qji[k] <= data.Imax[k] * data.Vmax[t]);
#         end
#     elseif formulation ∈ [CBRARVmodel, CBRAWVmodel] && flow
#         JuMP.@constraint(data.model, - data.Imax .<= Irij .<= data.Imax);
#         JuMP.@constraint(data.model, - data.Imax .<= Iiij .<= data.Imax);
#         JuMP.@constraint(data.model, - data.Imax .<= Irji .<= data.Imax);
#         JuMP.@constraint(data.model, - data.Imax .<= Iiji .<= data.Imax);
#     end
# end

function run_opf!(data::PowersenseData;
    formulation = PBRAPVmodel,
    create_model = true,
    solver = Powersense.Optimizer,
    initialize = true,
    FACTS_bShunt = true,
    FACTS_bSeries = false,
    box_constraints = false,
    solve = true,
    obj_type = "linear")    #obj_type = "linear" or "quadratic"

    if create_model
    	create_opf_model!(data, formulation = formulation, initialize = initialize, FACTS_bShunt = FACTS_bShunt,
    				FACTS_bSeries = FACTS_bSeries, box_constraints = box_constraints, obj_type = obj_type);
    end

    if solve
    	JuMP.set_optimizer(data.model, solver)
    	JuMP.optimize!(data.model);
    	if JuMP.has_values(data.model)
    		data.cost = JuMP.objective_value(data.model)
    	end
    end
end


text:
struct RawTime
    name::String
    stages::Dict{Period, Int}
    col_delims::Vector{Col}
    row_delims::Vector{Row}
end

function parse_tim(filename::String)
    # Initialize auxiliary variables
    name       = "SLP"
    mode       = :NAME
    stage      = 1
    stages     = Dict{Period, Int}()
    col_delims = Vector{Col}()
    row_delims = Vector{Row}()
    # Parse the file
    open(filename) do io
        firstline = split(readline(io))
        if Symbol(firstline[1]) == :TIME
            name = join(firstline[2:end], " ")
        else
            throw(ArgumentError("`TIME` field is expected on the first line."))
        end
        for line in eachline(io)
            if mode == END
                # Parse finished
                break
            end
            words = split(line)
            first_word = Symbol(words[1])
            if first_word in TIM_MODES
                mode = first_word
                continue
            end
            if mode == PERIODS
                push!(col_delims, Symbol(words[1]))
                push!(row_delims, Symbol(words[2]))
                stages[Symbol(words[3])] = stage
                stage += 1
            else
                throw(ArgumentError("$(mode) is not a valid tim file mode."))
            end
        end
    end
    # Return tim data
    return RawTime(name, stages, col_delims, row_delims)
end

num_stages(tim::RawTime) = length(tim.stages)


text:
using TropicalSemiring
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@test Trop(2.3) isa Trop{Max, Float64}
@test Trop(2) isa Trop{Max, Int64}
@test Trop{Min}(2) isa Trop{Min, Int64}
@test Trop{Min, Float64}(2) isa Trop{Min, Float64}

@test Trop{Max}(2.0, true) isa Trop{Max, Float64}
@test isinf(Trop{Max}(2.0, true))
@test isinf(inf(Max))
@test isinf(inf(Min))
@test inf() isa Trop{Max, Bool}
@test inf(Min) isa Trop{Min, Bool}
@test inf(Max) isa Trop{Max, Bool}

@test Trop(2) + Trop(3) == Trop(3)
@test Trop{Min}(2) + Trop{Min}(3) == Trop{Min}(2)
@test Trop{Min}(2) + inf(Min) == inf(Min)
@test Trop{Max}(2) + inf(Max) == inf(Max)

@test Trop{Min}(2) * Trop{Min}(3) == Trop{Min}(5)
@test Trop{Max}(2) * Trop{Max}(3) == Trop{Max}(5)
@test Trop{Max}(2) * inf() == inf()
@test Trop{Min}(2) * inf(Min) == inf(Min)

@test Trop(4)^4 == Trop(16)
@test -Trop(4) == Trop(4)

@test (Trop(3.1) == Trop(3.0)) == false
@test (Trop{Max}(3) == Trop{Min}(3)) == false

@test one(typeof(Trop(2))) == Trop(0)
@test one(Trop(3.0)) == Trop(0.0)

@test zero(typeof(Trop(3))) == Trop{Max}(zero(3), true)
@test zero(Trop(2.3)) == Trop{Max}(zero(3.2), true)

@test string(Trop(5.2)) == "5.2"
@test string(inf(Min)) == "-∞"
@test string(inf(Max)) == "∞"

@test promote_type(Trop{Max, Int}, Trop{Max, Float64}) == Trop{Max, Float64}


text:
using PyPlot, DelimitedFiles, Random, Statistics, Revise,
      transD_GP
## model fixed parts, i.e., air
Random.seed!(23)
zfixed   = [-1e5]
ρfixed   = [1e12]
nmax = 200
# Note that the receiver and transmitter need to be in layer 1
zstart = 0.0
extendfrac, dz = 1.03, 1.5
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=65)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
calcjacobian = true
##  geometry and modeling parameters
rRx = 13.
zRx = -42.0
zTx = -40.0
freqlow = 1e-3
include("electronics_halt.jl")
## LM operator
Flm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      lowpassfcs = lowpassfcs,
                      times  = LM_times,
                      ramp   = LM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      freqlow = freqlow,
                      calcjacobian = calcjacobian
                      )
## HM operator
Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      lowpassfcs = lowpassfcs,
                      times  = HM_times,
                      ramp   = HM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      freqlow = freqlow,
                      calcjacobian = calcjacobian
                      )
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<50)] .= 20.
ρ[(z.>=50) .&(z.<80)] .= 1
ρ[(z.>=80) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<200)] .= 50
ρ[(z.>=200) .&(z.<250)] .= 80
ρ[(z.>=250)]            .= 150
# add jitter to model in log10 domain
Random.seed!(11)
ρ = 10 .^(0.1*randn(length(ρ)) + log10.(ρ))
## add noise to data
transD_GP.plotmodelfield_skytem!(Flm, Fhm, z, ρ)
#remember to specify halt noise if it is in pV
dlow, dhigh, σlow, σhigh = transD_GP.addnoise_skytem(Flm, Fhm,
                z, ρ, noisefrac=0.03, halt_LM = LM_noise*1e-12, halt_HM = HM_noise*1e-12, rseed=21,
                dz=dz, extendfrac=extendfrac, nfixed=nfixed)
## create operator
# but with a coarser grid
extendfrac, dz = 1.06, 1.15
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=50)
zgrid, ρgrid, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
dlow, dhigh, σlow, σhigh = (dlow, dhigh, σlow, σhigh)./transD_GP.SkyTEM1DInversion.μ₀
aem = transD_GP.dBzdt(Flm, Fhm, dlow, dhigh, σlow, σhigh, z=zgrid, ρ=ρgrid, nfixed=nfixed)
## let's try gradient descent, all model values are in log10 conductivity
σstart, σ0 = map(x->zeros(length(aem.ρ)-1), 1:2)
σstart .= -2.
σ0 .= -2
regtype = :R1
## do it
m, χ², λ², idx = transD_GP.gradientinv(σstart, σ0, aem, nstepsmax=20,
                            regularizeupdate=false,
                            λ²min=0,
                            λ²max=7,
                            β² = 0.,
                            ntries=6,
                            knownvalue=0.7,
                            λ²frac=4,
                            αfrac=4,
                            dobo = true,
                            breakonknown = true,
                            firstvalue=:last,
                            κ = transD_GP.GP.Mat52(),
                            regtype = regtype);
## debug plots: all in each
alpha = 0.1
idxlast=length(m)
for (i, mi) in enumerate(m)
    if isempty(λ²[i])
        global idxlast = i-1
        break
    end
    f = figure(figsize=(7,6))
    s1 = subplot(131)
    for ii in 1:length(χ²[i])
        s1.step(-mi[ii], aem.z[2:end], alpha=alpha)
    end
    s1.step(-mi[idx[i]], aem.z[2:end], color="r", linewidth=2)
    s1.step(log10.(ρ[2:end]), z[2:end], color="k", linewidth=2, linestyle="--")
    s1.set_ylabel("depth m")
    s1.set_xlabel("Log₁₀ ρ")
    s1.invert_yaxis()
    s1.invert_xaxis()
    s2 = subplot(232)
    hps = hcat(λ²[i]...)'
    sortedλidx = sortperm(hps[:,1])
    sortedαidx = sortperm(hps[:,2])
    s2.loglog(hps[sortedλidx,1], χ²[i][sortedλidx], ".-", markersize=5)
    s2.plot(hps[idx[i],1], χ²[i][idx[i]], "o" )
    s2.set_ylabel("χ²")
    s3 = subplot(236)
    s3.loglog(χ²[i][sortedαidx], hps[sortedαidx,2], ".-", markersize=5)
    s3.plot(χ²[i][idx[i]], hps[idx[i],2], "o" )
    s3.set_xlabel("χ²")
    s4 = subplot(235, sharex=s2, sharey=s3)
    s4.scatter(hps[:,1], hps[:,2], c=log10.(χ²[i]), cmap="magma")
    s4.plot(hps[idx[i],1], hps[idx[i],2], "x", markersize=10)
    s4.set_xlabel("λ²"); s4.set_ylabel("step length α")
    plt.tight_layout()
end
## debug plots best in each
figure(figsize=(3,6))
alpha = 0.1
for (i, mi) in enumerate(m)
    step(-mi[idx[i]], aem.z[2:end], alpha=alpha)
    i == idxlast && break
end
step(-m[idxlast][idx[idxlast]], aem.z[2:end], color="r", linewidth=2)
step(log10.(ρ[2:end]), z[2:end], color="k", linewidth=2, linestyle="--")
gca().invert_yaxis()
gca().invert_xaxis()
ylabel("depth m")
xlabel("Log₁₀ ρ")
plt.tight_layout()
ax = gca()
## Compare with posterior model covariance
using LinearAlgebra
F = aem
R = transD_GP.makereg(regtype, F)
JtW, Wr = F.J'*F.W, F.W*F.res
H = JtW*(JtW)' + λ²[idxlast][idx[idxlast]][1]*R'R
Cpost = inv(Hermitian(H))
figure()
zplot = [zboundaries; zboundaries[end] + 5]
pcolormesh(zplot, zplot, Cpost)
xlabel("Depth m")
ylabel("Depth m")
colorbar()
sd = sqrt.(diag(Cpost))
ax.fill_betweenx(zall, -m[idxlast][idx[idxlast]] -sd, -m[idxlast][idx[idxlast]] +sd, alpha=0.2)
ax.set_xlim(3, -1)


text:
# start moving exports here and not all in IIF.jl
export blockCliqSiblingsParentNeedDown

function isCliqInitialized(cliq::TreeClique)::Bool
  return getCliqueData(cliq).initialized in [INITIALIZED; UPSOLVED]
end

function isCliqUpSolved(cliq::TreeClique)::Bool
  return getCliqueData(cliq).initialized == UPSOLVED
end

"""
    $SIGNATURES

Return the most likely  ordering for initializing factor (assuming up solve
sequence).

Notes:
- sorts id (label) for increasing number of connected factors using the clique subfg with messages already included.
"""
function getCliqVarInitOrderUp(subfg::AbstractDFG)
  # rules to explore dimension from one to the other?

  # get all variable ids and number of associated factors
  B, varLabels, facLabels = getBiadjacencyMatrix(subfg)
  nfcts = sum(B, dims=1)[:]

  # variables with priors
  varswithpriors = getNeighbors.(subfg, lsfPriors(subfg))
  singids = union(Symbol[], varswithpriors...)

  # sort permutation order for increasing number of factor association
  nfctsp = sortperm(nfcts)
  sortedids = varLabels[nfctsp]

  # organize the prior variables separately with asceding factor count
  initorder = intersect(sortedids, singids)
  # in ascending order of number of factors
  union!(initorder, sortedids)

  return initorder
end

"""
    $SIGNATURES

Special function to do initialization in downward direction, assuming that not all
variables can be initialized.  Relies on outside down messages.

Notes:
- assumed this `cliq` is being initialized from a previous `:needdownmsg` status.
- will use all possible local factors of cliquq in initilization process
- similar to upward initialization, but uses different message structure
  - first draft assumes upward messages will not be used,
  - full up solve still required which explicitly depends on upward messages.

Dev Notes
- Streamline add/delete msg priors from calling function and csm.
- TODO replace with nested 'minimum degree' type variable ordering.
"""
function getCliqInitVarOrderDown( dfg::AbstractDFG,
                                  cliq::TreeClique,
                                  dwnkeys::Vector{Symbol} )   # downmsgs
  #
  allsyms = getCliqAllVarIds(cliq)
  # convert input downmsg var symbols to integers (also assumed as prior beliefs)
  # make sure ids are in the clique set, since parent may have more variables.
  dwnmsgsym = intersect(dwnkeys, DFG.listVariables(dfg))
  dwnvarids = intersect(allsyms, dwnmsgsym)

  # find any other prior factors (might have partials)
  prvarids = getCliqVarIdsPriors(cliq, allsyms, true)
  hassinglids = union(dwnvarids, prvarids)

  # Get all other variable factor counts
  nfcts = getCliqNumAssocFactorsPerVar(cliq)
  # add msg marginal prior (singletons) to number of factors
  for msid in dwnmsgsym
    nfcts[msid .== allsyms] .+= 1
  end

  # sort permutation order for increasing number of factor association
  nfctsp = sortperm(nfcts)
  sortedids = allsyms[nfctsp]

  # all singleton variables
  singids = union(prvarids, dwnvarids)

  # organize the prior variables separately with asceding factor count
  initorder = Symbol[] #zeros(Int, 0)
  for id in sortedids
    if id in singids
      push!(initorder, id)
    end
  end
  # sort remaining variables for increasing associated factors
  for id in sortedids
    if !(id in initorder)
      push!(initorder, id)
    end
  end

  # return variable order
  return initorder::Vector{Symbol}
end

function _isInitializedOrInitSolveKey(var::DFGVariable, solveKey::Symbol=:default; N::Int=100)
  # TODO, this solveKey existence test should probably be removed?
  if !(solveKey in listSolveKeys(var))
    varType = getVariableType(var)
    setDefaultNodeData!(var, 0, N, getDimension(varType), solveKey=solveKey,
                        initialized=false, varType=varType, dontmargin=false)
    #
    # data = getSolverData(var, solveKey)
    # if data === nothing
    # end
    return false
  end
  # regular is initialized check, this is fine
  isinit = isInitialized(var, solveKey)
  return isinit
end

"""
    $SIGNATURES

Return true if clique has completed the local upward direction inference procedure.
"""
isUpInferenceComplete(cliq::TreeClique) = getCliqueData(cliq).upsolved

function areCliqVariablesAllInitialized(dfg::AbstractDFG,
                                        cliq::TreeClique,
                                        solveKey::Symbol=:default;
                                        N::Int=getSolverParams(dfg).N  )
  #
  allids = getCliqAllVarIds(cliq)
  isallinit = true
  for vid in allids
    var = DFG.getVariable(dfg, vid)
    isallinit &= _isInitializedOrInitSolveKey(var, solveKey, N=N)
    # isallinit &= isInitialized(var, solveKey)
  end
  isallinit
end

"""
    $SIGNATURES

Return true if all variables in clique are considered marginalized (and initialized).
"""
function areCliqVariablesAllMarginalized( subfg::AbstractDFG,
                                          cliq::TreeClique)
  for vsym in getCliqAllVarIds(cliq)
    vert = getVariable(subfg, vsym)
    if !isMarginalized(vert) || !isInitialized(vert)
      return false
    end
  end
  return true
end

function printCliqInitPartialInfo(subfg, cliq, solveKey::Symbol=:default, logger=ConsoleLogger())
  varids = getCliqAllVarIds(cliq)
  initstatus = Vector{Bool}(undef, length(varids))
  initpartial = Vector{Float64}(undef, length(varids))
  for i in 1:length(varids)
    initstatus[i] = isInitialized(subfg, varids[i], solveKey) # getSolverData(getVariable(subfg, varids[i]), solveKey).initialized
    initpartial[i] = -1 # getSolverData(getVariable(subfg, varids[i]), solveKey).inferdim
  end
  with_logger(logger) do
    tt = split(string(now()),'T')[end]
    @info "$tt, cliq $(cliq.id), PARINIT: $varids | $initstatus | $initpartial"
  end
end

#


text:
# Do Gauss quadrature by hand
# The routine that computes the weights are at
# https://github.com/JuliaLang/julia/blob/master/base/quadgk.jl
# note that the method is essentially the same as the one used
# in the Matlab example form class, but takes advantage of
# the special sparsity pattern and allows these points to be
# computed in arbitrary spaces
(gx,gw) = Base.gauss(Float64, 7)
do_quadrature(f,a,b) = sum([ w*f(x) for (w,x) in zip(gw,gx) ])
g(x) = cos(π*x)
data = Array(Array{Float64,1},4)
for p=1:4
    println("p: $p")
    # Pick the polynomial basis
    cheb = AC274.nodes(p+1)
    basis = AC274.phi(cheb)
    error = zeros(Float64,5)
    i=1
    for K = 10:10:50
        println("K: $K")
        # Regular Mesh of K cells
        mesh = generateMesh(Mesh1D,0.0,2.0,K)

        for cell in mesh
            # Interpolate at the nodes
            ℐg = AC274.evaluate(AC274.Coeffs{p+1}(map(n->g(𝜒⁻¹(mesh,cell,n)),cheb)),basis)

            # Do the quadrature on the current node

            error[i] += 0.5*(cell.right-cell.left)*AC274.integrate(y->(AC274.fapply(g,𝜒⁻¹(mesh,cell,y))-AC274.fapply(ℐg,y))^2,-1,1)
        end
        error[i] = sqrt(error[i])
        i += 1
    end
    data[p] = error
    #xs = 2./[10:10:50]
    #display(MIME("text/html"),plot(layer(h->h^(p+1),minimum(xs),maximum(xs)),
    #    layer(x=xs,y=error, Geom.point),Scale.x_log10,Scale.y_log10))
end


text:
"""
    mt_state(x) -> Bool
    mt_state!(x, val)

Specifies whether the magnetization transfer pulse is applied. This parameter
is REQUIRED by all the anatomical images grouped by MTR, MTS and MPM suffixes.
"""
@defprop MTState{:mt_state}::Bool

"""
    mt_offset_frequency(x)
    mt_offset_frequency!(x, val)

The frequency offset of the magnetization transfer pulse with respect to the
central H1 Larmor frequency in Hertz (Hz).
"""
@defprop MTOffsetFrequency{:mt_offset_frequency}::F64Hertz

"""
    mt_pulse_bandwidth(x)
    mt_pulse_bandwidth!(x, val)

The excitation bandwidth of the magnetization transfer pulse in Hertz (Hz).
"""
@defprop MTPulseBandwidth{:mt_pulse_bandwidth}::F64Hertz

"""
    mt_npulses(x)
    mt_npulses!(x, val)

Number of magnetization transfer RF pulses applied before the readout.
"""
@defprop MTNumberOfPulses{:mt_npulses}::Int

# TODO how should this be parameterized
"""
    mt_pulse_shape(x)
    mt_pulse_shape!(x, val)

Shape of the magnetization transfer RF pulse waveform. Accepted values:
* HARD
* GAUSSIAN
* GAUSSHANN (gaussian pulse with Hanning window)
* SINC
* SINCHANN (sinc pulse with Hanning window)
* SINCGAUSS (sinc pulse with Gaussian window)
* FERMI
"""
@defprop MTPulseShape{:mt_pulse_shape}

"""
    mt_pulse_duration(x)
    mt_pulse_duration!(x, val)

Duration of the magnetization transfer RF pulse in seconds.
"""
@defprop MTPulseDuration{:mt_pulse_duration}::F64Second

"""
    MagnetizationTransferMetadata

Metadata structure for information concerning MRI magnetization transfer pulse.

## Supported Properties

$(GETPROPERTY)

"""
struct MagnetizationTransferMetadata{H,S}
    mt_state::Bool
    mt_offset_frequency::H
    mt_pulse_bandwidth::H
    mt_npulses::Int
    mt_pulse_shape::Symbol
    mt_pulse_duration::S
end

@properties MagnetizationTransferMetadata begin
    mt_state(self) => :mt_state
    mt_offset_frequency(self) => :mt_offset_frequency
    mt_pulse_bandwidth(self) => :mt_pulse_bandwidth
    mt_npulses(self) => :mt_npulses
    mt_pulse_shape(self) => :mt_pulse_shape
    mt_pulse_duration(self) => :mt_pulse_duration
end


text:
@testset "dynamic programming" begin

    nominal_speed = 0.1
    walk_noise = 0.2
    prob_lag = 0.2
    prob_skip = 0.2
    prob_normal = 0.6
    walk_noise = prob_lag + prob_skip
    noise = 1.0
    obs_params = ObsModelParams(nominal_speed, walk_noise, noise)
    likelihood(a, b) = exp(Gen2DAgentMotion.noise_log_likelihood(obs_params, a, b))
    A = Point(0,0)
    B = Point(1,1)
    C = Point(2,2)
    obs1 = Point(0,0)
    obs2 = Point(1,1)
    (actual, _) = Gen2DAgentMotion.run_forward_backward(obs_params, [A, B, C], [obs1, obs2])
    first_prob_normal = prob_lag + prob_normal
    first_prob_skip = prob_skip
    expected = 0.0
    # 1-A, 2-A
    expected += first_prob_normal * prob_lag * likelihood(A, obs1) * likelihood(A, obs2)
    # 1-A, 2-B
    expected += first_prob_normal * prob_normal * likelihood(A, obs1) * likelihood(B, obs2)
    # 1-A, 2-C
    expected += first_prob_normal * prob_skip * likelihood(A, obs1) * likelihood(C, obs2)
    # 1-B, 2-B
    expected += first_prob_skip * prob_lag * likelihood(B, obs1) * likelihood(B, obs2)
    # 1-B, 2-C
    expected += first_prob_skip * prob_normal * likelihood(B, obs1) * likelihood(C, obs2)
    @test isapprox(actual, log(expected))

end


text:
module AdaMeanShift

    using Statistics, LinearAlgebra, SpecialFunctions, StaticArrays
    export atomic_meanshift, atomic_intensity, atomic_ratio, meanshift!, meanshift_nonadaptive!, extract_circ_region, extract_region
    @inline function Kern(x::AbstractVector{fT},h::AbstractVector{fT}) where fT
        max(zero(fT),one(fT) - norm(@inbounds x ./ max.(h,one(fT)))^2)
    end
    function clampdim(M,c::T,i)::T where T
        clamp(c,1,size(M,i))
    end
    function makerange(M,p,h,i)
        #round(Int,p[i]-h[i]-1):round(Int,p[i]+h[i]+1)
        @inbounds clampdim(M,round(Int,p[i]-h[i]-1),i):clampdim(M,round(Int,p[i]+h[i]+1),i)
    end

"""
        atomic_meanshift(M,p,h,isotropy,smoothing) -> NamedTuple(pnew,hnew,modeval,delta)

Performs a single iteration of meanshift over a single particle

# Arguments
- `M` density tensor
- `p` initial position for test particle
- `h` initial particle standard deviation
- `isotropy` scalar 0 -> 1, 0 is completely anisotropic, completely isotropic
- `smoothing` scalar 0 -> 1, regularization term for noisy density tensors

"""
    function atomic_meanshift(M,p,h,isotropy,smoothing)
        N=length(p)
        T=eltype(p)
        @assert N==length(h)
        @assert T==eltype(h)
        atomic_meanshift(M,SVector{length(p)}(p),SVector{length(h)}(h),T(isotropy),T(smoothing))
    end
    @generated function atomic_meanshift(M::Mty,p::K,h::K,isotropy::T,smoothing::T) where {T,N,K<:StaticArray,Mty<:AbstractArray{T,N}}
        D=(8//N^2 + 6//N +1)
        scale= sqrt(N * ((2//N + 1)*(4//N + 1)*(6//N + 1))//D)
        invvolume=gamma(N/2+1)*pi^(-N/2)
        modescale=invvolume*N^3*D
        Z=zero(T)
        O=one(T)
        quote
            @inbounds begin
                c=p.*$Z
                xsq=p.*$Z
                NW=IW=KW=W=$Z
                ls=$invvolume*inv(prod(max.(h,$O)))*smoothing

                Base.Cartesian.@nexprs $N i -> R_i = makerange(M,p,h,i)
                Base.Cartesian.@nloops $N i i->R_i begin
                    cp=$K(Base.Cartesian.@ntuple($N,k->i_k))
                    diff=p.-cp
                    keW=Kern(diff,h)
                    imW=Base.Cartesian.@nref($N,M,k->i_k) + ls
                    sW=keW*imW
                    isofactor = isotropy*norm(diff)^2/$N
                    xsq=xsq .+ (diff .^ 2 .*(1.0-isotropy) .+isofactor) .* sW
                    c=c.+cp.*sW
                    W+=sW
                    KW+=keW
                    IW+=imW
                    NW+=$O
                end
                if 2*NW<1
                    return (pnew=p,
                            hnew=h,
                            modeval=T(NaN),
                            delta=$Z)
                end
                xsq=xsq ./ W

                c=@. p+(c/W - p)/2
                return (pnew=K(Base.Cartesian.@ntuple($N,k->clampdim(M,c[k],k))),
                        hnew=max.(1,$scale .* sqrt.(xsq)),
                        modeval=W/prod(max.(h,$O))*$modescale,
                        delta=norm(c.-p))
            end
        end
    end

"""
        atomic_intensity(M,p,h) -> NamedTuple(intensity,numpoints,ndims)

Calculates the integrated intensity on the density tensor `M` over a ellipsoidal region caracterized by `p` region center, `h` semiaxis

# Arguments
- `M` density tensor (e.g. matrix N*N)
- `p` region center (e.g. 2D vector)
- `h` semiaxis lengths (e.g. 2D vector)

Note that both `p`,`h` should be of type Vector or StaticVector with the same length.
"""
    function atomic_intensity(M,p,h)
        N=length(p)
        T=eltype(p)
        @assert N==length(h)
        @assert T==eltype(h)
        atomic_intensity(M,SVector{length(p)}(p),SVector{length(h)}(h))
    end

    @generated function atomic_intensity(M::Mty,p::K,h::K) where {T,N,K<:StaticArray,Mty<:AbstractArray{T,N}}
        Z=zero(T)
        O=one(T)
        quote
            @inbounds begin
                NP=IP=$Z
                Base.Cartesian.@nexprs $N i -> R_i = makerange(M,p,h,i)
                Base.Cartesian.@nloops $N i i->R_i begin
                    cp=$K(Base.Cartesian.@ntuple($N,k->i_k))
                    norm((p.-cp)./h) > 1 && continue;
                    IP+=Base.Cartesian.@nref($N,M,k->i_k)
                    NP+=$O
                end
                return (intensity=IP,numpoints=NP,ndims=N)
            end
        end
    end

using ProgressMeter

    function atomic_expmeanlog(M,p,h)
        N=length(p)
        T=eltype(p)
        @assert N==length(h)
        @assert T==eltype(h)
        atomic_expmeanlog(M,SVector{length(p)}(p),SVector{length(h)}(h))
    end

    @generated function atomic_expmeanlog(M::Mty,p::K,h::K) where {T,N,K<:StaticArray,Mty<:AbstractArray{T,N}}
        Z=zero(T)
        O=one(T)
        quote
            @inbounds begin
                NP2=NP=IP=IP2=$Z
                maxnum=maximum(M)
                Base.Cartesian.@nexprs $N i -> R_i = makerange(M,p,h,i)
                Base.Cartesian.@nloops $N i i->R_i begin
                    cp=$K(Base.Cartesian.@ntuple($N,k->i_k))
                    norm((p.-cp)./h)>1 && continue;
                    I=Base.Cartesian.@nref($N,M,k->i_k)
                    I<=0 && continue;
                    wh=(I/maxnum)^2   # pow=2 true est, pow=1 seems good too
                    X=log(I)
                    IP+=X*wh
                    IP2+=X^2*wh
                    NP+=wh
                    NP2+=wh^2
                end
                #@show IP,IP2,NP,NP2
                IP/=NP
                IP2/=NP
                XV=exp(IP)
                XD=XV*sqrt(max(IP2-IP^2,$Z))*sqrt(NP2)/NP
                return (XV,XD)
            end
        end
    end

"""
        atomic_ratio(Mnum,Mden,p,h) -> (pointest=(ratio,stddev),ndims=N)

Calculates the pixel pair-wise ratio on for density tensor `Mnum/Mden` over a ellipsoidal region caracterized by `p` region center, `h` semiaxis

# Arguments
- `Mnum` density tensor (e.g. matrix N*N)
- `Mden` density tensor (e.g. matrix N*N)
- `p` region center (e.g. 2D vector)
- `h` semiaxis lengths (e.g. 2D vector)

Note that both `p`,`h` should be of type Vector or StaticVector with the same length.
"""
    function atomic_ratio(Mnum,Mden,p,h)
        n,σn=atomic_expmeanlog(Mnum,p,h)
        d,σd=atomic_expmeanlog(Mden,p,h)
        (n/d, sqrt(σn^2+(σd*n/d)^2)/d)
    end

    function meanshift_kernel!(::Val{isadaptive} ,M::Mty,P::AbstractVector{K},h::AbstractVector{K},w::AbstractVector{T}, hmax::T, isotropy::T,maxit::T,rtol::T,smoothing::T) where {isadaptive,T,N,K<:StaticArray,Mty<:AbstractArray{T,N}}
        pr = Progress(length(P));
        update!(pr,0)
        jj = Threads.Atomic{Int}(0)
        a=@elapsed Threads.@threads for i in eachindex(P)
            cnt=0
            delta=one(T)
            bufcnt=0
            while (cnt<maxit) & (bufcnt < 10)
                if norm(h[i])>hmax && isadaptive
                    h[i]=h[i] ./ norm(h[i]) .* hmax
                end
                data=atomic_meanshift(M,P[i],h[i],isotropy,smoothing)
                if isadaptive
                    @inbounds P[i],h[i],w[i],delta=data
                else
                    @inbounds P[i],_,w[i],delta=data
                end
                cnt+=1
                bufcnt=ifelse(delta<rtol,bufcnt+1,0)
            end
            Threads.atomic_add!(jj, 1)
            Threads.threadid() == 1 && update!(pr, jj[])
        end
        isadaptive && Threads.@threads for i in eachindex(P)
                if norm(h[i])>hmax
                    h[i]=h[i] ./ norm(h[i]) .* hmax
                end
        end
        finish!(pr)
        a
    end

"""
        meanshift!(M, P, h, w, hmax; isotropy=½, maxit=∞, rtol=√ϵ, smoothing=1)
        meanshift_nonadaptive!(M, P, h, w; maxit=∞, rtol=√ϵ, smoothing=1)

Performs MeanShift on a swarm of particles over a given density matrix locating all modes and their scale. The non adaptive method disables updating of the `h` vectors, and thus depends on fewer parameters.

# Arguments
- `M` is the density tensor over which particles evolve.
- `P` is a julia vector of StaticVectors describing all positions.
- `h` is a julia vector of StaticVectors describing all standard deviations.
- `w` is a julia vector which will contain the new modes intensity estimates.
- `hmax` is a real number for the maximum norm allowed to the `h` vectors.

# Keyword Arguments
- `isotropy` is a scalar [0,1] set to 1 for isotropic kernel, < 1 anisotropy along coordinate axis is allowed.
- `maxit` is a integer for the maximum number of meanshift iterations.
- `rtol` is the absolute tolerance to declare a particle as converged.
- `smoothing` is a regularization term for density tensors with noise.

"""
    function meanshift!(M,P,h,w, hmax; isotropy = 0.5,maxit=Inf,rtol=1e-8,smoothing=1.0)
        meanshift_kernel!(Val(true),M,P,h,w, promote(hmax, isotropy,maxit,rtol,smoothing)...)
    end
    function meanshift_nonadaptive!(M,P,h,w; maxit=Inf,rtol=1e-8,smoothing=1.0)
        meanshift_kernel!(Val(false),M,P,h,w, promote(0, 0,maxit,rtol,smoothing)...)
    end
"""
    extract_region(M,p,h) -> M',p',h
Extracts a rectangular section with center "p" and axes "h" from a multidimensional array and returns the region with the new center coordinates.
"""
    function extract_region(M::Mty,p::Vty,h::Vty) where {T,N,Mty<:AbstractArray{T,N},Vty<:AbstractVector{T}}
        rg=[makerange(M,p,h,i) for i in 1:N]
        redM=collect(view(M,rg...))
        startind=minimum.(rg)
        redp=p.-startind.+1
        redM,redp,h
    end
"""
    extract_circ_region(M,p,h) -> M',p',h
Extracts a spherical section with center "p" and axes "h" from a multidimensional array and returns the region with the new center coordinates.
"""
    function extract_circ_region(M::Mty,p::Vty,h::Vty) where {T,N,Mty<:AbstractArray{T,N},Vty<:AbstractVector{T}}
        redM,redp,redh=extract_region(M,p,h)
        for i in CartesianIndices(redM)
            if norm((Tuple(i).-redp)./redh) > 1
                redM[i] =0.0
            end
        end
        redM,redp,redh
    end
end


text:
# Script for creating Figure 2 and 3 in scale manuscript

using IcanProj
using NetCDF
using DataFrames
using ProgressMeter
using Statistics
using CSV
using PyCall

function project_results(df::DataFrame, on::Symbol)

  map_values = fill(NaN, 1550, 1195)

  map_values[df[:ind_julia]] = df[on]

  return map_values

end

# Settings

figpath = joinpath(dirname(pathof(IcanProj)), "..", "plots", "forcings")

# Load data

df_all = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "forcings_summary.txt")) |> DataFrame

df_all.prec_mean = 365*8*(df_all.rainf_mean .+ df_all.snowf_mean)

# Plot averages

map_prec = project_results(df_all, :prec_mean)

map_wind = project_results(df_all, :wind_mean)

map_tair = project_results(df_all, :tair_mean)

file = joinpath(figpath, "forcings_mean.pdf")

py"""
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np

fig, ax = plt.subplots(1, 3)

fig.set_size_inches(12, 8)

# Precipitation

im = ax[0].imshow($(map_prec), cmap = 'jet', vmin=-0.1, vmax=4000)

ax[0].grid(linestyle='dotted')
ax[0].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[0].annotate("(A)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[0])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Precipitation (mm/year)")

# Wind speed

im = ax[1].imshow($(map_wind), cmap = 'jet', vmin=-0.1, vmax=9)

ax[1].grid(linestyle='dotted')
ax[1].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[1].annotate("(B)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[1])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Wind speed (m/s)")

# Air temperature

im = ax[2].imshow($(map_tair), cmap = 'jet', vmin=-9, vmax=9)

ax[2].grid(linestyle='dotted')
ax[2].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[2].annotate("(C)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[2])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Air temperature (°C)")

# Save plot

fig.savefig($(file), dpi = 600, bbox_inches='tight')

plt.close()
"""

# Plot standard deviations

df_tmp = by(df_all, :ind_10km, d -> DataFrame(tair_std = std(d[:tair_mean]), wind_std = std(d[:wind_mean]), prec_std = std(d[:prec_mean])))

df_10km = join(df_all, df_tmp, on = :ind_10km)

df_tmp = by(df_all, :ind_50km, d -> DataFrame(tair_std = std(d[:tair_mean]), wind_std = std(d[:wind_mean]), prec_std = std(d[:prec_mean])))

df_50km = join(df_all, df_tmp, on = :ind_50km)

prec_10km = project_results(df_10km, :prec_std)

wind_10km = project_results(df_10km, :wind_std)

tair_10km = project_results(df_10km, :tair_std)

prec_50km = project_results(df_50km, :prec_std)

wind_50km = project_results(df_50km, :wind_std)

tair_50km = project_results(df_50km, :tair_std)

prec_max = maximum(df_50km.prec_std[.!isnan.(df_50km.prec_std)])

wind_max = maximum(df_50km.wind_std[.!isnan.(df_50km.wind_std)])

tair_max = maximum(df_50km.tair_std[.!isnan.(df_50km.tair_std)])

file = joinpath(figpath, "forcings_std.pdf")

py"""
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np

fig, ax = plt.subplots(2, 3)

fig.set_size_inches(18, 12)

# Precipitation - 10 km

im = ax[0,0].imshow($(prec_10km), cmap = 'jet', vmin=-0.1, vmax=$(prec_max))

ax[0,0].grid(linestyle='dotted')
ax[0,0].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[0,0].annotate("(A)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[0,0])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
cb = fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Precipitation (mm/year)")

cb.remove()

# Precipitation - 50 km

im = ax[1,0].imshow($(prec_50km), cmap = 'jet', vmin=-0.1, vmax=$(prec_max))

ax[1,0].grid(linestyle='dotted')
ax[1,0].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[1,0].annotate("(D)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[1,0])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Precipitation (mm/year)")

# Wind speed - 10 km

im = ax[0,1].imshow($(wind_10km), cmap = 'jet', vmin=-0.1, vmax=$(wind_max))

ax[0,1].grid(linestyle='dotted')
ax[0,1].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[0,1].annotate("(B)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[0,1])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
cb = fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Wind speed (m/s)")
cb.remove()

# Wind speed - 50 km

im = ax[1,1].imshow($(wind_50km), cmap = 'jet', vmin=-0.1, vmax=$(wind_max))

ax[1,1].grid(linestyle='dotted')
ax[1,1].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[1,1].annotate("(E)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[1,1])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Wind speed (m/s)")

# Air temperature - 10 km

im = ax[0,2].imshow($(tair_10km), cmap = 'jet', vmin=-0.1, vmax=$(tair_max))

ax[0,2].grid(linestyle='dotted')
ax[0,2].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[0,2].annotate("(C)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[0,2])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
cb = fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Air temperature (°C)")
cb.remove()

# Air temperature - 50 km

im = ax[1,2].imshow($(tair_50km), cmap = 'jet', vmin=-0.1, vmax=$(tair_max))

ax[1,2].grid(linestyle='dotted')
ax[1,2].tick_params(axis='both', left=False, top=False, right=False, bottom=False, labelleft=False, labeltop=False, labelright=False, labelbottom=False)
ax[1,2].annotate("(F)", xy=(0.05,0.90), xycoords="axes fraction", fontsize = 14)

divider = make_axes_locatable(ax[1,2])
cax = divider.new_vertical(size="5%", pad=0.05, pack_start=True)
fig.add_axes(cax)
fig.colorbar(im, cax=cax, orientation="horizontal", label=u"Air temperature (°C)")

# Save plot

fig.savefig($(file), dpi = 200, bbox_inches='tight')

plt.close()
"""


text:
using Pkg
Pkg.activate(".")

using Images, CausticsEngineering

image = Images.load("./input/logos_patente512x512.jpg") # Check current working directory with pwd()
engineer_caustics(image);

max_stars_repo_name: timtomch/causticsEngineering

max_issues_repo_name: timtomch/causticsEngineering

max_forks_repo_name: timtomch/causticsEngineering