text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for magnetic fields.
"""
abstract type MagneticFields{T} end

"""
$(SIGNATURES)

Return the magnetic field strength H [T/μ₀] (μ₀ being the vacuum permeability)
as `SVector{3,T}` for a given location r [m] and time t [s].
"""
function magneticFieldStrength(magneticField::MagneticFields{T},r::SVector{3,T},t::T) where {T}
    error("magneticFieldStrength is not implemented for magneticField::$(typeof(magneticField)).")
end

"""
$(TYPEDEF)

Abstract supertype for static magnetic fields.
"""
abstract type StaticMagneticFields{T} <: MagneticFields{T} end

"""
$(SIGNATURES)

Return the magnetic field strength H [T/μ₀] (μ₀ being the vacuum permeability)
as `SVector{3,T}` for a given location r [m].
"""
function magneticFieldStrength(magneticField::StaticMagneticFields{T},r::SVector{3,T}) where {T}
    error("magneticFieldStrength is not implemented for magneticField::$(typeof(magneticField)).")
end

function magneticFieldStrength(magneticField::StaticMagneticFields{T},r::SVector{3,T},t::T) where {T}
    return magneticFieldStrength(magneticField,r)
end

# Define specific subtypes
## ideal gradient field
struct StaticGradientField{T} <: StaticMagneticFields{T}
    center::SVector{3,T}
    gradient::SMatrix{3,3,T}

    StaticGradientField(center::AbstractVector{T},gradient::AbstractMatrix{T}) where {T} = new{T}(SVector{3,T}(center),SMatrix{3,3,T}(gradient))
end

function magneticFieldStrength(magneticField::StaticGradientField{T},r::SVector{3,T}) where {T}
    return magneticField.gradient*(r-magneticField.center)
end

## ideal homogenous field
struct StaticHomogeneousField{T} <: StaticMagneticFields{T}
    magneticFieldStrength::SVector{3,T}

    StaticHomogeneousField(magneticFieldStrength::Vector{T}) where {T} = new{T}(SVector{3,T}(magneticFieldStrength))
end

function magneticFieldStrength(magneticField::StaticHomogeneousField{T},r::SVector{3,T}) where {T}
    return magneticField.magneticFieldStrength
end

## ideal homogenous field with sinusoidal 1D excitation
struct SinusoidalField{T} <: MagneticFields{T}
    spatialFieldProfile::StaticMagneticFields{T}
    baseFrequency::T
    divider::Int
    doubleFrequency::T
    relativePhase::T # = phase/π
end

function SinusoidalField(spatialFieldProfile::StaticMagneticFields{T},baseFrequency::Real,divider::Real,relativePhase::Real=0.0) where {T}
    doubleFrequency = 2//divider * baseFrequency
    return SinusoidalField{T}(spatialFieldProfile,T(baseFrequency),T(divider),T(doubleFrequency),T(relativePhase))
end

function magneticFieldStrength(magneticField::SinusoidalField{T},r::SVector{3,T},t::T) where {T}
    return magneticFieldStrength(magneticField.spatialFieldProfile,r)*sinpi(magneticField.doubleFrequency*t+magneticField.relativePhase)
end


text:
module Model
using Random
using Distributions:Normal, Bernoulli, logpdf

using Flux

function conv_MINST_model(latent_size::Int)
  """
  Convolutional Encoder (g) and Decoder (f)
  which have the following properties:
  X ~ f(reparameterize.(split_encoder_result(g(X))))

  """
  g = Chain(
    Conv((3, 3), 1 => 32, relu, stride = (2, 2)),
    Conv((3, 3), 32 => 64, relu, stride = (2, 2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(6 * 6 * 64, latent_size * 2)
  )
  f = Chain(
    Dense(latent_size, 7 * 7 * 32, relu),
    x -> reshape(x, 7, 7, 32, size(x, 2)),
    ConvTranspose((3, 3), 32 => 64, relu, stride = (2, 2), pad = (0,0)),
    ConvTranspose((3, 3), 64 => 32, relu, stride = (2,2), pad = (0,0)),
    ConvTranspose((3, 3), 32 => 1, stride = (1, 1), pad = (2,3,2,3)),
  )
  return f, g
end

function model_sample(f)
  """
  model_sample(f)
  where f is a Flux Chain of operations
  and the first operation is a Dense Layer
  """
  latent_size = size(f[1].W)[2]
  zero_input = zeros(Float32, latent_size)
  rand.(Bernoulli.(sigmoid.(f(reparameterize.(zero_input, zero_input)))))
end

# overloaded version encodes/decodes X_samples
function model_sample(g, f, X_samples)
  xs = reshape(X_samples[1:28, 1:28, 1, 1], 28, 28, 1, 1)
  (μ̂, logσ̂) = split_encoder_result(g(xs))
  #rand.(Bernoulli.(sigmoid.(f(reparameterize.(μ̂, logσ̂)))))
  sigmoid.(f(reparameterize.(μ̂, logσ̂)))
end

function split_encoder_result(X, n_latent::Int64)
  μ = X[1:n_latent, :]
  logσ = X[(n_latent + 1):(n_latent * 2), :]
  return μ, logσ
end

function split_encoder_result(X)
  n_latent = convert(Int, floor(size(X,1) / 2))
  μ = X[1:n_latent, :]
  logσ = X[(n_latent + 1):(n_latent * 2), :]
  return μ, logσ
end

function reparameterize(μ :: T, logσ :: T) where {T}
  """
  reparameterize the results of 'encoder'
  onto a Normal(0,1) distribution,
  we could also pass "Float32" here instead, which would
  provide a uniform Distribution
  """
  return rand(Normal(0,1)) * exp(logσ * 0.5f0) + μ
end

function kl_q_p(μ :: T, logσ :: T) where {T}
  """"
  KL-divergence divergence, between approximate posterior/prior
  This is the KLD between a
  1) distribution of with means μ, variance logσ, and importanty, an
  identity covariance matrix, and
  2) A standard normal distribution Normal(μ = 0,σ = 1)
  """
  return 0.5f0 * sum(exp.(2f0 .* logσ) + μ.^2 .- 1 .- (2 .* logσ))
end

function sigmoid(z)
  return 1.0f0 ./ (1.0f0 .+ exp(-z))
end

function logp_x_z(x, z, f)
  """
  calculates log  P( X = x | Z = z ), reconstruction error
  or: the conditional probability of data(digit image) given latent(z).
  Our output is an image, but we are creating a distribution
  using the sigmoid function to normalize decoder outputs
  then passing it into a Bernoulli distribution.
  logpdf(Distribution from decodor, image pixel) ...
  => logpdf(Bernoulli(0.5f0),{0,1}) -> {-0.693, -0.693}
  => summate over all of these values
 """
  return sum(logpdf.(Bernoulli.(sigmoid.(f(z))), x))
end

# Monte Carlo estimator of mean ELBO using M samples.
function L̄(X, g, f, M :: Int64)
  (μ̂, logσ̂) = split_encoder_result(g(X))
  return (logp_x_z(X, reparameterize.(μ̂, logσ̂), f) - kl_q_p(μ̂, logσ̂)) * 1 // M
end

# returns loss function for X with input g/f M
function build_loss(g, f, M)
  return X -> (-L̄(X, g, f, M) + 0.01f0 * sum(x->sum(x.^2), params(f)))
end

## Interface
function create_vae(latent_size :: Int64, M :: Int64)
  f, g = conv_MINST_model(latent_size)
  ps = params(f, g)
  loss = build_loss(g, f, M)
  return ps, loss, f, g
end

end

ext: jl

max_stars_repo_name: aniketjivani/variational-autoencoders

max_issues_repo_name: aniketjivani/variational-autoencoders

max_forks_repo_name: aniketjivani/variational-autoencoders


text:
###########################################################################
#                 Extensions of BranchAndBound Functions                  #
#                   General Algorithm Flow Goes here                      #
###########################################################################

# return a "deepcopy" of SpatialBCBranch for everything except for variable references
# this is needed to ensure variable references in new branch copy are still valid
# the mod_branch will not be duplicated here if it is not nothing
function branch_copy(branch::SpatialBCBranch)
    return SpatialBCBranch(branch.i, branch.j, branch.wii, branch.wjj, branch.wr, branch.wi,
                           branch.mod_branch, deepcopy(branch.bounds), deepcopy(branch.valid_ineq_coeffs))
end

# this uses trace minimization to test if the given incomplete wr, wi can recover a rank-1 complex matrix
# if yes, wr and wi will be updated with the complete rank-1 solution, and the function returns true
# if no, wr and wi will remain unchanged, and the function returns false
# TODO: modify this so that it takes into account cut lines if necessary
#=
function min_trace_solution!(wr::Matrix{Float64}, wi::Matrix{Float64}, pm::AbstractPowerModel)::Bool
    bus_ids = ids(pm, :bus)
    n_bus = length(bus_ids)
    lookup_w_index = Dict((bi, i) for (i, bi) in enumerate(bus_ids))

    @info "--- Eigenvalues of W solution from SDP relaxation: " svdvals(wr + im .* wi)[1:3]

    # set up trace minimization
    m = Model(MOSEK_OPTIMIZER)
    @variable(m, XR[1:n_bus, 1:n_bus])
    @variable(m, XI[1:n_bus, 1:n_bus])
    @constraint(m, [XR -XI; XI XR] in PSDCone())
    for i in bus_ids
        @constraint(m, XR[i,i] == wr[i,i])
    end
    for (i,j) in ids(pm, :buspairs)
        windex_i = lookup_w_index[i]
        windex_j = lookup_w_index[j]
        @constraint(m, XR[i,j] == wr[i,j])
        @constraint(m, XI[i,j] == wi[i,j])
        @constraint(m, XR[j,i] == wr[j,i])
        @constraint(m, XI[j,i] == wi[j,i])
    end
    # TODO: add codes here to take into account cut lines if necessary (for network decomposition later)
    @objective(m, Min, tr([XR -XI; XI XR]))
    optimize!(m)
    @info "--- Trace minimization termination status: $(JuMP.termination_status(m)) ---"

    # eigenvalue test
    sol = value.(XR) + im * value.(XI)
    # @info "--- Trace minimization solution: " sol
    rk = rank(sol, rtol = 1e-4)
    @info "--- Eigenvalues of W solution after trace minimization: " svdvals(sol)[1:3]
    @info "--- The rank of the trace minimization solution: $(rk) ---"
    if rk == 1
        wr = value.(XR)
        wi = value.(XI)
        return true
    end
    return false
end
=#
# simply check the w values returned by the solver, no special handling of free variables
function min_trace_solution!(wr::Matrix{Float64}, wi::Matrix{Float64}, pm::PM.AbstractPowerModel)::Bool
    solution = wr + im .* wi
    # @info "--- Eigenvalues: " svdvals(solution)
    return rank(solution, rtol = 1e-4) == 1
end

function BB.branch!(tree::BB.AbstractTree, node::BB.AbstractNode)
    @info " Node id $(node.id), status $(node.solution_status), bound $(node.bound)"
    root = find_root(node)
    model = root.model
    if node.bound >= tree.best_incumbent
        if isapprox(node.bound, tree.best_incumbent, rtol = 1e-5) && node.bound > root.auxiliary_data["best_bound"]
            root.auxiliary_data["best_bound"] = node.bound
            root.auxiliary_data["best_id"] = node.id
        end
        @info " Fathomed by bound"
    elseif node.depth >= 100
        if node.bound > root.auxiliary_data["best_bound"]
            root.auxiliary_data["best_bound"] = node.bound
            root.auxiliary_data["best_id"] = node.id
        end
        @info " Fathomed by maximum depth"
    elseif node.solution_status in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.SLOW_PROGRESS, MOI.ALMOST_OPTIMAL]
        # determine the complex entry to branch on based on solution
        pm = root.auxiliary_data["PM"]
        (i,j) = find_min_eigen(node)
        WR = PM.var(pm, :WR)
        WI = PM.var(pm, :WI)
        wr = JuMP.value.(WR)
        wi = JuMP.value.(WI)
        is_rank1 = min_trace_solution!(wr, wi, root.auxiliary_data["PM"]) # this implements trace minimization for rank-1 check and recovery
        # eigvs = eigvals(w)
        # println(eigvs)
        # if node.auxiliary_data["eigenvalues"][(i,j)] <= 1e-8
        # if eigvs[end-1] <= 5e-5
        if is_rank1
            # update the node with best bound (including processed ones), recorded in root
            root.auxiliary_data["rank1"][node.id] = node.bound
            @info " Fathomed by reaching rank-1 solution"
        else
            new_sbc_branch = create_sbc_branch(i, j, node)
            (new_i,new_j) = find_branching_entry(new_sbc_branch, node)

            # create branches and child nodes accordingly
            up_bounds_arr = [k for k in new_sbc_branch.bounds]
            down_bounds_arr = [k for k in new_sbc_branch.bounds]
            bus_ids = PM.ids(pm, :bus)
            lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))
            widx_new_i = lookup_w_index[new_i]
            widx_new_j = lookup_w_index[new_j]
            if new_i != new_j # branch on Wij
                up_bounds_arr[5] = (up_bounds_arr[5] + up_bounds_arr[6]) / 2
                down_bounds_arr[6] = (down_bounds_arr[5] + down_bounds_arr[6]) / 2
                vpair = (PM.var(pm, :WR)[widx_new_i,widx_new_j], PM.var(pm, :WI)[widx_new_i,widx_new_j])
                up_mod_branch = ComplexVariableBranch(Dict(vpair => up_bounds_arr[5]), Dict(vpair => up_bounds_arr[6]))
                down_mod_branch = ComplexVariableBranch(Dict(vpair => down_bounds_arr[5]), Dict(vpair => down_bounds_arr[6]))
                @info " Branch at W$(new_i)$(new_j), [L, U] breaks into by [$(down_bounds_arr[5]),$(up_bounds_arr[5]),$(up_bounds_arr[6])]."
            elseif new_i == i # branch on Wii
                up_bounds_arr[1] = (up_bounds_arr[1] + up_bounds_arr[2]) / 2
                down_bounds_arr[2] = (down_bounds_arr[1] + down_bounds_arr[2]) / 2
                v = PM.var(pm, :WR)[widx_new_i, widx_new_j]
                up_mod_branch = BB.VariableBranch(Dict(v => up_bounds_arr[1]), Dict(v => up_bounds_arr[2]))
                down_mod_branch = BB.VariableBranch(Dict(v => down_bounds_arr[1]), Dict(v => down_bounds_arr[2]))
                @info " Branch at W$(new_i)$(new_j), [L, U] breaks into by [$(down_bounds_arr[1]),$(up_bounds_arr[1]),$(up_bounds_arr[2])]."
            else # branch on Wjj
                up_bounds_arr[3] = (up_bounds_arr[3] + up_bounds_arr[4]) / 2
                down_bounds_arr[4] = (down_bounds_arr[3] + down_bounds_arr[4]) / 2
                v = PM.var(pm, :WR)[widx_new_i, widx_new_j]
                up_mod_branch = BB.VariableBranch(Dict(v => up_bounds_arr[3]), Dict(v => up_bounds_arr[4]))
                down_mod_branch = BB.VariableBranch(Dict(v => down_bounds_arr[3]), Dict(v => down_bounds_arr[4]))
                @info " Branch at W$(new_i)$(new_j), [L, U] breaks into by [$(down_bounds_arr[3]),$(up_bounds_arr[3]),$(up_bounds_arr[4])]."
            end
            next_branch_up = branch_copy(new_sbc_branch)
            next_branch_down = new_sbc_branch
            next_branch_up.bounds = Tuple(up_bounds_arr)
            next_branch_down.bounds = Tuple(down_bounds_arr)
            next_branch_up.mod_branch = up_mod_branch
            next_branch_down.mod_branch = down_mod_branch
            child_up = BB.create_child_node(node, next_branch_up)
            child_down = BB.create_child_node(node, next_branch_down)
            BB.push!(tree, child_up)
            BB.push!(tree, child_down)
        end
    else
        @info " Fathomed by solution status: $(node.solution_status)"
    end
    push!(root.auxiliary_data["best_bounds"], root.auxiliary_data["best_bound"])
    delete_prev_branch_constr!(model, node)
end

# implement depth first rule
function BB.next_node(tree::BB.AbstractTree)
    # sort! is redefined to find the node with maximum depth (consistent with the paper's implementation)
    # sort!(tree::BB.AbstractTree) = Base.sort!(tree.nodes, by=x->x.depth)
    BB.sort!(tree)
    node = Base.pop!(tree.nodes)
    return node
end

function BB.termination(tree::BB.AbstractTree)
    @info "Tree nodes: processed $(length(tree.processed)), left $(length(tree.nodes)), total $(tree.node_counter), best bound $(tree.best_bound), best incumbent $(tree.best_incumbent)"
    if BB.isempty(tree)
        # @info "Completed the tree search"
        return true
    end
    if length(tree.processed) >= MAX_NODE_ALLOWED
        # @info "Reached node limit"
        return true
    end
    return false
end

# This makes original adjust_branch! invalid
BB.adjust_branch!(branch_objects::Array{SpatialBCBranch,1}) = nothing

# This makes original apply_changes! invalid
BB.apply_changes!(node::BB.JuMPNode) = nothing

# This makes original processed! invalid
function BB.processed!(tree::BB.AbstractTree, node::BB.JuMPNode)
    Base.push!(tree.processed, node)
end

function BB.bound!(node::BB.JuMPNode)
    model = find_root(node).model
    backtracking!(model, node)
    # The commented part is used for debugging JuMP.optimize!(model) in environment with MosekTools@0.9.4
    # This bug does not appear with MosekTools@0.9.3
    # try
    #     JuMP.optimize!(model)
    # catch e
    #     println(Mosek.getbarxj(model.moi_backend.optimizer.model.task, Mosek.MSK_SOL_ITR, 1)) # MSK_SOL_ITR gives back the solution ... somehow MOSEK goes through MSK_SOL_ITG which looks for the integer solution
    #     println(Mosek.solutiondef(model.moi_backend.optimizer.model.task, Mosek.MSK_SOL_ITR))
    #     println(Mosek.solutiondef(model.moi_backend.optimizer.model.task, Mosek.MSK_SOL_ITG)) # This unexpectedly returns true @0.9.4
    #     println(Mosek.solutiondef(model.moi_backend.optimizer.model.task, Mosek.MSK_SOL_BAS))
    #     JuMP.write_to_file(model, "/home/weiqizhang/anl/mosek_test_model", format = MOI.FileFormats.FORMAT_SDPA)
    # end
    JuMP.optimize!(model)
    node.solution_status = JuMP.termination_status(model)
    if node.solution_status == MOI.INFEASIBLE || JuMP.dual_status(model) in [MOI.INFEASIBILITY_CERTIFICATE]
        node.bound = Inf
    elseif node.solution_status == MOI.DUAL_INFEASIBLE || JuMP.primal_status(model) in [MOI.INFEASIBILITY_CERTIFICATE]
        node.bound = -Inf
    elseif node.solution_status in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.SLOW_PROGRESS, MOI.ALMOST_OPTIMAL]
        node.bound = max(node.bound, JuMP.objective_value(model))
        vrefs = JuMP.all_variables(model)
        for v in vrefs
            node.solution[v] = JuMP.value(v)
        end
    else
        @warn "Unexpected node solution status: $(node.solution_status)"
        node.bound = -Inf
    end

    # store the solution in node
    solution = Dict{String, Float64}() # name => value
    for var in JuMP.all_variables(model)
        solution[JuMP.name(var)] = JuMP.value(var)
    end
    node.auxiliary_data["node solution"] = solution

end

# This finds an incumbent for each node, by converting SDP problem into rank-1 form and solving with Ipopt
#=
function BB.heuristics!(node::BB.JuMPNode)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
    root = find_root(node)
    tree = root.auxiliary_data["tree"]
    pm = root.auxiliary_data["PM"]
    exact_model = copy(root.model)
    buses = collect(keys(ref(pm, :bus)))
    JuMP.@variable(exact_model, vr[buses])
    JuMP.@variable(exact_model, vi[buses])

    # get rid of PSDCone constraint (for Ipopt)
    psdcone = JuMP.all_constraints(exact_model, Array{GenericAffExpr{Float64, VariableRef}, 1}, MOI.PositiveSemidefiniteConeSquare)[1]
    JuMP.delete(exact_model, psdcone)

    # add rank 1 constraints
    bus_ids = ids(pm, :bus)
    lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))
    for i in buses, j in buses
        wi = lookup_w_index[i]
        wj = lookup_w_index[j]
        if i == j
            Wi = JuMP.variable_by_name(exact_model, "0_WR[$(i),$(i)]")
            Wj = JuMP.variable_by_name(exact_model, "0_WR[$(j),$(j)]")
            JuMP.@constraint(exact_model, Wi == vr[i]^2 + vi[i]^2)
            JuMP.@constraint(exact_model, Wj == vr[j]^2 + vi[j]^2)
        else
            if !isnothing(JuMP.variable_by_name(exact_model, "0_WR[$(i),$(j)]"))
                WR = JuMP.variable_by_name(exact_model, "0_WR[$(i),$(j)]")
            else
                WR = JuMP.variable_by_name(exact_model, "0_WR[$(j),$(i)]")
            end
            WI = JuMP.variable_by_name(exact_model, "0_WI[$(i),$(j)]")
            JuMP.@constraint(exact_model, WR == vr[i] * vr[j] + vi[i] * vi[j])
            JuMP.@constraint(exact_model, WI == vr[j] * vi[i] - vi[j] * vr[i])
        end
    end

    # convert second order cone constraints into quadratic constraints
    soc_crefs = JuMP.all_constraints(exact_model, Array{GenericAffExpr{Float64, VariableRef}, 1}, MOI.SecondOrderCone)
    for cref in soc_crefs
        cobj = JuMP.constraint_object(cref)
        terms = cobj.func
        JuMP.@constraint(exact_model, sum(terms[i]^2 for i in 2:length(terms)) <= terms[1]^2)
        JuMP.delete(exact_model, cref)
    end

    # convert rotated second order cone constraints into quadratic constraints (this is needed if objective is quadratic)
    rsoc_crefs = JuMP.all_constraints(exact_model, Array{GenericAffExpr{Float64, VariableRef}, 1}, MOI.RotatedSecondOrderCone)
    for cref in rsoc_crefs
        cobj = JuMP.constraint_object(cref)
        terms = cobj.func
        JuMP.@constraint(exact_model, sum(terms[i]^2 for i in 3:length(terms)) <= 2 * terms[1].constant * terms[2]) # hacky way to avoid higher order (nonlinear) terms
        JuMP.delete(exact_model, cref)
    end

    # initialize with sdp solution
    # sdp_sol = node.auxiliary_data["node solution"]
    # for (name, val) in sdp_sol
    #     var = JuMP.variable_by_name(exact_model, name)
    #     JuMP.set_start_value(var, val)
    # end

    # initialize with PowerModels default for ACRPowerModel
    JuMP.set_start_value.(vr, 1.0)

    JuMP.set_optimizer(exact_model, ipopt_optimizer)
    JuMP.optimize!(exact_model)
    @info " Heuristics result: $(objective_value(exact_model))"
    exact_term_status = JuMP.termination_status(exact_model)

    if exact_term_status != MOI.LOCALLY_SOLVED
        @info "  Warning: node exact model is not solved to local optimality, only to $(exact_term_status). "
    end

    if exact_term_status == MOI.LOCALLY_SOLVED && tree.best_incumbent >= JuMP.objective_value(exact_model)
        @info "  Updating tree best incumbent. Exact model termination status: $(termination_status(exact_model))"
        tree.best_incumbent = JuMP.objective_value(exact_model)
    end
end
=#

function BB.heuristics!(node::BB.JuMPNode)
    ipopt_optimizer = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
    root = find_root(node)
    tree = root.auxiliary_data["tree"]
    pm = root.auxiliary_data["PM"]
    exact_model = copy(root.model)
    buses = collect(keys(PM.ref(pm, :bus)))
    JuMP.@variable(exact_model, vr[buses])
    JuMP.@variable(exact_model, vi[buses])

    # get rid of PSDCone constraint
    psdcone = JuMP.all_constraints(exact_model, Array{JuMP.GenericAffExpr{Float64, JuMP.VariableRef}, 1}, MOI.PositiveSemidefiniteConeSquare)[1]
    JuMP.delete(exact_model, psdcone)

    # replace all w variables with voltage product pairs
    wr = PM.var(pm, :wr)
    wr_names = Dict(JuMP.name(var) => i for (i, var) in wr)
    for (i, var) in PM.var(pm, :w)
        wr_names[JuMP.name(var)] = (i, i)
    end
    wi = PM.var(pm, :wi)
    wi_names = Dict(JuMP.name(var) => i for (i, var) in wi)
    # bus_ids = ids(pm, :bus)
    # lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))
    for (functype, settype) in [ (JuMP.GenericAffExpr{Float64,JuMP.VariableRef}, MOI.EqualTo{Float64}),
                                 (JuMP.GenericAffExpr{Float64,JuMP.VariableRef}, MOI.GreaterThan{Float64}),
                                 (JuMP.GenericAffExpr{Float64,JuMP.VariableRef}, MOI.LessThan{Float64})]
        #
        crefs = JuMP.all_constraints(exact_model, functype, settype)
        for cref in crefs
            cobj = JuMP.constraint_object(cref)
            original_lhs = cobj.func
            new_lhs = zero(JuMP.QuadExpr)
            cref_modified = false
            for (var, coeff) in cobj.func.terms
                var_name = JuMP.name(var)
                if var_name in keys(wr_names)
                    i, j = wr_names[var_name]
                    new_lhs += coeff * (vr[i] * vr[j] + vi[i] * vi[j])
                    original_lhs -= coeff * var
                    cref_modified = true
                elseif var_name in keys(wi_names)
                    i, j = wi_names[var_name]
                    new_lhs += coeff * (vi[i] * vr[j] - vr[i] * vi[j])
                    original_lhs -= coeff * var
                    cref_modified = true
                end
            end
            if cref_modified
                JuMP.drop_zeros!(original_lhs)
                JuMP.@constraint(exact_model, original_lhs + new_lhs in cobj.set)
                JuMP.delete(exact_model, cref)
            end
        end
    end
    for (functype, settype) in [ (JuMP.VariableRef, MOI.GreaterThan{Float64}),
                                 (JuMP.VariableRef, MOI.LessThan{Float64})]
        crefs = JuMP.all_constraints(exact_model, functype, settype)
        for cref in crefs
            cobj = JuMP.constraint_object(cref)
            new_lhs = zero(JuMP.QuadExpr)
            cref_modified = false
            var = cobj.func
            var_name = JuMP.name(var)
            if var_name in PM.keys(wr_names)
                i, j = wr_names[var_name]
                new_lhs = (vr[i] * vr[j] + vi[i] * vi[j])
                cref_modified = true
            elseif var_name in PM.keys(wi_names)
                i, j = wi_names[var_name]
                new_lhs = (vi[i] * vr[j] - vr[i] * vi[j])
                cref_modified = true
            end
            if cref_modified
                JuMP.@constraint(exact_model, new_lhs in cobj.set)
                JuMP.delete(exact_model, cref)
            end
        end
    end

    # delete W variables
    # for var in all_variables(exact_model)
    #     if occursin("WR", JuMP.name(var)) || occursin("WI", JuMP.name(var))
    #         delete(exact_model, var)
    #     end
    # end

    # convert second order cone constraints into quadratic constraints
    soc_crefs = JuMP.all_constraints(exact_model, Array{JuMP.GenericAffExpr{Float64, JuMP.VariableRef}, 1}, MOI.SecondOrderCone)
    for cref in soc_crefs
        cobj = JuMP.constraint_object(cref)
        terms = cobj.func
        JuMP.@constraint(exact_model, sum(terms[i]^2 for i in 2:length(terms)) <= terms[1]^2)
        JuMP.delete(exact_model, cref)
    end

    # convert rotated second order cone constraints into quadratic constraints (this is needed if objective is quadratic)
    rsoc_crefs = JuMP.all_constraints(exact_model, Array{JuMP.GenericAffExpr{Float64, JuMP.VariableRef}, 1}, MOI.RotatedSecondOrderCone)
    for cref in rsoc_crefs
        cobj = JuMP.constraint_object(cref)
        terms = cobj.func
        JuMP.@constraint(exact_model, sum(terms[i]^2 for i in 3:length(terms)) <= 2 * terms[1].constant * terms[2]) # hacky way to avoid higher order (nonlinear) terms
        JuMP.delete(exact_model, cref)
    end

    # initialize with sdp solution
    # sdp_sol = node.auxiliary_data["node solution"]
    # for (name, val) in sdp_sol
    #     var = JuMP.variable_by_name(exact_model, name)
    #     if !isnothing(var)
    #         JuMP.set_start_value(var, val)
    #     end
    # end

    # initialize with PowerModels default for ACRPowerModel
    JuMP.set_start_value.(vr, 1.0)

    JuMP.set_optimizer(exact_model, ipopt_optimizer)
    JuMP.optimize!(exact_model)
    @info " Heuristics result: $(JuMP.objective_value(exact_model))"
    exact_term_status = JuMP.termination_status(exact_model)

    if exact_term_status != MOI.LOCALLY_SOLVED
        @info "  Warning: node exact model is not solved to local optimality, only to $(exact_term_status). "
    end

    if exact_term_status == MOI.LOCALLY_SOLVED && tree.best_incumbent >= JuMP.objective_value(exact_model)
        @info "  Updating tree best incumbent. Exact model termination status: $(JuMP.termination_status(exact_model))"
        tree.best_incumbent = JuMP.objective_value(exact_model)
    end
end

max_stars_repo_name: kibaekkim/ComplexBC

max_issues_repo_name: kibaekkim/ComplexBC

max_forks_repo_name: kibaekkim/ComplexBC


text:
export CyclotomicExt, cyclotomic_extension

mutable struct CyclotomicExt
  k::AnticNumberField
  n::Int
  Kr::Hecke.NfRel
  Ka::AnticNumberField
  mp::Any
  function CyclotomicExt()
    return new()
  end
end

function Base.show(io::IO, c::CyclotomicExt)
  print(io, "Cyclotomic extension by zeta_$(c.n) of degree $(degree(c.Ka))")
end

function cyclotomic_extension(k::AnticNumberField, n::Int)
  Ac = CyclotomicExt[]
  try
    Ac = Hecke._get_cyclotomic_ext_nf(k)
    for i = Ac
      if i.n == n
        return i
      end
    end
  catch e
    if e != AbstractAlgebra.AccessorNotSetError()
      rethrow(e)
    end
  end

  ZX, X = PolynomialRing(FlintZZ)
  Qx = parent(k.pol)
  kt, t = PolynomialRing(k, "t", cached = false)
  f = cyclotomic(n, X)
  fq = Qx(f)
  fk = evaluate(fq, t)
  lf = factor(fk)
  fk = first(keys(lf.fac))

  Kr, Kr_gen = number_field(fk, "z_$n", false)
  Ka, rel2abs, small2abs = Hecke.absolute_field(Kr, false)

  c = CyclotomicExt()
  c.k = k
  c.n = n
  c.Kr = Kr
  c.Ka = Ka
  c.mp = (rel2abs, small2abs)

  b_k = basis(maximal_order(k), k)
  B_k = [small2abs(x) for x = b_k]
  g = rel2abs(Kr_gen)
  for i=1:degree(fk)-1
    for j=1:degree(k)
      push!(B_k, B_k[j]*g^i)
    end
  end

  # we argued that an equation order defined from a factor of a
  # cyclotomic polynomial should always be maximal by Dedekind
  # hence the product basis should be maximal...
  # we can later implement proper relative maximal orders
  # wrong! for Q(sqrt(10))(zeta_5) this is not 5 maximal.
  # so we need to become p-maximal for all p dividing the gcd()

  ZKa = Hecke.NfOrd(B_k)
  for (p,v) = factor(gcd(discriminant(maximal_order(k)), fmpz(n))).fac
    ZKa = pmaximal_overorder(ZKa, p)
  end
  ZKa = lll(ZKa)
  ZKa.ismaximal = 1
  Hecke._set_maximal_order_of_nf(Ka, ZKa)
  push!(Ac, c)
  Hecke._set_cyclotomic_ext_nf(k, Ac)
  return c
end


text:
"""
Nodal admittance matrix (Ybus) is an N x N matrix describing a power system with N buses. It represents the nodal admittance of the buses in a power system.

