text:
module SubsetSelection
using Compat, LinearAlgebra, Statistics, Random

import Compat.String

export LossFunction, Regression, Classification, OLS, L1SVR, L2SVR, LogReg, L1SVM, L2SVM
export Sparsity, Constraint, Penalty
export SparseEstimator, subsetSelection

include("types.jl")
include("stepsize.jl")
include("recover_primal.jl")

# Type to hold preallocated memory
struct Cache
  g::Vector{Float64}
  ax::Vector{Float64}
  sortperm::Vector{Int}

  function Cache(n::Int, p::Int)
    new(
      Vector{Float64}(undef,n),
      Vector{Float64}(undef,p),
      Vector{Int}(undef,p),
    )
  end
end

mutable struct bestSupport
  indices::Vector{Int}
  nindices::Int
  w::Vector{Float64}
  ub::Float64 #UB corresponds to the upper bound on the objective value of the relaxation. obj ≧ UB, but sometimes obj > UB when relaxation not tight
  obj::Float64
  lb::Float64
end
##############################################
##DUAL SUB-GRADIENT ALGORITHM
##############################################
""" Function to compute a sparse regressor/classifier. Solve an optimization problem of the form
        min_s max_α f(α, s)
by gradient ascent on α:        α_{t+1} ← α_t + δ ∇_α f(α_t, s_t)
and partial minimization on s:  s_{t+1} ← argmin_s f(α_t, s)

INPUTS
- ℓ           Loss function used
- Card        Model to enforce sparsity (constraint or penalty)
- Y (n×1)     Vector of outputs. For classification, use ±1 labels
- X (n×p)     Array of inputs.
- indInit     (optional) Initial subset of features s
- αInit       (optional) Initial dual variable α
- γ           (optional) ℓ2 regularization penalty
- intercept   (optional) Boolean. If true, an intercept term is computed as well

- stepRule    (optional) Stepsize rule to be used in the subgradient algorithm. Default Poliak's
- maxIter     (optional) Total number of iterations. Default 500
- numberRestarts     (optional) Number of restarting. Default 1
Stopping criterion
- noImprov_threshold      (optional) Algorithm stops after noImprov_threshold iterations without improvement
- dGap (optional) Algorithm stops if the duality + integrality gap drops below dGap

OUTPUT
- SparseEstimator """
function subsetSelection(ℓ::LossFunction, Card::Sparsity, Y, X;
    indInit = ind_init(Card, size(X,2)), αInit=alpha_init(ℓ, Y),
    γ = 1/sqrt(size(X,1)),  intercept = false,
    stepRule::StepSizeRule = poliakRule(1),
    maxIter = 200, numberRestarts = 1,
    noImprov_threshold = maxIter, dGap = 1e-4)

  n,p = size(X)
  cache = Cache(n, p)

  #Add sanity checks
  if size(Y,1) != n
    throw(DimensionMismatch("X and Y must have the same number of rows"))
  end
  if isa(ℓ, SubsetSelection.Classification)
    levels = sort(unique(Y))
    if length(levels) != 2
      throw(ArgumentError("subsetSelection only supports two-class classification"))
    elseif (levels[1] != -1) || (levels[2] != 1)
      throw(ArgumentError("Class labels must be ±1's"))
    end
  end

  #Initialize variables
  indices = indInit #Current support
  n_indices = length(indices) #Number of non-zero elements

  n_indices_max = max_index_size(Card, p) #Max size of the support
  resize!(indices, n_indices_max+1)

  α = αInit[:]  #Dual variable α
  a = αInit[:]  #Past average of α

  w = recover_primal(ℓ, Y, X[:,indices[1:n_indices]], γ)  #Current regressor / primal variable
  resize!(w, n_indices_max+1)

  obj_init = value_primal(ℓ, Y, X[:,indices[1:n_indices]], w[1:n_indices], γ, cache) #Current primal objective

  #Initialize memory
  subsetPool = Set() #Memory: set of visited supports
  push!(subsetPool, indices[1:n_indices])
  best_s = bestSupport(indices[:], n_indices, w[:], obj_init, obj_init, -Inf) #Memory: best solution so far
  consecutive_noimprov = 0

  ##Dual Sub-gradient Algorithm
  # @showprogress 2 "Feature selection in progress... "
  effectiveEpoch = 0
  for restart in 1:numberRestarts
    α[:] = a[:] #Set initial point to the average
    effectiveEpoch += 1

    for iter in 2:round(Int, maxIter/numberRestarts)
      effectiveEpoch += 1
      #Minimization w.r.t. s
      n_indices, tie = partial_min!(indices, Card, X, α, γ, cache)

      if tie
        #If one tie in the k best features, UB for the relaxation is obtained by including all (k+1) features
        #Warning: in this case, we generate an infeasible solution
        #TO DO: take into account multiple ties
        w[1:(n_indices+1)] = recover_primal(ℓ, Y, X[:,indices[1:(n_indices+1)]], γ)
        upper_bound = value_primal(ℓ, Y, X[:,indices[1:(n_indices+1)]], w[1:(n_indices+1)], γ, cache)
        best_s.ub = min(best_s.ub, upper_bound)
      else
        if indices[1:n_indices] ∉ subsetPool #if support has never been seen
          push!(subsetPool, indices[1:n_indices])

          w[1:n_indices] = recover_primal(ℓ, Y, X[:,indices[1:n_indices]], γ)
          upper_bound = value_primal(ℓ, Y, X[:,indices[1:n_indices]], w[1:n_indices], γ, cache)
          print("upper bound", upper_bound)

          if upper_bound < best_s.obj
              consecutive_noimprov = 0

              best_s.indices[1:n_indices] = indices[1:n_indices]
              best_s.nindices = n_indices
              best_s.w[1:n_indices] = w[1:n_indices]

              best_s.obj = upper_bound
              best_s.ub = min(best_s.ub, upper_bound)
          else
              consecutive_noimprov += 1
          end
        else
          consecutive_noimprov += 1
        end
      end

      best_s.lb = max(best_s.lb, value_dual(ℓ, Y, X, α, indices, n_indices, γ))

      #Gradient ascent on α
      ∇ = grad_dual(ℓ, Y, X, α, indices, n_indices, γ, cache)

      #Stop if small gradient
      if norm(∇, 1) <= 1e-14
        maxIter = iter
        break
      end

      #For numeric stability
      if isa(ℓ, SubsetSelection.Classification) && norm(∇, 1) == Inf
        pInfIndex = findall(∇ .== Inf)
        nInfIndex = findall(∇ .== -Inf)

        α[pInfIndex] .= -Y[pInfIndex]*1e-14
        α[nInfIndex] .= -Y[nInfIndex]*(1-1e-14)

        ∇[pInfIndex] .= 0.
        ∇[nInfIndex] .= 0.
      end

      δ =  computeStepSize(stepRule, iter, ∇, best_s.ub, best_s.lb)
      α .+= δ*∇

      proj_dual!(ℓ, Y, α)
      proj_intercept!(intercept, α)

      #Update average a
      @__dot__ a = (iter - 1) / iter * a + 1 / iter * α

      #Duality gap rule
      if (best_s.ub - best_s.lb) / abs(best_s.lb) < dGap
        break
      end

      #No improvement rule
      if consecutive_noimprov >= noImprov_threshold
        break
      end
    end

    slashStep!(stepRule, 2.)
  end

  ##Compute sparse estimator
  #Subset of relevant features
  n_indices, = partial_min!(indices, Card, X, a, γ, cache)
  #Regressor
  w[1:n_indices] = recover_primal(ℓ, Y, X[:,indices[1:n_indices]], γ)

  objvalue = value_primal(ℓ, Y, X[:,indices[1:n_indices]], w[1:n_indices], γ, cache)
  if objvalue > best_s.obj #If last solution worse than the best found
      n_indices = best_s.nindices
      indices[1:n_indices] = best_s.indices[1:n_indices]
      w[1:n_indices] = best_s.w[1:n_indices]
  end

  #Bias
  b = compute_bias(ℓ, Y, X, a, indices, n_indices, γ, intercept, cache)

  #Resize final indices vector to only have relevant entries
  resize!(indices, n_indices)

  return SparseEstimator(ℓ, Card, indices, w[1:n_indices], a, b, effectiveEpoch)