The Ybus Struct is indexed using the Bus Numbers, no need for them to be sequential
"""
struct Ybus{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF64}
    data::SparseArrays.SparseMatrixCSC{ComplexF64, Int}
    axes::Ax
    lookup::L
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    br::ACBranch,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(br)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_l = (1 / (get_r(br) + get_x(br) * 1im))
    Y11 = Y_l + (1im * get_b(br).from)
    if !isfinite(Y11) || !isfinite(Y_l)
        error("Data in $(get_name(br)) is incorrect. r = $(get_r(br)), x = $(get_x(br))")
    end
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = -Y_l
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_l + (1im * get_b(br).to)
    ybus[bus_to_no, bus_to_no] += Y22
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    br::Transformer2W,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(br)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_t = 1 / (get_r(br) + get_x(br) * 1im)
    Y11 = Y_t
    b = get_primary_shunt(br)
    if !isfinite(Y11) || !isfinite(Y_t) || !isfinite(b)
        error("Data in $(get_name(br)) is incorrect. r = $(get_r(br)), x = $(get_x(br))")
    end
    ybus[bus_from_no, bus_from_no] += Y11
    ybus[bus_from_no, bus_to_no] += -Y_t
    ybus[bus_to_no, bus_from_no] += -Y_t
    ybus[bus_to_no, bus_to_no] += Y_t + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    br::TapTransformer,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(br)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (get_r(br) + get_x(br) * 1im)
    c = 1 / get_tap(br)
    b = get_primary_shunt(br)

    Y11 = (Y_t * c^2)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t * c)
    if !isfinite(Y11) || !isfinite(Y12) || !isfinite(b)
        error("Data in $(get_name(br)) is incorrect. r = $(get_r(br)), x = $(get_x(br))")
    end
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    br::PhaseShiftingTransformer,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(br)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_t = 1 / (get_r(br) + get_x(br) * 1im)
    tap = (get_tap(br) * exp(get_α(br) * 1im))
    c_tap = (get_tap(br) * exp(-1 * get_α(br) * 1im))
    b = get_primary_shunt(br)
    Y11 = (Y_t / abs(tap)^2)
    if !isfinite(Y11) || !isfinite(Y_t) || !isfinite(b * c_tap)
        error("Data in $(get_name(br)) is incorrect. r = $(get_r(br)), x = $(get_x(br))")
    end
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t / c_tap)
    ybus[bus_from_no, bus_to_no] += Y12
    Y21 = (-Y_t / tap)
    ybus[bus_to_no, bus_from_no] += Y21
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    fa::FixedAdmittance,
    num_bus::Dict{Int, Int},
)
    bus = get_bus(fa)
    bus_no = num_bus[get_number(bus)]
    ybus[bus_no, bus_no] += fa.Y
    return
end

function _buildybus(branches, nodes, fixed_admittances)
    buscount = length(nodes)
    num_bus = Dict{Int, Int}()

    for (ix, b) in enumerate(nodes)
        num_bus[get_number(b)] = ix
    end
    ybus = SparseArrays.spzeros(ComplexF64, buscount, buscount)
    for b in branches
        if get_name(b) == "init"
            throw(DataFormatError("The data in Branch is invalid"))
        end
        get_available(b) && _ybus!(ybus, b, num_bus)
    end
    for fa in fixed_admittances
        get_available(fa) && _ybus!(ybus, fa, num_bus)
    end
    return SparseArrays.dropzeros!(ybus)
end

"""
Builds a Ybus from a collection of buses and branches. The return is a Ybus Array indexed with the bus numbers and the branch names.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
- `connectivity_method::Function = goderya_connectivity`: method (`goderya_connectivity` or `dfs_connectivity`) for connectivity validation
"""
function Ybus(
    branches,
    nodes,
    fixed_admittances = Vector{FixedAdmittance}();
    check_connectivity::Bool = true,
    kwargs...,
)
    nodes = sort!(collect(nodes), by = x -> get_number(x))
    bus_ax = get_number.(nodes)
    axes = (bus_ax, bus_ax)
    bus_lookup = _make_ax_ref(bus_ax)
    look_up = (bus_lookup, bus_lookup)
    ybus = _buildybus(branches, nodes, fixed_admittances)
    if check_connectivity
        connected = validate_connectivity(ybus, nodes, bus_lookup; kwargs...)
        !connected && throw(DataFormatError("Network not connected"))
    end
    return Ybus(ybus, axes, look_up)
end

"""
Builds a Ybus from the system. The return is a Ybus Array indexed with the bus numbers and the branch names.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
- `connectivity_method::Function = goderya_connectivity`: method (`goderya_connectivity` or `dfs_connectivity`) for connectivity validation
"""
function Ybus(sys::System; check_connectivity::Bool = true, kwargs...)
    branches = get_components(ACBranch, sys)
    nodes = get_components(Bus, sys)
    fixed_admittances = get_components(FixedAdmittance, sys)
    return Ybus(
        branches,
        nodes,
        fixed_admittances;
        check_connectivity = check_connectivity,
        kwargs...,
    )
end

"""
Nodal incidence matrix (Adjacency) is an N x N matrix describing a power system with N buses. It represents the directed connectivity of the buses in a power system.

The Adjacency Struct is indexed using the Bus Numbers, no need for them to be sequential
"""
struct Adjacency{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Int}
    data::SparseArrays.SparseMatrixCSC{Int, Int}
    axes::Ax
    lookup::L
end

"""
Builds a Adjacency from a collection of buses and branches. The return is an N x N Adjacency Array indexed with the bus numbers.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
"""
function Adjacency(branches, nodes; check_connectivity::Bool = true, kwargs...)
    buscount = length(nodes)
    bus_ax = get_number.(nodes)
    axes = (bus_ax, bus_ax)
    bus_lookup = _make_ax_ref(bus_ax)
    look_up = (bus_lookup, bus_lookup)

    a = SparseArrays.spzeros(Int, buscount, buscount)

    for b in branches
        (fr_b, to_b) = get_bus_indices(b, bus_lookup)
        a[fr_b, to_b] = 1
        a[to_b, fr_b] = -1
        a[fr_b, fr_b] = 1
        a[to_b, to_b] = 1
    end

    if check_connectivity
        connected = validate_connectivity(a, nodes, bus_lookup; kwargs...)
        !connected && throw(DataFormatError("Network not connected"))
    end

    return Adjacency(a, axes, look_up)
end

"""
Builds a Adjacency from the system. The return is an N x N Adjacency Array indexed with the bus numbers.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
- `connectivity_method::Function = goderya_connectivity`: method (`goderya_connectivity` or `dfs_connectivity`) for connectivity validation
"""
function Adjacency(sys::System; check_connectivity::Bool = true, kwargs...)
    nodes = sort!(
        collect(get_components(Bus, sys, x -> get_bustype(x) != BusTypes.ISOLATED)),
        by = x -> get_number(x),
    )
    branches = get_components(Branch, sys, get_available)
    return Adjacency(branches, nodes; check_connectivity = check_connectivity, kwargs...)
end

function _goderya(ybus::SparseArrays.SparseMatrixCSC)
    node_count = size(ybus)[1]
    max_I = node_count^2
    I, J, val = SparseArrays.findnz(ybus)
    T = SparseArrays.sparse(I, J, ones(Int, length(val)))
    T_ = T * T
    for n in 1:(node_count - 1)
        I, _, _ = SparseArrays.findnz(T_)
        if length(I) == max_I
            @info "The System has no islands"
            break
        elseif length(I) < max_I
            temp = T_ * T
            I_temp, _, _ = SparseArrays.findnz(temp)
            if all(I_temp == I)
                @warn "The system contains islands" maxlog = 1
            end
            T_ = temp
        else
            @assert false
        end
        #@assert n < node_count - 1
    end
    return I
end

"""
Checks the network connectivity of the system.