end

# Helper to check if indices and indices_old are the same
# Equivalent to `indices[1:n_indices] == indices_old[1:n_indices]` but without
# allocating temp arrays
function indices_same(indices, indices_old, n_indices)
  for j = 1:n_indices
    if indices[j] != indices_old[j]
      return false
    end
  end
  true
end

##############################################
##AUXILIARY FUNCTIONS
##############################################

##Default Initialization of prinal variables s depending on the model used
function ind_init(Card::Constraint, p::Integer)
  indices0 = findall(x-> x<Card.k/p, rand(p))
  while !(Card.k >= length(indices0) >= 1)
    indices0 = findall(x-> x<Card.k/p, rand(p))
  end
  return indices0
end
function ind_init(Card::Penalty, p::Integer)
  indices0 = findall(x-> x<1/2, rand(p))
  while !(length(indices0)>=1)
    indices0 = findall(x-> x<1/2, rand(p))
  end
  return indices0
end

##Default Initialization of dual variables α depending on the model used
function alpha_init(ℓ::OLS, Y)
  return -Y
end
function alpha_init(ℓ::L1SVR, Y)
  return max.(-1,min.(1, -Y))
end
function alpha_init(ℓ::L2SVR, Y)
  return -Y
end
function alpha_init(ℓ::Classification, Y)
  return -Y./2
end

##Dual objective function value for a given dual variable α
function value_dual(ℓ::LossFunction, Y, X, α, indices, n_indices, γ)
  v = - sum([fenchel(ℓ, Y[i], α[i]) for i in 1:size(X, 1)])
  for j in 1:n_indices
    v -= γ/2*(dot(X[:, indices[j]], α)^2)
  end
  return v
end

##Gradient of f w.r.t. the dual variable α
function grad_dual(ℓ::LossFunction, Y, X, α, indices, n_indices, γ, cache::Cache)
  g = cache.g
  for i in 1:size(X, 1)
    g[i] = -grad_fenchel(ℓ, Y[i], α[i])
  end
  for j in 1:n_indices
    x = @view(X[:, indices[j]])
    # @__dot__ g -= γ * dot(x, α) * x
    g -= γ*dot(x, α)*x
  end
  g
end

##Projection of α on the feasible set of the Fenchel conjugate
function proj_dual!(ℓ::OLS, Y, α)
end
function proj_dual!(ℓ::L1SVR, Y, α)
  α[:] = max.(-1,min.(1, α))
end
function proj_dual!(ℓ::L2SVR, Y, α)
end
function proj_dual!(ℓ::LogReg, Y, α)
  α[:] = Y.*max.(-1,min.(0, Y.*α))
end
function proj_dual!(ℓ::L1SVM, Y, α)
  α[:] = Y.*max.(-1,min.(0, Y.*α))
end
function proj_dual!(ℓ::L2SVM, Y, α)
  α[:] = Y.*min.(0, Y.*α)
end

##Projection of α on e^T α = 0 (if intercept)
function proj_intercept!(intercept::Bool, α)
  if intercept
    α[:] .-= mean(α)
  end
  # return α
end