# Keyword arguments
- `connectivity_method::Function = goderya_connectivity`: Specifies the method used as Goderya's algorithm (`goderya_connectivity`) or depth first search/network traversal (`dfs_connectivity`)
* Note that the default Goderya method is more efficient, but is resource intensive and may not scale well on large networks.
"""
function validate_connectivity(
    sys::System;
    connectivity_method::Function = goderya_connectivity,
)
    nodes = sort!(
        collect(get_components(Bus, sys, x -> get_bustype(x) != BusTypes.ISOLATED)),
        by = x -> get_number(x),
    )
    branches = get_components(Branch, sys, get_available)
    a = Adjacency(branches, nodes; check_connectivity = false)

    return validate_connectivity(
        a.data,
        nodes,
        a.lookup[1];
        connectivity_method = connectivity_method,
    )
end

function validate_connectivity(
    M,
    nodes::Vector{Bus},
    bus_lookup::Dict{Int64, Int64};
    connectivity_method::Function = goderya_connectivity,
)
    connected = connectivity_method(M, nodes, bus_lookup)
    return connected
end

function goderya_connectivity(M, nodes::Vector{Bus}, bus_lookup::Dict{Int64, Int64})
    @info "Validating connectivity with Goderya algorithm"
    length(nodes) > 15_000 &&
        @warn "The Goderya algorithm is memory intensive on large networks and may not scale well, try `connectivity_method = PowerSystems.dfs_connectivity"

    I = _goderya(M)

    node_count = length(nodes)
    connections = Dict([i => count(x -> x == i, I) for i in Set(I)])

    if length(Set(I)) == node_count
        connected = true
        if any(values(connections) .!= node_count)
            cc = Set(values(connections))
            @warn "Network has at least $(length(cc)) connected components with $cc nodes"
            connected = false
        end
    else
        disconnected_nodes = get_name.(nodes[setdiff(values(bus_lookup), I)])
        @warn "Principal connected component does not contain:" disconnected_nodes
        connected = false
    end
    return connected
end

"""
Finds the set of bus numbers that belong to each connnected component in the System
"""
# this function extends the PowerModels.jl implementation to accept a System
function find_connected_components(sys::System)
    a = Adjacency(sys; check_connectivity = false)
    return find_connected_components(a.data, a.lookup[1])
end

# this function extends the PowerModels.jl implementation to accept an adjacency matrix and bus lookup
function find_connected_components(M, bus_lookup::Dict{Int64, Int64})
    pm_buses = Dict([i => Dict("bus_type" => 1, "bus_i" => b) for (i, b) in bus_lookup])

    arcs = findall((LinearAlgebra.UpperTriangular(M) - LinearAlgebra.I) .!= 0)
    pm_branches = Dict([
        i => Dict("f_bus" => a[1], "t_bus" => a[2], "br_status" => 1) for
        (i, a) in enumerate(arcs)
    ],)

    data = Dict("bus" => pm_buses, "branch" => pm_branches)
    cc = calc_connected_components(data)
    bus_decode = Dict(value => key for (key, value) in bus_lookup)
    connected_components = Vector{Set{Int64}}()
    for c in cc
        push!(connected_components, Set([bus_decode[b] for b in c]))
    end
    return Set(connected_components)
end

function dfs_connectivity(M, _::Vector{Bus}, bus_lookup::Dict{Int64, Int64})
    @info "Validating connectivity with depth first search (network traversal)"
    cc = find_connected_components(M, bus_lookup)
    if length(cc) != 1
        @warn "Network has at least $(length(cc)) connected components with $(length.(cc)) nodes"
        connected = false
    else
        connected = true
    end
    return connected
end


text:
mutable struct QM_FloatData_ref{T<:Real} <: Abstract_QM_FloatData{T}
    Q                   :: SparseMatrixCSC{T,Int}
    AT                  :: SparseMatrixCSC{T,Int} # using AT is easier to form systems
    b                   :: Vector{T}
    c                   :: Vector{T}
    c0                  :: T
    lvar                :: Vector{T}
    uvar                :: Vector{T}
    Δref                :: T
    x_approx            :: Vector{T}
    x_approxQx_approx_2 :: T
    b_init              :: Vector{T}
    bTy_approx          :: T
    c_init              :: Vector{T}
    pri_obj_approx      :: T
end

function fd_refinement(fd :: QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T}, Δxy :: Vector{T}, pt :: Point{T},
                       itd :: IterData{T}, ϵ :: Tolerances{T}, dda :: DescentDirectionAllocs{T}, pad :: PreallocatedData{T},
                       cnts :: Counters, T0 :: DataType, refinement :: Symbol; centering :: Bool = false) where {T<:Real}

    # center Points before zoom
    if centering
        # if pad.fact_fail # re-factorize if it failed in a lower precision system
        #     out = factorize_K2!(pad.K, pad.K_fact, pad.D, pad.diag_Q, pad.diagind_K, pad.regu,
        #                         pt.s_l, pt.s_u, itd.x_m_lvar, itd.uvar_m_x, id.ilow, id.iupp,
        #                         id.ncon, id.nvar, cnts, itd.qp, T, T0)
        # end
        # dda.rxs_l .= -itd.μ
        # dda.rxs_u .= itd.μ
        # solve_augmented_system_cc!(pad.K_fact, itd.Δxy, itd.Δs_l, itd.Δs_u, itd.x_m_lvar, itd.uvar_m_x,
        #                         pad.rxs_l, pad.rxs_u, pt.s_l, pt.s_u, id.ilow, id.iupp)
        itd.Δxy .= 0
        itd.Δxy[id.ilow] .-= itd.μ  ./ itd.x_m_lvar
        itd.Δxy[id.iupp] .+= itd.μ ./ itd.uvar_m_x
        if pad.fact_fail
            out = update_pad!(pad, dda, pt, itd, fd, id, res, cnts, T0)
            out = solver!(pad, dda, pt, itd, fd, id, res, cnts, T0, :IPF)
        else
            out = solver!(pad, dda, pt, itd, fd, id, res, cnts, T0, :cc)
        end
        itd.Δs_l .= @views .-(.-itd.μ .+ pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar
        itd.Δs_u .= @views (itd.μ .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x

        α_pri, α_dual = compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)
        update_data!(pt, α_pri, α_dual, itd, pad, res, fd, id)
    end

    c_ref = fd.c .- itd.ATy .+ itd.Qx
    if refinement == :zoom || refinement == :multizoom
        # zoom parameter
        Δref = one(T) / res.rbNorm
    elseif refinement == :ref || refinement == :multiref
        Δref = one(T)
        αref = T(1.0e12)
        δd = norm(c_ref, Inf)
        if id.nlow == 0 && id.nupp > 0
            δp  = max(res.rbNorm, maximum(itd.uvar_m_x))
        elseif id.nlow > 0 && id.nupp == 0
            δp  = max(res.rbNorm, maximum(itd.x_m_lvar))
        elseif id.nlow == 0 && id.nupp == 0
            δp  = max(res.rbNorm)
        else
            δp  = max(res.rbNorm, maximum(itd.x_m_lvar), maximum(itd.uvar_m_x))
        end
        Δref = max(αref / Δref, one(T) / δp, one(T) / δd)
    end

    c_ref .*= Δref
    fd_ref = QM_FloatData_ref(fd.Q, fd.AT, .-res.rb .* Δref, c_ref, fd.c0,
                             Δref .* (fd.lvar .- pt.x), Δref .* (fd.uvar .- pt.x), Δref, pt.x, copy(itd.xTQx_2), fd.b, dot(fd.b, pt.y),
                             fd.c, copy(itd.pri_obj))

    # init zoom Point
    pt_z = Point(zeros(T, id.nvar), zeros(T, id.ncon), max.(abs.(c_ref[id.ilow]), eps(T)), max.(abs.(c_ref[id.iupp]), eps(T)))
    starting_points!(pt_z, fd_ref, id, itd)

    # update residuals
    res.rb .= itd.Ax .- fd_ref.b
    res.rc .= itd.ATy .- itd.Qx .- fd_ref.c
    res.rc[id.ilow] .+= pt_z.s_l
    res.rc[id.iupp] .-= pt_z.s_u
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)

    return fd_ref, pt_z
end

function update_pt_ref!(Δref :: T, pt :: Point{T}, pt_z :: Point{T}, res :: Residuals{T}, id :: QM_IntData,
                        fd :: QM_FloatData{T}, itd :: IterData{T}) where {T<:Real}

    # update Point
    pt.x .+= pt_z.x ./ Δref
    pt.y .+= pt_z.y ./ Δref
    pt.s_l .= pt_z.s_l ./ Δref
    pt.s_u .= pt_z.s_u ./ Δref

    # update IterData
    itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, :U), pt.x)
    itd.xTQx_2 = dot(pt.x, itd.Qx) / 2
    itd.ATy = mul!(itd.ATy, fd.AT, pt.y)
    itd.Ax = mul!(itd.Ax, fd.AT', pt.x)
    itd.cTx = dot(fd.c, pt.x)
    itd.pri_obj = itd.xTQx_2 + itd.cTx + fd.c0
    itd.dual_obj = @views dot(fd.b, pt.y) - itd.xTQx_2 + dot(pt.s_l, fd.lvar[id.ilow]) - dot(pt.s_u, fd.uvar[id.iupp]) + fd.c0
    itd.pdd = abs(itd.pri_obj - itd.dual_obj ) / (one(T) + abs(itd.pri_obj))

    # update Residuals
    res.rb .= itd.Ax .- fd.b
    res.rc .= itd.ATy .- itd.Qx .- fd.c
    res.rc[id.ilow] .+= pt.s_l
    res.rc[id.iupp] .-= pt.s_u
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
end

function update_data!(pt :: Point{T}, α_pri :: T, α_dual :: T, itd :: IterData{T}, pad :: PreallocatedData{T},
                      res :: Residuals{T}, fd :: QM_FloatData_ref{T}, id :: QM_IntData) where {T<:Real}

    # (x, y, s_l, s_u) += α * Δ
    update_pt!(pt.x, pt.y, pt.s_l, pt.s_u, α_pri, α_dual, itd.Δxy, itd.Δs_l, itd.Δs_u, id.ncon, id.nvar)

    # update IterData
    itd.x_m_lvar .= @views pt.x[id.ilow] .- fd.lvar[id.ilow]
    itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt.x[id.iupp]
    boundary_safety!(itd.x_m_lvar, itd.uvar_m_x, id.nlow, id.nupp, T)

    itd.μ = compute_μ(itd.x_m_lvar, itd.uvar_m_x, pt.s_l, pt.s_u, id.nlow, id.nupp)
    itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, :U), pt.x)
    x_approxQx = dot(fd.x_approx, itd.Qx)
    itd.xTQx_2 =  dot(pt.x, itd.Qx) / 2
    itd.ATy = mul!(itd.ATy, fd.AT, pt.y)
    itd.Ax = mul!(itd.Ax, fd.AT', pt.x)
    itd.cTx = dot(fd.c_init, pt.x)
    itd.pri_obj = fd.pri_obj_approx + x_approxQx / fd.Δref + itd.xTQx_2 / fd.Δref^2 + dot(fd.c_init, pt.x) / fd.Δref
    itd.dual_obj = @views fd.bTy_approx + dot(fd.b_init, pt.y) / fd.Δref - fd.x_approxQx_approx_2 - dot(fd.x_approx, itd.Qx) / fd.Δref -
                        itd.xTQx_2 / fd.Δref^2 +
                        dot(pt.s_l, fd.lvar[id.ilow]) / fd.Δref^2 + dot(pt.s_l, fd.x_approx[id.ilow]) / fd.Δref -
                        dot(pt.s_u, fd.uvar[id.iupp]) / fd.Δref^2 - dot(pt.s_u, fd.x_approx[id.iupp]) / fd.Δref + fd.c0
    itd.pdd = abs(itd.pri_obj - itd.dual_obj ) / (one(T) + abs(itd.pri_obj))

    #update Residuals
    res.n_Δx = @views α_pri * norm(itd.Δxy[1:id.nvar])
    res.rb .= itd.Ax .- fd.b
    res.rc .= itd.ATy .- itd.Qx .- fd.c
    res.rc[id.ilow] .+= pt.s_l
    res.rc[id.iupp] .-= pt.s_u
    res.rcNorm, res.rbNorm = norm(res.rc, Inf) / fd.Δref, norm(res.rb, Inf) / fd.Δref
end


text:
include("quality.jl")
include("models.jl")

using Random
rvec = rand(10)
println( quality.rankdata( rvec ) )
println( quality.fractional_ranks( rvec ) )
println( quality.qvalues(models.pvalue_list) )

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
#=##############################################################################################
Filename: set_wing_reference.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: define an `Action` object to solve a CCBlade rotor
=###############################################################################################

"""
    set_wing_reference(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

Action function.

# Arguments:

* `aircraft::Aircraft`: `Aircraft` system object
* `parameters<:Parameters`: `Parameters` struct
* `freestream::Freestream`: `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray`: array of times for which the simulation is run
* `stepi::Int`: index of the current step
* `step_symbol::String`: defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step

"""
function set_wing_reference(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    # interpret freestream
    vlmfreestream = VL.Freestream(freestream)
    vwake = parameters.wake_function[stepi]

    # read current reference
    reference = aircraft.wing_system.system.reference[1]
    S = reference.S
    c = reference.c
    b = reference.b
    r = reference.r
    V = reference.V

    # get new V
    newV = freestream.vinf

    # build new reference
    newreference = VL.Reference(S,c,b,r,newV)

    # set
    aircraft.wing_system.system.reference[1] = newreference

    return false
end

"""
    set_wing_reference(aircraft, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft` : aircraft system to be simulated
* `step_range::AbstractArray` : defines each step of the simulation

# Returns:

* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step

"""
function set_wing_reference(aircraft, step_range)

    wake_function = Function[(x) -> [0.0, 0.0, 0.0] for i in 1:length(step_range)]

    return wake_function
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using RegularizedLeastSquares, LinearAlgebra
# Packages for tetsing only
using Random, Test

include("Kaczmarz.jl")
include("testProxMaps.jl")
include("testSolvers.jl")


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x1'(t) = -(a21 + a31 + a01) * x1(t) + a12 * x2(t) + a13 * x3(t) + u(t),
  x2'(t) = a21 * x1(t) - a12 * x2(t),
  x3'(t) = a31 * x1(t) - a13 * x3(t),
  y(t) = x1(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using BinDeps
using CxxWrap

@BinDeps.setup

build_type  = "Release"

jlcxx_cmake_dir   = joinpath(CxxWrap.prefix_path(), "lib", "cmake", "JlCxx")
xtl_cmake_dir     = joinpath(dirname(@__FILE__), "usr", "lib", "cmake", "xtl")
xtensor_cmake_dir = joinpath(dirname(@__FILE__), "usr", "lib", "cmake", "xtensor")
julia_bindir      = Sys.BINDIR

prefix                    = joinpath(dirname(@__FILE__), "usr")
xtl_srcdir                = joinpath(dirname(@__FILE__), "xtl")
xtensor_srcdir            = joinpath(dirname(@__FILE__), "xtensor")
xtensor_julia_srcdir      = joinpath(dirname(@__FILE__), "xtensor-julia")
xtensor_examples_srcdir   = joinpath(dirname(@__FILE__), "xtensor-julia-examples")

xtl_builddir              = joinpath(dirname(@__FILE__), "..", "builds", "xtl")
xtensor_builddir          = joinpath(dirname(@__FILE__), "..", "builds", "xtensor")
xtensor_julia_builddir    = joinpath(dirname(@__FILE__), "..", "builds", "xtensor-julia")
xtensor_examples_builddir = joinpath(dirname(@__FILE__), "..", "builds", "xtensor-julia-examples")

# Set generator if on windows
@static if Sys.iswindows()
    genopt = "MinGW Makefiles"
else
    genopt = "Unix Makefiles"
end

# Build on windows: push BuildProcess into BinDeps defaults
@static if Sys.iswindows()
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    saved_defaults = deepcopy(BinDeps.defaults)
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, [BuildProcess])
  end
end

# Functions library for testing
example_labels = [:tensors]
xtensor_examples = BinDeps.LibraryDependency[]
for l in example_labels
   @eval $l = $(library_dependency(string(l), aliases=["lib" * string(l)]))
   push!(xtensor_examples, eval(:($l)))
end

# Version of xtl to vendor
xtl_version = "0.7.4"

# Version of xtensor to vendor
xtensor_version = "0.24.0"

# Version of xtensor-julia to vendor
xtensor_julia_version = "0.10.1"

xtl_steps = @build_steps begin
  `git clone -b $xtl_version --single-branch https://github.com/QuantStack/xtl $xtl_srcdir`
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DBUILD_TESTS=OFF -DCMAKE_INSTALL_LIBDIR=lib $xtl_srcdir`
  `cmake --build . --config $build_type --target install`
end

xtensor_steps = @build_steps begin
  `git clone -b $xtensor_version --single-branch https://github.com/QuantStack/xtensor $xtensor_srcdir`
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DBUILD_TESTS=OFF -DCMAKE_INSTALL_LIBDIR=lib $xtensor_srcdir`
  `cmake --build . --config $build_type --target install`
end

xtensor_julia_steps = @build_steps begin
  `git clone -b $xtensor_julia_version --single-branch https://github.com/QuantStack/xtensor-julia $xtensor_julia_srcdir`
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DJlCxx_DIR=$jlcxx_cmake_dir -Dxtensor_DIR=$xtensor_cmake_dir -DCMAKE_PROGRAM_PATH=$julia_bindir -DCMAKE_INSTALL_LIBDIR=lib $xtensor_julia_srcdir`
  `cmake --build . --config $build_type --target install`
end

xtensor_examples_steps = @build_steps begin
  `cmake -G "$genopt" -DCMAKE_PREFIX_PATH=$prefix -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_BUILD_TYPE="$build_type" -DJlCxx_DIR=$jlcxx_cmake_dir -Dxtensor_DIR=$xtensor_cmake_dir -DCMAKE_PROGRAM_PATH=$julia_bindir -DCMAKE_INSTALL_LIBDIR=lib $xtensor_examples_srcdir`
  `cmake --build . --config $build_type --target install`
end

provides(BuildProcess,
  (@build_steps begin

    println("Building xtl")
    CreateDirectory(xtl_builddir)
    @build_steps begin
      ChangeDirectory(xtl_builddir)
      xtl_steps
    end

    println("Building xtensor")
    CreateDirectory(xtensor_builddir)
    @build_steps begin
      ChangeDirectory(xtensor_builddir)
      xtensor_steps
    end

    println("Building xtensor-julia")
    CreateDirectory(xtensor_julia_builddir)
    @build_steps begin
      ChangeDirectory(xtensor_julia_builddir)
      xtensor_julia_steps
    end

    println("Building xtensor-julia-examples")
    CreateDirectory(xtensor_examples_builddir)
    @build_steps begin
      ChangeDirectory(xtensor_examples_builddir)
      xtensor_examples_steps
    end
  end), xtensor_examples)

@BinDeps.install Dict([
    (:tensors, :_l_tensors)
])

# Build on windows: pop BuildProcess from BinDeps defaults
@static if Sys.iswindows()
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, saved_defaults)
  end
end


text:
module ComplexWrapper

using JavaCall

const CalculusFieldElement = @jimport org.hipparchus.CalculusFieldElement
const Class = @jimport java.lang.Class
const Comparator = @jimport java.util.Comparator
const ComplexComparator = @jimport org.hipparchus.complex.ComplexComparator
const ComplexField = @jimport org.hipparchus.complex.ComplexField
const ComplexFormat = @jimport org.hipparchus.complex.ComplexFormat
const ComplexUtils = @jimport org.hipparchus.complex.ComplexUtils
const Double = @jimport java.lang.Double
const FieldPosition = @jimport java.text.FieldPosition
const FieldSinCos = @jimport org.hipparchus.util.FieldSinCos
const Function = @jimport java.util.function.Function
const JComplex = @jimport org.hipparchus.complex.Complex
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Locale = @jimport java.util.Locale
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const ParsePosition = @jimport java.text.ParsePosition
const Quaternion = @jimport org.hipparchus.complex.Quaternion
const RootsOfUnity = @jimport org.hipparchus.complex.RootsOfUnity
const StringBuffer = @jimport java.lang.StringBuffer
const ToDoubleFunction = @jimport java.util.function.ToDoubleFunction
const ToIntFunction = @jimport java.util.function.ToIntFunction
const ToLongFunction = @jimport java.util.function.ToLongFunction

export ComplexComparator
export ComplexField
export ComplexFormat
export ComplexUtils
export JComplex
export Quaternion
export RootsOfUnity
export abs
export acos
export acosh
export add
export asin
export asinh
export atan
export atan2
export atanh
export cbrt
export ceil
export compare
export compute_roots
export conjugate
export convert_to_complex
export copy_sign
export cos
export cosh
export divide
export dot_product
export equals
export equals_with_relative_tolerance
export exp
export expm1
export floor
export format
export get_argument
export get_available_locales
export get_class
export get_complex_format
export get_conjugate
export get_exponent
export get_field
export get_imaginary
export get_imaginary_character
export get_imaginary_format
export get_instance
export get_inverse
export get_norm
export get_number_of_roots
export get_one
export get_positive_polar_form
export get_q0
export get_q1
export get_q2
export get_q3
export get_real
export get_real_format
export get_runtime_class
export get_scalar_part
export get_vector_part
export get_zero
export hash_code
export hypot
export is_counter_clock_wise
export is_infinite
export is_mathematical_integer
export is_na_n
export is_pure_quaternion
export is_real
export is_unit_quaternion
export linear_combination
export log
export log10
export log1p
export multiply
export negate
export new_instance
export norm
export normalize
export notify
export notify_all
export nth_root
export parse
export polar2_complex
export pow
export reciprocal
export remainder
export reversed
export rint
export root_n
export scalb
export signum
export sin
export sin_cos
export sinh
export sqrt
export sqrt1z
export subtract
export tan
export tanh
export then_comparing
export then_comparing_double
export then_comparing_int
export then_comparing_long
export to_degrees
export to_radians
export to_string
export value_of
export wait

include("complex_comparator.jl")
include("complex_field.jl")
include("complex_format.jl")
include("complex_utils.jl")
include("j_complex.jl")
include("quaternion.jl")
include("roots_of_unity.jl")

end


text:
@testset "QuasiSteady + Section" begin

    coupling = Coupling(QuasiSteady(), Section())

    # test provided jacobians
    coupling_jacobian_tests(coupling)

    # test consistency of input/output functions
    coupling_io_tests(coupling)
end


text:
module Presentations

  using ..DB
  import ..DB: @db_schema

  using Catlab
  using Catlab.Theories
  using Catlab.Graphics
  using Catlab.Present: Presentation
  using Catlab.Theories.FreeSchema: Attr
  using AlgebraicPetri
  import Catlab.Theories: FreeSymmetricMonoidalCategory, ⊗
  import Catlab.Programs: @program

  export present_to_schema, @program, draw_workflow, FreeSymmetricMonoidalCategory,
         add_types!, add_type!, add_process!, add_processes!, Presentation, ⊗,
         draw_schema, @db_schema, @present_to_schema

  hasprop(o, p) = p in propertynames(o)

  Presentation() = Presentation(FreeSymmetricMonoidalCategory)

  function add_type!(p::Presentation, field::Tuple{Symbol, <:Type})
    ob = Ob(FreeSymmetricMonoidalCategory, field[1])
    push!(ob.args, Symbol(field[2]))
    add_generator!(p, ob)
    return ob
  end

  function add_types!(p::Presentation, fields::Array{<:Tuple{Symbol, Type}})
    return map(field->add_type!(p, field), fields)
  end

  function add_process!(p::Presentation, hom::Tuple{Symbol, GATExpr, GATExpr})
    h = Hom(hom...)
    add_generator!(p, h)
    return h
  end

  function add_processes!(p::Presentation, homs::Array{<:Tuple{Symbol, <:GATExpr, <:GATExpr},1})
    return map(hom->add_process!(p, hom), homs)
  end

  macro present_to_schema(head)
    present_to_schema(head.args[1], head.args[2])
  end

  function present_to_schema(sch_name::Symbol, wf::Symbol)
    quote
      gens = Array{GATExpr, 1}()
      tables = Dict{Symbol, GATExpr}()
      sym_app(s::Symbol, suffix::String) = Symbol(string(s, suffix))
      get_syms(g::Union{GATExpr, Array}) = begin
        if g isa Array
          if eltype(g) == Symbol
            return [g]
          end
          sym_array = Array{Array{Symbol, 1}, 1}()
          for i in g
            append!(sym_array, get_syms(i))
          end
          return sym_array
        elseif hasprop(g, :args)
          return get_syms(g.args)
        end
      end

      # Evaluate homs to purely data-connected tables
      for g in generators($(esc(wf)), :Hom)
        g_name = g.args[1]
        table = Ob(FreeSchema, g_name)
        tables[g_name] = table
        push!(gens, table)
        append!(gens, map(enumerate(get_syms(g.type_args))) do (i, sym)
                  Attr(sym_app(g_name, "_$(i)_$(sym[1])$i"),
                       table, generator(TheorySQL, sym[2]))
        end)
      end

      @present p <: TheorySQL begin end
      add_generators!(p, gens)
      $(esc(:eval))(:(@db_schema $($(Meta.quot(sch_name)))($(p))))
    end
  end

  function draw_schema(p::Presentation; kw...)
    ob_names = Symbol.(generators(p, :Ob))
    hom_names = Symbol.(generators(p, :Hom))
    hom_dict = Dict{Symbol, Tuple}()
    for hom in generators(p, :Hom)
      # Evaluate Dom
      dom_v = if eltype(dom(hom).args) <: GATExpr
        Tuple(Symbol.(dom(hom).args))
      else
        Symbol(dom(hom))
      end
      # Operate on Codom
      codom_v = if eltype(codom(hom).args) <: GATExpr
        Tuple(Symbol.(codom(hom).args))
      else
        Symbol(codom(hom))
      end

      hom_dict[Symbol(hom)] = (dom_v, codom_v)
    end
    schema_ap = LabelledPetriNet(ob_names, hom_dict...)
    Graph(schema_ap)
  end

  function draw_workflow(p; kw...)
    to_graphviz(p; orientation=LeftToRight, kw...)
  end
end


text:
@doc """
	density(pts::IndexablePoints2D) -> Real

Here density is defined as the minimum radius of a covering of
the GeometricalPredicates region with equal-sized balls centered at the
points in `pts`.

By definition of the Voronoi tesselation this radius is the maximum
distance from a Voronoi cell vertix to its generator.
"""->
function density(generators::IndexablePoints2D)
	corn = voronoicells(generators)
	density(generators, corn)
end

function density(generators::IndexablePoints2D, corners::Tessellation)
	sort!(generators, by=getindex)
	dens = 0.0

	for idx in keys(corners)
		gen = generators[idx]
		for C in corners[idx]
			dens = max(dens, dist_squared(gen,C))
		end
	end

	return sqrt(dens)
end

@doc """
	density(x::Vector, y::Vector; rw) -> Real

Compute the density for points with coordinates `x` and `y` in the window `rw`.

The vector `rw` specifies the boundary rectangle as `[xmin, xmax, ymin, ymax]`.
By default, `rw` is the unit rectangle.
"""->
function density{T<:Real}(x::AbstractVector{T}, y::AbstractVector{T}, rw::Vector{Float64}=[0.0;1.0;0.0;1.0])
	pts, SCALEX, SCALEY = fit2boundingbox(x, y, rw)

	density(pts) * sqrt(SCALEX*SCALEY)
end


text:
# load the examples array
include("hclust-generated-examples.jl")

# test to make sure many random examples match R's implementation
for example in examples
    h = hclust(example["D"], example["method"])
    for i in 1:size(h.merge)[2]
        for j in 1:size(h.merge)[1]
            @test h.merge[j,i] == example["merge"][j,i]
        end
    end
    for i in 1:length(h.height)
        @test_approx_eq_eps h.height[i] example["height"][i] 1e-5
    end
    for i in 1:length(h.order)
        @test h.order[i] == example["order"][i]
    end
end


text:
"""
    j_integral(f, λ)

Compute the integral

```math
j(f, \\lambda) = \\lambda \\frac{\\int_0^{T_\\text{echo}} F(t) \\int_0^t
\\mathrm{e}^{-\\lambda(t - s)} f(s) \\, \\mathrm{d} s \\, \\mathrm{d}
t}{\\int_0^{T_\\text{echo}} F^2(t) \\, \\mathrm{d} t}
```

for the time profile `f` and Laplace eigenvalue `λ`.
"""
function j_integral end

function j_integral(f::PGSE, λ)
    (; δ, Δ) = f

    if λ < 1e-7
        # Use Taylor expansion when λ is close to 0
        # to improve numerical stability
        λ - λ^2 * Δ^2 / 2(Δ - δ / 3) + λ^3 * (10Δ^3 + 5Δ * δ^2 - δ^3) / 20(3 * Δ - δ) -
        λ^4 * (Δ^4 + Δ^2 * δ^2) / 8(3Δ - δ) +
        λ^5 * (21Δ^5 + 35Δ^3 * δ^2 + 7Δ * δ^4 - δ^5) / 840(3Δ - δ)
    else
        -(
            exp(-λ * (Δ + δ)) + exp(-λ * (Δ - δ)) - 2exp(-λ * δ) - 2exp(-λ * Δ) +
            2(1 - λ * δ)
        ) / (λ^2 * δ^2 * (Δ - δ / 3))
    end
end

function j_integral(f::DoublePGSE, λ)
    #J Compute the quantity J(λ) for the sequence
    #   An analytical expression is available for the DoublePGSE sequence
    (; δ, Δ, tpause) = f
    tm = tpause + δ

    if λ < 1e-7
        # Use Taylor expansion when λ is close to 0
        # to improve numerical stability
        λ - λ^2 * 3Δ^2 / (3Δ - δ) +
        λ^3 * (40Δ^3 + 5Δ * δ^2 - δ^3 + 30Δ^2 * tm) / 20(3Δ - δ) -
        λ^4 * (4Δ^4 + Δ^2 * δ^2 + 6Δ^3 * tm + 3Δ^2 * tm^2) / 4(3Δ - δ) +
        λ^5 * (
            336Δ^5 + 140Δ^3 * δ^2 + 7Δ * δ^4 - δ^5 +
            735Δ^4 * tm +
            105Δ^2 * δ^2 * tm +
            630Δ^3 * tm^2 +
            210Δ^2 * tm^3
        ) / 840(3Δ - δ)
    else
        -(
            exp(-λ * (2Δ + tm - δ)) +
            exp(-λ * (2Δ + tm + δ)) +
            exp(-λ * (tm - δ)) +
            exp(-λ * (tm + δ)) - 2exp(-λ * (Δ + tm - δ)) - 2exp(-λ * (Δ + tm + δ)) -
            2exp(-λ * (2Δ + tm)) +
            2exp(-λ * (Δ - δ)) +
            2exp(-λ * (Δ + δ)) - 2exp(-λ * tm) +
            4 +
            4exp(-λ * (Δ + tm)) - 4exp(-λ * Δ) - 4exp(-λ * δ) - 4λ * δ
        ) / (2λ^2 * δ^2 * (Δ - δ / 3))
    end
end

function j_integral(f::CosOGSE, λ)
    #J Compute the quantity J(λ) for the sequence
    #   An analytical expression is available for the CosOGSE sequence
    (; δ, Δ, nperiod) = f
    n = nperiod

    if λ < 1e-7
        # Use Taylor expansion when λ is close to 0
        # to improve numerical stability
        λ - λ^3 * δ^2 * 3 / 4(n * π)^2 +
        λ^5 * (12Δ * n^2 * π^2 * δ^3 + 15δ^4 - 4n^2 * π^2 * δ^4) / (48n^4 * π^4) -
        λ^6 * Δ^2 * δ^3 / 8(n * π)^2
    else
        4n^2 *
        π^2 *
        λ *
        (
            -exp(-λ * (Δ + δ)) * λ * δ - exp(-λ * (Δ - δ)) * λ * δ +
            2exp(-λ * Δ) * λ * δ +
            2exp(-λ * δ) * λ * δ +
            4n^2 * π^2 +
            (λ * δ - 2) * λ * δ
        ) / (4n^2 * π^2 + λ^2 * δ^2)^2
    end
end

function j_integral(f::SinOGSE, λ)
    #J Compute the quantity J(λ) for the sequence
    #   An analytical expression is available for the SinOGSE sequence
    (; δ, Δ, nperiod) = f
    n = nperiod

    if λ < 1e-7
        # Use Taylor expansion when λ is close to 0
        # to improve numerical stability
        λ +
        λ^3 * δ * (4n^2 * π^2 * δ - 12Δ * n^2 * π^2 - 15δ) / (36n^2 * π^2) +
        λ^4 * δ * Δ^2 / 6 +
        λ^5 *
        δ *
        (
            120Δ * n^2 * π^2 * δ^2 + 4n^4 * π^4 * δ^3 - 20Δ * n^4 * π^4 * δ^2 + 105δ^3 -
            40n^2 * π^2 * δ^3 - 40Δ^3 * n^4 * π^4
        ) / (720n^4 * π^4) +
        λ^6 * ((Δ^4 * δ^2 + Δ^2 * δ^4) / 24 - (Δ^2 * δ^4) / (4n^2 * π^2)) / 3δ
    else
        4π^2 *
        n^2 *
        (
            (λ * δ)^3 +
            4n^2 *
            π^2 *
            (
                exp(-λ * (Δ + δ)) + exp(-λ * (Δ - δ)) - 2exp(-λ * δ) - 2exp(-λ * Δ) +
                λ * δ +
                2
            )
        ) / (3δ * (4n^2 * π^2 + λ^2 * δ^2)^2)
    end
end


text:
using CLI
using Base.Test

@testset "CLI Elements" begin

@testset "basic flags" begin

# NOTE: isoption only captures prefix
@test CLI.isoption("-h") == true
@test CLI.isoption("--help") == true
@test CLI.isoption("cmd") == false
@test CLI.isoption("@file") == false
@test CLI.isoption(nothing) == false

@test CLI.islong("-h") == false
@test CLI.islong("--help") == true
@test CLI.islong("cmd") == false
@test CLI.islong("@file") == false
@test CLI.islong(nothing) == false

end # basic flags

@testset "Option Compare" begin
    testop = Option("help", :help, "No Doc")
    @test "--help" == testop
    @test "-h" != testop

    testop = ShortOption('h', "help", :help, "No Doc")
    @test "--help" == testop
    @test "-h" == testop
end # Option Compare

# @testset "Option Search" begin
#     options = AbstractOption[
#         Option("foo", :foo, "No Doc"),
#         ShortOption('c', "cool", :cool, "No Doc"),
#         Option("aaa-aaa", :aaa_aaa, "No Doc")
#     ]

#     @test CLI.findoption("--foo", options) == options[1]
#     @test CLI.findoption("-f", options) == nothing

#     @test CLI.findoption("-c", options) == options[2]
#     @test CLI.findoption("-cool", options) == nothing
#     @test CLI.findoption("--cool", options) == options[2]

#     @test CLI.findoption("aaa-aaa", options) == options[3]
#     @test CLI.findoption("--aaa-aaa", options) == options[3]
#     @test CLI.findoption("aaa_aaa", options) == nothing
#     @test CLI.findoption("--aaa_aaa", options) == nothing

# end # Option Search

end # CLI Elements


text:
using Trapla

Tf = 10.0 # Final time
Δt = 0.1  # Time discretization
tspan = 0.0 : Δt : Tf # Time span

w = 1.5 # Steepness

y = smoothjump.(tspan, Tf; ω = w); # Planned trajectory

# using Plots
# plot(tspan, y)


text:
module CrystalTest
using Crystals
using DataFrames: nrow, DataFrame, ncol, deleterows!, missing, ismissing, eachrow
using Base.Test
using Unitful
using Documenter

contains(x) = y -> x ∈ y
all_integers(x::Array, ε::AbstractFloat=1e-8) = all(abs(x - round(Integer, x)) .< ε)
all_integers(ε::AbstractFloat=1e-8) = y -> all_integers(y, ε)
is_subtype(x::Type) = y -> (y <: x)

@testset "Crystal" begin
    include("Crystal.jl")
end
@testset "Atoms and iteration" begin
    include("CrystalAtoms.jl")
end
@testset "SNF" begin
    include("SNF.jl")
end
@testset "Utilities" begin
    include("utilities.jl")
end
@testset "SpaceGroup" begin
    include("SpaceGroup.jl")
end
@testset "Gruber" begin
    include("Gruber.jl")
end

@testset "Docstests" begin
    mktempdir() do path
        makedocs(modules = [Crystals], clean = true,
             build = joinpath(path, "build"),
             pages = Crystals._doc_pages(),
             root=joinpath(dirname(@__DIR__), "docs"))
    end
end

end


text:
using Base.Test
using PGM

@testset "DMRF" begin

@testset "DPMRF.Basics" begin

    vlst = VarList([dvar("x1", 3),
                    dvar("x2", 3),
                    dvar("x3", 2),
                    dvar("x4", 2),
                    dvar("x5", 3),
                    dvar("x6", 3)])
    elists = [[(1, 2), (5, 6)],
              [(3, 4)],
              [(1, 3), (2, 4), (5, 3), (6, 4)]]
    gr = ugraph_with_tedges(6, elists)
    mrf0 = DiscretePairwiseMRF{Float64}(vlst, gr)

    @test nvars(mrf0) == 6
    @test nvertices(mrf0) == 6
    @test nedgetypes(mrf0) == 3
    @test nedges(mrf0) == 7
    @test [edges(mrf0, t) for t in 1:3] == elists
    @test vars(mrf0) === vlst
    @test [var_cardinalities(mrf0, t) for t in 1:6] == [3, 3, 2, 2, 3, 3]
    @test [edge_pdims(mrf0, t) for t in 1:3] == [(3, 3), (2, 2), (3, 2)]
    @test !has_vertex_potentials(mrf0)
    @test !has_edge_potentials(mrf0)

    vpots = [floor.(rand(d) * 100) for d in [3, 3, 2, 2, 3, 3]]
    epots = [floor.(rand(d1, d2) * 100) for (d1, d2) in [(3, 3), (2, 2), (3, 2)]]

    mrf1 = pwmrf(vlst, gr, vpots, epots)
    @test nvars(mrf1) == 6
    @test nvertices(mrf1) == 6
    @test nedgetypes(mrf1) == 3
    @test nedges(mrf1) == 7
    @test [var_cardinalities(mrf1, t) for t in 1:6] == [3, 3, 2, 2, 3, 3]
    @test [edge_pdims(mrf1, t) for t in 1:3] == [(3, 3), (2, 2), (3, 2)]
    @test has_vertex_potentials(mrf1)
    @test has_edge_potentials(mrf1)
    @test vertex_potentials(mrf1) === vpots
    @test edge_potentials(mrf1) === epots
end

function collect_tpotentials(mrf::DiscretePairwiseMRF)
    r = zeros(mrf.cards...)
    for ci in Base.CartesianRange(size(r))
        r[ci] = tpotential(mrf, collect(ci.I))
    end
    r
end

@testset "DPMRF.Potentials" begin

    vlst = VarList([dvar("x1", 3),
                    dvar("x2", 3),
                    dvar("x3", 2),
                    dvar("x4", 2),
                    dvar("x5", 3),
                    dvar("x6", 3)])
    elists = [[(1, 2), (5, 6)],
              [(3, 4)],
              [(1, 3), (2, 4), (5, 3), (6, 4)]]
    gr = ugraph_with_tedges(6, elists)
    mrf = DiscretePairwiseMRF{Float64}(vlst, gr)

    r0_a = collect_tpotentials(mrf)
    r0_b = tpotentials(mrf)

    z = zeros(3, 3, 2, 2, 3, 3)
    @test r0_a == z
    @test r0_b == z

    vpots = [floor.(rand(d) * 100) for d in [3, 3, 2, 2, 3, 3]]
    epots = [floor.(rand(d1, d2) * 100) for (d1, d2) in [(3, 3), (2, 2), (3, 2)]]
    set_vertex_potentials!(mrf, vpots)
    set_edge_potentials!(mrf, epots)

    vp1, vp2, vp3, vp4, vp5, vp6 = vpots
    ep1, ep2, ep3 = epots

    s1 = [1, 1, 1, 1, 1, 1]
    p1 = vp1[1] + vp2[1] + vp3[1] + vp4[1] + vp5[1] + vp6[1] +
         ep1[1,1] * 2 + ep2[1,1] + ep3[1,1] * 4
    @test tpotential(mrf, s1) == p1

    s2 = [2, 1, 1, 2, 3, 2]
    p2 = vp1[2] + vp2[1] + vp3[1] + vp4[2] + vp5[3] + vp6[2] +
         ep1[2,1] + ep1[3,2] + ep2[1,2] + ep3[2,1] + ep3[1,2] + ep3[3,1] + ep3[2,2]
    @test tpotential(mrf, s2) == p2

    r1_a = collect_tpotentials(mrf)
    r1_b = tpotentials(mrf)
    @test r1_a == r1_b
end

end # DMRF


text:
#= using UnicodePlots =#

# Set up model
function build_model()
    F = diagm([1.0])
    V = diagm([2.0])
    G = reshape([1, 2, -0.5], 3, 1)
    W = diagm([8.0, 2.5, 4.0])
    x0 = randn(1)
    P0 = diagm([1e7])
    StateSpaceModel(F, V, G, W, x0, P0)
end

# Time varying model
function sinusoid_model(omega::Real; fs::Int=256, x0=[0.5, -0.5], W::AbstractFloat=0.1)
    F(n)  = [1.0 0; 0 1.0]
    V(n)  = diagm([1e-10, 1e-10])
    G(n)  = [cos(2*pi*omega*(1/fs)*n) -sin(2*pi*omega*(1/fs)*n)]
    w(n)  = diagm([W])
    P0    = diagm([1e-1, 1e-1])
    StateSpaceModel(F, V, G, w, x0, P0)
end

facts("Kalman Filter") do

    srand(1)

    context("Time invariant models") do

        context("Building model") do
            build_model()
        end

        context("Simulations") do
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
        end

        context("Filtering") do
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
            filt = kalman_filter(y, mod1)
        end

        context("Smoothing") do
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)

            smooth = kalman_filter(y, mod1) |> kalman_smooth
            smooth2 = kalman_smooth(y, mod1)

            @fact smooth.smoothed --> roughly(smooth2.smoothed, atol=1e-2)
        end

        context("Missing data") do
            mod1 = build_model()
            x, y = simulate(mod1, 100)
            y[1:9:end] = NaN
            y[100] = NaN
            filt = kalman_filter(y, mod1)
            smooth = kalman_smooth(filt)
            smooth2 = kalman_smooth(y, mod1)
            @fact any(isnan(filt.filtered)) --> false
            @fact any(isnan(smooth.filtered)) --> false
            @fact any(isnan(smooth2.smoothed)) --> false
        end

        context("Return sizes") do
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
            filt = kalman_filter(y, mod1)
            smooth = kalman_smooth(filt)
            smooth2 = kalman_smooth(y, mod1)

            @fact size(filt.filtered) --> size(filt.predicted)
            @fact size(filt.filtered) --> size(smooth.filtered)
            @fact size(filt.filtered) --> size(smooth.smoothed)
            @fact size(filt.error_cov) --> size(smooth.error_cov)
            @fact size(filt.filtered) --> size(smooth2.smoothed)
            @fact size(filt.error_cov) --> size(smooth2.error_cov)
        end

        context("Linear regression test") do
            m, b, s, dt = 5, 2, 2, .1
            t = 0:dt:10
            y_true = m*t + b
            input = 100*[sin(t/2) sin(t/4) cos(t/2) cos(t/4)] + 10
            y_noisy = [y_true zeros(length(t)) -y_true] +
                        100*[sin(t/2)+sin(t/4) sin(t/2)+cos(t/2) cos(t/2)+cos(t/4)] + 10 + randn(length(t), 3)
            lm = StateSpaceModel([1 dt; 0 1], zeros(2,2),
                                  [1. 0; 0 0; -1 0], s*eye(3),
                                  zeros(2), 100*eye(2),
                                  B=zeros(2, 4), D=[1. 1 0 0; 1 0 1 0; 0 0 1 1])
            lm_filt = kalman_filter(y_noisy, lm, u=input)
            @fact lm_filt.filtered[end,1] --> roughly(y_true[end, 1], atol=3*sqrt(lm_filt.error_cov[1,1,end]))

            lm_smooth = kalman_smooth(lm_filt)
            stderr = sqrt(lm_smooth.error_cov[1,1,:][:])
            @fact lm_filt.filtered[end,:] --> lm_smooth.smoothed[end,:]
            @fact all(abs(y_true - lm_smooth.smoothed[:,1]) .< 3*stderr) --> true
            @fact ones(t) * lm_smooth.smoothed[1,2] --> roughly(lm_smooth.smoothed[:, 2], atol=1e-12)

            # Repeat with DK smoother
            lm_smooth = kalman_smooth(y_noisy, lm, u=input)
            stderr = sqrt(lm_smooth.error_cov[1,1,:][:])
            @fact all(abs(y_true - lm_smooth.smoothed[:,1]) .< 3*stderr) --> true
            @fact ones(t) * lm_smooth.smoothed[1,2] --> roughly(lm_smooth.smoothed[:, 2], atol=1e-12)

        end

    end

    context("Time varying models") do

        context("Building Model") do
            sinusoid_model(40)
        end

        context("Simulations") do
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 256)
            x, y = TimeModels.simulate(mod2, fs*2)

            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(y), width = 120, title = "Original Data")) =#
        end

        context("Filtering") do
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 256)
            x, y = TimeModels.simulate(mod2, fs*2)

            context("Correct initial guess") do
                filt = kalman_filter(y, mod2)
                @fact filt.predicted[end, :] --> roughly([0.5 -0.5]; atol= 0.3)

                #= x_est = round(filt.predicted[end, :], 3) =#
                #= display(lineplot(collect(1:size(x, 1)) / fs, vec(filt.predicted[:, 1]), width = 120, title="Filtered State 1: $(x_est[1])")) =#
                #= display(lineplot(collect(1:size(x, 1)) / fs, vec(filt.predicted[:, 2]), width = 120, title="Filtered State 2: $(x_est[2])")) =#
            end

            context("Incorrect initial guess") do
                mod3 = sinusoid_model(4, fs = 256, x0=[1.7, -0.2])
                filt = kalman_filter(y, mod3)
                @fact filt.predicted[end, :] --> roughly([0.5 -0.5]; atol= 0.3)

            end

            context("Model error") do
                mod4 = sinusoid_model(4, fs = 256, x0=[1.7, -0.2], W=3.0)
                filt = kalman_filter(y, mod4)
                @fact filt.predicted[end, :] --> roughly([0.5 -0.5]; atol= 0.3)
            end

        end

        context("Smoothing") do
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 8192)
            x, y = TimeModels.simulate(mod2, fs*10)
            smooth = kalman_smooth(y, sinusoid_model(4, fs = 8192, x0=[1.7, -0.2]) )
            @fact mean(smooth.smoothed, 1) --> roughly([0.5 -0.5]; atol= 0.1)

            #= x_est = round(smooth.smoothed[end, :], 3) =#
            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(smooth.smoothed[1:end, 1]), width = 120, title="Smoothed State 1: $(x_est[1])")) =#
            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(smooth.smoothed[1:end, 2]), width = 120, title="Smoothed State 2: $(x_est[2])")) =#
        end

    end

    context("Correct failure") do
        @fact_throws StateSpaceModel([1 0.1; 0 1], zeros(2,3), zeros(2,2),
                [1. 0; 0 0; -1 0], [1. 1 0 0; 1 0 1 0; 0 0 1 1], s*eye(3), zeros(2), 100*eye(2))
    end

end