##Minimization w.r.t. s
function partial_min!(indices, Card::Constraint, X, α, γ, cache::Cache)
  ax = cache.ax
  perm = cache.sortperm
  p = size(X,2)
  n_indices = max_index_size(Card, p)

  # compute α'*X into pre-allocated scratch space
  mul!(ax, X', α)
  # take the k largest (absolute) values of ax
  map!(abs, ax, ax)

  sortperm!(perm, ax, rev=true)

  indices[1:n_indices] = perm[1:n_indices]

  sort!(@view(indices[1:n_indices]))

  tie = false
  if n_indices < p
    if abs(ax[perm[n_indices+1]] - ax[perm[n_indices]])/ax[perm[n_indices]] <= 1e-6
      tie = true
      indices[n_indices+1] = perm[n_indices+1]
    end
  end
  # Return the updated size of indices
  return n_indices, tie
end

function partial_min!(indices, Card::Penalty, X, α, γ, cache::Cache)
  ax = cache.ax

  # compute (α'*X).^2 into pre-allocated scratch space
  mul!(ax, X', α)
  map!(abs2, ax, ax)

  # find indices with `λ - γ / 2 * (a'X_j)^2 < 0`
  n_indices = 0
  for j = 1:size(X,2)
    if Card.λ - γ / 2 * ax[j] < 0
      n_indices += 1
      indices[n_indices] = j
    end
  end
  tie = false #TO DO
  return n_indices, tie
end

##Bias term
function compute_bias(ℓ::LossFunction, Y, X, α, indices, n_indices, γ,
                      intercept::Bool, cache::Cache)
  if intercept
    g = grad_dual(ℓ, Y, X, α, indices, n_indices, γ, cache)
    return (minimum(g[α != 0.]) + maximum(g[α != 0.]))/2
  else
    return 0.
  end
end

##Primal objective function value for a given primal variable w
function value_primal(ℓ::LossFunction, Y, X, w, γ, cache::Cache)
  g = cache.g
  for i in 1:size(X, 1)
    g[i] = loss(ℓ, Y[i], dot(X[i,:], w))
  end
  return sum(g) + dot(w,w)/2/γ
  # return sum(loss(ℓ, Y, X*w)) + dot(w,w)/2/γ
  # v = sum([loss(ℓ, Y[i], dot(X[i,:], w)) for i in 1:size(X, 1)])
  # v += dot(w,w)/2/γ
  # return v
end

end #module


text:
module DERs
using HCEstimator
export add_ders!

struct Config
    active_dispatchable::Bool
    reactive_dispatchable::Bool
    power::Float64
    energy::Float64
    α::Float64
    β::Float64
end

function configure_limits!(MW::Vector{Float64}, MVAr::Vector{Float64}, configs::Config)::Tuple{Vector{Float64},Vector{Float64}}
    if !configs.active_dispatchable
        MW = [0.0, 0.0]
    end

    if !configs.reactive_dispatchable
        MVAr = [0.0, 0.0]
    end

    return MW, MVAr
end

function dg(configs::Config, bus)
    power = configs.power
    MW = [0.0, power]
    MVAr = [-power, power]

    MW, MVAr = configure_limits!(MW, MVAr, configs)

    if isapprox(configs.α, 1.0; atol=1e-3)
        oper = [0.0]
    else
        oper = [0.0 , 1]
    end

    return DistSystem.DER(
        bus,                  # Bus
        power,                # Sᴰᴱᴿ (MVA)
        configs.energy,                # Eᴰᴱᴿ (MWh)
        configs.α,                # αᴰᴱᴿ
        configs.β,                # βᴰᴱᴿ
        MW,         # [Pᴰᴱᴿ_low, Pᴰᴱᴿ_upper]  (MW)
        MVAr,       # [Qᴰᴱᴿ_low, Qᴰᴱᴿ_upper]  (MVAr)
        oper,          # Possible DER's Operation Scenarios ≠ μᴰᴱᴿ
        [0.0] # Costs (Not used)
    )
end

function ess(configs::Config, bus)
    power = configs.power
    MW = [-power, power]
    MVAr = [-power, power]

    MW, MVAr = configure_limits!(MW, MVAr, configs)

    if isapprox(configs.α, 1.0; atol=1e-3)
        oper = [0.0]
    else
        oper = [-1.0, 0.0, 1.0]
    end

    return DistSystem.DER(
        bus,                  # Bus
        power,                # Sᴰᴱᴿ (MVA)
        configs.energy,                # Eᴰᴱᴿ (MWh)
        configs.α,                # αᴰᴱᴿ
        configs.β,                # βᴰᴱᴿ
        MW,         # [Pᴰᴱᴿ_low, Pᴰᴱᴿ_upper]  (MW)
        MVAr,       # [Qᴰᴱᴿ_low, Qᴰᴱᴿ_upper]  (MVAr)
        oper,          # Possible DER's Operation Scenarios ≠ μᴰᴱᴿ
        [0.0] # Costs (Not used)
    )
end

function ev_charger(configs::Config, bus)
    power = configs.power
    MW = [-power, 0.0]
    MVAr = [-power, power]
    MW, MVAr = configure_limits!(MW, MVAr, configs)
    return DistSystem.DER(
        bus,                  # Bus
        power,                # Sᴰᴱᴿ (MVA)
        configs.energy,                # Eᴰᴱᴿ (MWh)
        configs.α,                # αᴰᴱᴿ
        configs.β,                # βᴰᴱᴿ
        MW,         # [Pᴰᴱᴿ_low, Pᴰᴱᴿ_upper]  (MW)
        MVAr,       # [Qᴰᴱᴿ_low, Qᴰᴱᴿ_upper]  (MVAr)
        [-1.0, 0.0],          # Possible DER's Operation Scenarios ≠ μᴰᴱᴿ
        [0.0] # Costs (Not used)
    )
end

end

ext: jl


text:
module NormalizeType
export normalize_typevars, split_typevar

using WhereTraits.Utils
using ExprParsers

const VectorLike = Union{Vector, Tuple, Core.SimpleVector}
const True = Val{true}
const False = Val{false}
const TrueFalse = Union{True, False}
Base.:(!)(::True) = False()
Base.:(!)(::False) = True()

"""
when normalizing a type signature, we have need to distinguish TypeVar according to three different context
- whether we need to track the concrete normalization-renaming (normalize_track)
- whether we just need to normalize (normalize_notrack)
- or whether a typevar is not to be normalized (neither of the other)
"""
Base.@kwdef struct TypeVarContext
  normalize_notrack::Vector{TypeVar} = TypeVar[]
  normalize_track::Vector{TypeVar} = TypeVar[]
end

"""
traverses given types, renaming all typeparameters to normalized name

returns 3-tuple
- type with substituted typevars
- list of function level typevars (includable into where-statement)
- dict mapping new typevar to old once
"""
function normalize_typevars(types::VectorLike, typevars_track_normalization)
  tuple_parents_only = True()
  countfrom_typevars! = Base.Iterators.Stateful(Base.Iterators.countfrom(1))
  typevar_old_to_new! = Dict{TypeVar, Symbol}()
  typevar_context = TypeVarContext(normalize_track = typevars_track_normalization)
  typeexprs, typevars = _normalize_typevars(types, typevar_context, typevar_old_to_new!, countfrom_typevars!, tuple_parents_only)
  typevar_new_to_old = Dict(v => k.name for (k, v) in typevar_old_to_new!)
  typeexprs, typevars, typevar_new_to_old
end

function _normalize_typevars(types::VectorLike, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only::TrueFalse)
  typeexprs = []
  typevars = []
  for p in types
    _typeexpr, _typevars = _normalize_typevars(p, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only)
    push!(typeexprs, _typeexpr)
    append!(typevars, _typevars)
  end
  typeexprs, typevars
end

function _simplify_type(varargexpr::Expr)
  parsed = parse_expr(EP.Type(), varargexpr)

  isdroppable(other) = false
  isdroppable(sym::Symbol) = sym in parsed.wheres

  droppable_curlies = []
  for sym in reverse(parsed.curlies)
    isdroppable(sym) || break
    push!(droppable_curlies, sym)
  end
  parsed.curlies = filter(x -> x ∉ droppable_curlies, parsed.curlies)
  parsed.wheres = filter(x -> x ∉ droppable_curlies, parsed.wheres)
  to_expr(parsed)
end

# UnionAlls and abstracttypes are treated differently in that they will create function level typevariables
function _normalize_typevars(type::Type, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only::True)
  typebase, extra_typevars = split_typevar(type)
  name = _name(typebase)

  if name === :Vararg
    # Vararg behave different than all other types when used as TypeVar for Tuple, namely
    # Tuple{Vararg{T} where T} != (Tuple{Vararg{T}} where T)
    # which is the same as if tuple_parents_only == False
    varargexpr, typevars = _normalize_typevars(type, typevar_context, typevar_old_to_new!, countfrom!, False())

    # we need to simplify the type because of a bug in julia
    # https://discourse.julialang.org/t/inconsistent-behaviour-of-vararg-unionall/31644
    _simplify_type(varargexpr), typevars

  elseif isabstracttype(type) || type isa UnionAll || name === :Union
    # each abstract or unionall type or Union type can have different concrete types
    # hence it makes sense to give it a typevariable like
    # Tuple{Vector} <: Tuple{V} where {V <: Vector}
    # more concretely, the typeparameter representing this type must
    # come after the typeparameters of the type
    # Tuple{Vector{T} where T} <: Tuple{V} where {T, V <: Vector{T}}

    tuple_parents_only = name === :Tuple ? True() : False()
    # extra_typevars need to get normalized too
    typevar_context′ = TypeVarContext(
      normalize_track = typevar_context.normalize_track,
      normalize_notrack = [typevar_context.normalize_notrack; extra_typevars],
    )
    typeexpr, typevars = _normalize_typevars(typebase, typevar_context′, typevar_old_to_new!, countfrom!, tuple_parents_only)

    # add type variable representing the entire type
    i = next!(countfrom!)
    new = normalized_typevar_by_position(i, typeexpr)
    # we replace the given type with the new typevariablename
    # and add the constraint to the function level typevars
    new.name, [typevars; new]

  elseif isempty(type.parameters)  # concrete DataType without typeparameters
    type, []

  else
    # recurse into typeparameters of concrete DataType
    tuple_parents_only = name === :Tuple ? True() : False()
    typeexprs, typevars = _normalize_typevars(typebase.parameters, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only)
    Expr(:curly, name, typeexprs...), typevars
  end
end
_name(type::Type) = type.name.name
_name(type::Union) = :Union

function _normalize_typevars(tv::TypeVar, typevar_context, typevar_old_to_new!, countfrom!::Base.Iterators.Stateful, tuple_parents_only::TrueFalse)
  if tv in keys(typevar_old_to_new!)
    # we may have complex interactions like diagonal types and other reuses of typevariables
    # hence we make sure we replace a variable only once
    new_name = typevar_old_to_new![tv]
    # we replace the typevar with new_name, but there is no new typevar
    new_name, []
  elseif tv in typevar_context.normalize_track
    i = next!(countfrom!)
    lowerbound, lb_vars = _normalize_typevars(tv.lb, typevar_context, typevar_old_to_new!, countfrom!, False())
    upperbound, ub_vars = _normalize_typevars(tv.ub, typevar_context, typevar_old_to_new!, countfrom!, False())
    new = normalized_typevar_by_position(i, lowerbound, upperbound)
    typevar_old_to_new![tv] = new.name
    new.name, [lb_vars; ub_vars; new]
  elseif tv in typevar_context.normalize_notrack
    # notrack, means this is implicit
    # hence all these typevars are independent and can just be replaced
    i = next!(countfrom!)
    lowerbound, lb_vars = _normalize_typevars(tv.lb, typevar_context, typevar_old_to_new!, countfrom!, False())
    upperbound, ub_vars = _normalize_typevars(tv.ub, typevar_context, typevar_old_to_new!, countfrom!, False())
    new = normalized_typevar_by_position(i, lowerbound, upperbound)
    new.name, [lb_vars; ub_vars; new]
  else  # this is a nested UnionAll somewhere
    tv.name, []
  end
end

# we need to solve ambiguities with type::Union (vs type::Type) and tuple_parents_only::TrueFalse (vs ::True)
function _normalize_typevars(type::Union, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only::True)
  _normalize_typevars_unions(type, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only)
end
function _normalize_typevars(type::Union, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only::False)
  _normalize_typevars_unions(type, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only)
end
function _normalize_typevars_unions(type::Union, typevar_context, typevar_old_to_new!, countfrom!, tuple_parents_only::TrueFalse)
  typeexprs, typevars = _normalize_typevars(_collect_union_types(type), typevar_context, typevar_old_to_new!, countfrom!, False())
  Expr(:curly, :Union, typeexprs...), typevars
end
_collect_union_types(any) = [any]
_collect_union_types(type::Union) = [type.a; _collect_union_types(type.b)]

function _normalize_typevars(type::DataType, typevar_context, typevar_old_to_new!, countfrom!, ::False)
  # if tuple_parents_only == False, `type` cannot be a Vararg, as Varargs are only allowed in Tuples
  if isempty(type.parameters)  # concrete DataType without typeparameters
    type, []
  else
    # recurse into typeparameters of concrete DataType
    typeexprs, typevars = _normalize_typevars(type.parameters, typevar_context, typevar_old_to_new!, countfrom!, False())
    Expr(:curly, type.name.name, typeexprs...), typevars
  end
end

function _normalize_typevars(type::UnionAll, typevar_context, typevar_old_to_new!, countfrom!, ::False)
  # if there is not everywhere tuple types above, this UnionAll stays as it is, only renaming the inner Type
  typebase, typevars_unionall = split_typevar(type)
  typevars_unionall_names = [tv.name for tv in typevars_unionall]
  typevars_track_unshadowed = [tv for tv ∈ typevar_context.normalize_track if tv.name ∉ typevars_unionall_names]
  typevars_notrack_unshadowed = [tv for tv ∈ typevar_context.normalize_notrack if tv.name ∉ typevars_unionall_names]
  typevar_context′ = TypeVarContext(
    normalize_track = typevars_track_unshadowed,
    normalize_notrack = typevars_notrack_unshadowed,
  )
  typebase_expr_new, typevars_new = _normalize_typevars(typebase, typevar_context′, typevar_old_to_new!, countfrom!, False())
  typevars_unionall_new, typevars_new′ = _normalize_unionall(typevars_unionall, typevar_context′, typevar_old_to_new!, countfrom!, False())
  _expr_rewrap_typevars(typebase_expr_new, typevars_unionall_new), [typevars_new; typevars_new′]
end

function _normalize_unionall(types::VectorLike, typevar_context, typevar_old_to_new!, countfrom!::Base.Iterators.Stateful, tuple_parents_only::TrueFalse)
  typeexprs = []
  typevars = []
  for p in types
    _typeexpr, _typevars = _normalize_unionall(p, typevar_context, typevar_old_to_new!, countfrom!, False())
    push!(typeexprs, _typeexpr)
    append!(typevars, _typevars)
  end
  typeexprs, typevars
end
function _normalize_unionall(tv::TypeVar, typevar_context, typevar_old_to_new!, countfrom!::Base.Iterators.Stateful, tuple_parents_only::TrueFalse)
  lowerbound, lb_vars = _normalize_typevars(tv.lb, typevar_context, typevar_old_to_new!, countfrom!, False())
  upperbound, ub_vars = _normalize_typevars(tv.ub, typevar_context, typevar_old_to_new!, countfrom!, False())
  to_expr(EP.TypeRange_Parsed(lowerbound, tv.name, upperbound)), [lb_vars; ub_vars]
end

# anything else, like Int or Symbol TypeVars
_normalize_typevars(any, typevars_old, typevar_old_to_new!, countfrom!, ::TrueFalse) = any, []
_normalize_typevars(bounds::Type{Union{}}, _, _, _, ::True) = bounds, []
_normalize_typevars(bounds::Type{Union{}}, _, _, _, ::False) = bounds, []

normalized_typevar_by_position(position::Int) = EP.TypeRange_Parsed(Union{}, Symbol("T", position), Any)
normalized_typevar_by_position(position::Int, ub) = EP.TypeRange_Parsed(Union{}, Symbol("T", position), ub)
normalized_typevar_by_position(position::Int, lb, ub) = EP.TypeRange_Parsed(lb, Symbol("T", position), ub)

split_typevar(base) = base, TypeVar[]
function split_typevar(t::UnionAll)
  base, typevars = split_typevar(t.body)
  base, [t.var; typevars...]
end

function _expr_rewrap_typevars(typeexpr::Expr, typevars)
  Expr(:where, typeexpr, typevars...)
end
function _expr_rewrap_typevars(typeexpr::Expr, typevars::Vector{TypeVar})
  _expr_rewrap_typevars(typeexpr, to_expr.(typevars))
end

end  # module


text:
immutable NoncentralChisq <: ContinuousUnivariateDistribution
    df::Float64
    ncp::Float64
    function NoncentralChisq(d::Real, nc::Real)
    	d >= zero(d) && nc >= zero(nc) || error("df and ncp must be non-negative")
    	new(float64(d), float64(nc))
    end
end

mean(d::NoncentralChisq) = d.df + d.ncp
var(d::NoncentralChisq) = 2.0*(d.df + 2.0*d.ncp)
skewness(d::NoncentralChisq) = 2.0*sqrt2*(d.df + 3.0*d.ncp)/sqrt(d.df + 2.0*d.ncp)^3
kurtosis(d::NoncentralChisq) = 12.0*(d.df + 4.0*d.ncp)/(d.df + 2.0*d.ncp)^2

function mgf(d::NoncentralChisq, t::Real)
    k = d.df
    exp(d.ncp * t/(1.0 - 2.0 * t))*(1.0 - 2.0 * t)^(-k / 2.0)
end

cf(d::NoncentralChisq, t::Real) = exp(im * d.ncp * t/(1.0 - 2.0 * im * t))*(1.0 - 2.0 * im * t)^(-d.df / 2.0)

@_jl_dist_2p NoncentralChisq nchisq

@continuous_distr_support NoncentralChisq 0.0 Inf


text:
# Outputs a compressed integer array
function pack{T<:Integer}(x::Vector{T})
    # Length of x
    n = length(x)
    # Number of 128 integer chunks in compressed array
    nchunks = cld(n, 128)
    # Stores chunked compressed data
    data = Vector{Vector{m128i}}(0)
    # Bit width of compressed integer chunks
    b = Vector{UInt8}(nchunks)
    # Offset of each chunk
    offset = Vector{UInt32}(nchunks)
    # A temporary, guaranteed contiguous array to store chunks of x
    xtemp = Vector{UInt32}(128)
    xptr = pointer(xtemp)
    # Maximum bit width of chunk
    bmax = Vector{UInt32}(1)
    # Minimum bit width of chunk
    bmin = Vector{UInt32}(1)
    bmaxptr = pointer(bmax)
    bminptr = pointer(bmin)
    for chunk = 1:nchunks
        if chunk != nchunks
            ntrail = 0
        else
            # Pad the last chunk with zeros so it fits 128 integers
            ntrail = nchunks*128 - n
            fill!(xtemp, 0)
        end
        xstart = (chunk - 1)*128
        for i = 1:128 - ntrail
            @inbounds xtemp[i] = x[xstart + i]
        end
        # Evaluate and store offset
        offsetval = minimum(xtemp)
        offset[chunk] = offsetval
        # Evaluate and store bit width
        ccall(simdmaxmin, Void, (Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), xptr, bminptr, bmaxptr)
        bval = ccall(simdbits, UInt32, (UInt32,), bmax[1] - bmin[1])
        b[chunk] = bval
        # A compressed data chunk
        dataentry = Vector{m128i}(bval)
        dataentryptr = pointer(dataentry)
        # Compress
        ccall(simdpackFOR, Void, (UInt32, Ptr{UInt32}, Ptr{m128i}, UInt32), offsetval, xptr, dataentryptr, bval)
        # Add compressed chunk to data
        push!(data, dataentry)
    end
    return TinyIntVec(data, b, n, offset)
end

# Outputs an uncompressed array
function unpack(tinyx::TinyIntVec)
    # Outputs an uncompressed array
    x = zeros(UInt32, length(tinyx.b)*128)
    nchunks = length(tinyx.b)
    for chunk = 1:nchunks
        xptr = pointer(x, (chunk-1)*128 + 1)
        dataentryptr = pointer(tinyx.data[chunk])
        ccall(simdunpackFOR, Void, (UInt32, Ptr{m128i}, Ptr{UInt32}, UInt32), tinyx.offset[chunk], dataentryptr, xptr, tinyx.b[chunk])
    end
    # Remove padding
    resize!(x, tinyx.n)
    return x
end

# Overwrite x with an uncompressed array
function unpack!{T<:Integer}(tinyx::TinyIntVec, x::Vector{T})
    # Guaranteed contiguous temporary array
    xtemp = Vector{UInt32}(128)
    xptr = pointer(xtemp)
    nchunks = length(tinyx.b)
    for chunk = 1:nchunks
        dataentryptr = pointer(tinyx.data[chunk])
        ccall(simdunpackFOR, Void, (UInt32, Ptr{m128i}, Ptr{UInt32}, UInt32), tinyx.offset[chunk], dataentryptr, xptr, tinyx.b[chunk])
        ntrail = ifelse(chunk != nchunks, 0, nchunks*128 - tinyx.n)
        xstart = (chunk - 1)*128
        for i = 1:128 - ntrail
            @inbounds x[xstart + i] = xtemp[i]
        end
    end
end

# Outputs an uncompressed chunk
function unpack(tinyx::TinyIntVec, chunk::Integer)
    # Outputs an uncompressed array
    x = Vector{UInt32}(128)
    xptr = pointer(x, (chunk-1)*128 + 1)
    dataentryptr = pointer(tinyx.data[chunk])
    ccall(simdunpackFOR, Void, (UInt32, Ptr{m128i}, Ptr{UInt32}, UInt32), tinyx.offset[chunk], dataentryptr, xptr, tinyx.b[chunk])
    return x
end

# Overwrite x with an uncompressed chunk
function unpack!(tinyx::TinyIntVec, x::Vector{UInt32}, chunk::Integer)
    ccall(simdunpackFOR, Void, (UInt32, Ptr{m128i}, Ptr{UInt32}, UInt32), tinyx.offset[chunk], pointer(tinyx.data[chunk]), pointer(x), tinyx.b[chunk])
end

# Finds the location within a chunk where the i'th element resides
i2chunk_i(i::Integer) = (i - 1) & 127

# Finds the chunk where the i'th element resides
i2chunk(i::Integer) = (i - 1) >> 7 + 1


text:
using Test
using LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators
using ClimaCore.Geometry

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())
function hvspace_3D(
    xlim = (-π, π),
    ylim = (-π, π),
    zlim = (0, 4π),
    xelem = 4,
    yelem = 4,
    zelem = 16,
    npoly = 3,
)
    FT = Float64

    xdomain = Domains.IntervalDomain(
        Geometry.XPoint{FT}(xlim[1]),
        Geometry.XPoint{FT}(xlim[2]),
        periodic = true,
    )
    ydomain = Domains.IntervalDomain(
        Geometry.YPoint{FT}(ylim[1]),
        Geometry.YPoint{FT}(ylim[2]),
        periodic = true,
    )

    horzdomain = Domains.RectangleDomain(xdomain, ydomain)
    horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem)
    horztopology = Topologies.Topology2D(horzmesh)

    zdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_names = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(zdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

# set up 3D domain - doubly periodic box
hv_center_space, hv_face_space = hvspace_3D((-500, 500), (-500, 500), (0, 1000))

const MSLP = 1e5 # mean sea level pressure
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const γ = 1.4 # heat capacity ratio
const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
const C_v = R_d / (γ - 1) # heat capacity at constant volume
const R_m = R_d # moist R, assumed to be dry

function pressure(ρθ)
    if ρθ >= 0
        return MSLP * (R_d * ρθ / MSLP)^γ
    else
        return NaN
    end
end

Φ(z) = grav * z

# Reference: https://journals.ametsoc.org/view/journals/mwre/140/4/mwr-d-10-05073.1.xml, Section 5a
function init_dry_rising_bubble_3d(x, y, z)
    x_c = 0.0
    y_c = 0.0
    z_c = 350.0
    r_c = 250.0
    θ_b = 300.0
    θ_c = 0.5
    cp_d = C_p
    cv_d = C_v
    p_0 = MSLP
    g = grav

    # auxiliary quantities
    r = sqrt((x - x_c)^2 + (y - y_c)^2 + (z - z_c)^2)
    θ_p = r < r_c ? 0.5 * θ_c * (1.0 + cospi(r / r_c)) : 0.0 # potential temperature perturbation

    θ = θ_b + θ_p # potential temperature
    π_exn = 1.0 - g * z / cp_d / θ # exner function
    T = π_exn * θ # temperature
    p = p_0 * π_exn^(cp_d / R_d) # pressure
    ρ = p / R_d / T # density
    ρθ = ρ * θ # potential temperature density

    return (ρ = ρ, ρθ = ρθ)
end

# initial conditions
coords = Fields.coordinate_field(hv_center_space)
face_coords = Fields.coordinate_field(hv_face_space)

Yc = map(coord -> init_dry_rising_bubble_3d(coord.x, coord.y, coord.z), coords)
uₕ = map(_ -> Geometry.Covariant12Vector(0.0, 0.0), coords)
w = map(_ -> Geometry.Covariant3Vector(0.0), face_coords)
Y = Fields.FieldVector(Yc = Yc, uₕ = uₕ, w = w)

function energy(ρ, ρθ, ρu, z)
    u = ρu ./ ρ
    kinetic = ρ .* norm(u)^2 ./ 2
    potential = z .* grav .* ρ
    internal = C_v .* pressure.(ρθ) ./ R_d
    return kinetic .+ potential .+ internal
end
function combine_momentum(ρuₕ, ρw)
    Geometry.transform(Geometry.Covariant123Axis(), ρuₕ) +
    Geometry.transform(Geometry.Covariant123Axis(), ρw)
end
function center_momentum(Y)
    If2c = Operators.InterpolateF2C()
    ρ = Y.Yc.ρ
    ρuₕ = ρ .* Y.uₕ
    ρw = ρ .* If2c.(Y.w)
    combine_momentum.(ρuₕ, ρw)
end
function total_energy(Y)
    ρ = Y.Yc.ρ
    ρu = center_momentum(Y)
    ρθ = Y.Yc.ρθ
    z = Fields.coordinate_field(axes(ρ)).z
    sum(energy(ρ, ρθ, ρu, z))
end

energy_0 = total_energy(Y)
mass_0 = sum(Yc.ρ) # Computes ∫ρ∂Ω such that quadrature weighting is accounted for.
theta_0 = sum(Yc.ρθ) # Computes ∫ρ∂Ω such that quadrature weighting is accounted for.

function rhs_invariant!(dY, Y, _, t)

    cρ = Y.Yc.ρ # scalar on centers
    fw = Y.w # Covariant3Vector on faces
    cuₕ = Y.uₕ # Covariant12Vector on centers
    cρθ = Y.Yc.ρθ

    dρ = dY.Yc.ρ
    dw = dY.w
    duₕ = dY.uₕ
    dρθ = dY.Yc.ρθ

    # 0) update w at the bottom
    # fw = -g^31 cuₕ/ g^33

    hdiv = Operators.Divergence()
    hwdiv = Operators.Divergence()
    hgrad = Operators.Gradient()
    hwgrad = Operators.Gradient()
    hcurl = Operators.Curl()
    hwcurl = Operators.Curl()

    dρ .= 0 .* cρ

    ### HYPERVISCOSITY
    # 1) compute hyperviscosity coefficients

    χθ = @. dρθ = hwdiv(hgrad(cρθ / cρ)) # we store χθ in dρθ
    χuₕ = @. duₕ =
        hwgrad(hdiv(cuₕ)) - Geometry.Covariant12Vector(
            hwcurl(Geometry.Covariant3Vector(hcurl(cuₕ))),
        )

    Spaces.weighted_dss!(dρθ)
    Spaces.weighted_dss!(duₕ)

    κ₄ = 100.0 # m^4/s
    @. dρθ = -κ₄ * hwdiv(cρ * hgrad(χθ))
    @. duₕ =
        -κ₄ * (
            hwgrad(hdiv(χuₕ)) - Geometry.Covariant12Vector(
                hwcurl(Geometry.Covariant3Vector(hcurl(χuₕ))),
            )
        )

    # 1) Mass conservation
    If2c = Operators.InterpolateF2C()
    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    cw = If2c.(fw)
    fuₕ = Ic2f.(cuₕ)
    cuvw = Geometry.Covariant123Vector.(cuₕ) .+ Geometry.Covariant123Vector.(cw)

    dw .= fw .* 0

    # 1.a) horizontal divergence
    dρ .-= hdiv.(cρ .* (cuvw))

    # 1.b) vertical divergence
    vdivf2c = Operators.DivergenceF2C(
        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
    )
    # we want the total u³ at the boundary to be zero: we can either constrain
    # both to be zero, or allow one to be non-zero and set the other to be its
    # negation

    # explicit part
    dρ .-= vdivf2c.(Ic2f.(cρ .* cuₕ))
    # implicit part
    dρ .-= vdivf2c.(Ic2f.(cρ) .* fw)

    # 2) Momentum equation

    # curl term
    hcurl = Operators.Curl()
    # effectively a homogeneous Dirichlet condition on u₁ at the boundary
    vcurlc2f = Operators.CurlC2F(
        bottom = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
        top = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
    )
    cω³ = hcurl.(cuₕ) # Contravariant3Vector
    fω¹² = hcurl.(fw) # Contravariant12Vector
    fω¹² .+= vcurlc2f.(cuₕ) # Contravariant12Vector

    # cross product
    # convert to contravariant
    # these will need to be modified with topography
    fu¹² =
        Geometry.Contravariant12Vector.(
            Geometry.Covariant123Vector.(Ic2f.(cuₕ)),
        ) # Contravariant12Vector in 3D
    fu³ = Geometry.Contravariant3Vector.(Geometry.Covariant123Vector.(fw))
    @. dw -= fω¹² × fu¹² # Covariant3Vector on faces
    @. duₕ -= If2c(fω¹² × fu³)

    # Needed for 3D:
    @. duₕ -=
        cω³ × Geometry.Contravariant12Vector(Geometry.Covariant123Vector(cuₕ))

    cp = @. pressure(cρθ)
    @. duₕ -= hgrad(cp) / cρ
    vgradc2f = Operators.GradientC2F(
        bottom = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
        top = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
    )
    @. dw -= vgradc2f(cp) / Ic2f(cρ)

    cE = @. (norm(cuvw)^2) / 2 + Φ(coords.z)
    @. duₕ -= hgrad(cE)
    @. dw -= vgradc2f(cE)

    # 3) potential temperature

    @. dρθ -= hdiv(cuvw * cρθ)
    @. dρθ -= vdivf2c(fw * Ic2f(cρθ))
    @. dρθ -= vdivf2c(Ic2f(cuₕ * cρθ))

    fcc = Operators.FluxCorrectionC2C(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )

    @. dρ += fcc(fw, cρ)
    @. dρθ += fcc(fw, cρθ)
    # dYc.ρuₕ += fcc(w, Yc.ρuₕ)

    Spaces.weighted_dss!(dY.Yc)
    Spaces.weighted_dss!(dY.uₕ)
    Spaces.weighted_dss!(dY.w)

    return dY
end

dYdt = similar(Y);
rhs_invariant!(dYdt, Y, nothing, 0.0);

# run!
using OrdinaryDiffEq
Δt = 0.050
prob = ODEProblem(rhs_invariant!, Y, (0.0, 1.0))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 1.0,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "bubble_3d_invariant"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

# post-processing
Es = [total_energy(u) for u in sol.u]
Mass = [sum(u.Yc.ρ) for u in sol.u]
Theta = [sum(u.Yc.ρθ) for u in sol.u]

Plots.png(
    Plots.plot((Es .- energy_0) ./ energy_0),
    joinpath(path, "energy.png"),
)

Plots.png(Plots.plot((Mass .- mass_0) ./ mass_0), joinpath(path, "mass.png"))
Plots.png(Plots.plot((Theta .- theta_0) ./ theta_0), joinpath(path, ".png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(joinpath(path, "energy.png"), joinpath(@__DIR__, "../..")),
    "Total Energy",
)
linkfig(
    relpath(joinpath(path, "rho_theta.png"), joinpath(@__DIR__, "../..")),
    "Potential Temperature",
)
linkfig(
    relpath(joinpath(path, "mass.png"), joinpath(@__DIR__, "../..")),
    "Mass",
)


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::FiniteVariableIndex
    )::FiniteVariableRef
    return FiniteVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::VariableData{<:JuMP.ScalarVariable}
    )::FiniteVariableIndex
    return MOIUC.add_item(model.finite_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(
    model::InfiniteModel,
    ::Type{FiniteVariable}
    )::MOIUC.CleverDict{FiniteVariableIndex, VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}}
    return model.finite_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    vref::FiniteVariableRef
    )::MOIUC.CleverDict{FiniteVariableIndex, VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}}
    return JuMP.owner_model(vref).finite_vars
end

# Extend _data_object
function _data_object(
    vref::FiniteVariableRef
    )::VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}
  object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
  object === nothing && error("Invalid finite variable reference, cannot find " *
                        "corresponding variable in the model. This is likely " *
                        "caused by using the reference of a deleted variable.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(
    vref::FiniteVariableRef
    )::JuMP.ScalarVariable{Float64, Float64, Float64, Float64}
    return _data_object(vref).variable
end

################################################################################
#                             DEFINTION METHODS
################################################################################
## Process the variable to use the correct info
# All Floats
function _process_scalar_var(
    v::V
    )::V where {V <: JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}
    return v
end

# Other
function _process_scalar_var(
    v::JuMP.ScalarVariable
    )::JuMP.ScalarVariable{Float64, Float64, Float64, Float64}
    return JuMP.ScalarVariable(_make_float_info(v.info))
end

"""
    JuMP.add_variable(model::InfiniteModel, var::JuMP.ScalarVariable,
                      [name::String = ""])::GeneralVariableRef

Extend the `JuMP.add_variable` function to accomodate finite (scalar) variable
types. Adds a variable to an infinite model `model` and returns a
[`GeneralVariableRef`](@ref). Primarily intended to be an internal function of
the constructor macro `@variable`. However, it can be used in combination with
`JuMP.build_variable` to add finite variables to an infinite model object.

**Examples**
```julia-repl
julia> f_var = build_variable(error, info);

julia> fvref = add_variable(m, f_var, "var_name")
var_name
```
"""
function JuMP.add_variable(
    model::InfiniteModel,
    var::JuMP.ScalarVariable,
    name::String = ""
    )::GeneralVariableRef
    new_var = _process_scalar_var(var)
    data_object = VariableData(new_var, name)
    vindex = _add_data_object(model, data_object)
    vref = FiniteVariableRef(model, vindex)
    gvref = _make_variable_ref(model, vindex)
    _set_info_constraints(new_var.info, gvref, vref)
    model.name_to_var = nothing
    return gvref
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Set info for finite variables
function _update_variable_info(
    vref::FiniteVariableRef,
    info::JuMP.VariableInfo
    )::Nothing
    _set_core_variable_object(vref, JuMP.ScalarVariable(info))
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::FiniteVariableRef)::Nothing
    # remove variable info constraints associated with vref
    _delete_info_constraints(vref)
    return
end