text:
import Robotlib: ad, adi
using Random
"""
    Tn0, xi, et, er = calibLPOE(xin, Tn0in, Ta, q; maxiter=10, λ=1.0)

Performs kinematic calibration using the local POE formulation of kinematics.
- `Tn0` Nominal forward kinematics
- `xi` Twists
- `et` Translational errors as function of iteration of algorithm
- `er` Rotational errors as function of iteration of algorithm
- `λ` Regularization parameter

See `runtests.jl` for simulated usage. Some functions to generate simulation data for calibration purposes are provided in $(@__FILE__)
"""
function calibLPOE(xin,Tn0in,Ta,q;maxiter=10, λ=1.0)
    xi      = deepcopy(xin)
    Tn0     = deepcopy(Tn0in)
    n       = size(xi,2)-1
    N       = size(q,1)
    y       = zeros(6N)
    A       = zeros(6N,6*(n+1))
    xini    = deepcopy(xin)
    Tn0cand = deepcopy(Tn0)
    v       = zeros(3) # only revolute joints
    et      = zeros(maxiter+1)
    er      = zeros(maxiter+1)
    et[1],er[1]   = evalErrorLPOE(xi,Tn0,Ta,q)
    println("Error: ",et[1], " λ: ", λ)
    @assert size(Ta,3) == N

    for iter = 1:maxiter # do a few iterations of the calibration
        for i = 1:N
            Tfull = fkineLPOE(Tn0,xi,q[i,:])
            # populate the matrices of the linear estimation problem
            y[(i-1)*6+1:6i] = twistcoords(logT(Ta[:,:,i]*trinv(Tfull)))
            A[(i-1)*6+1:6*i,1:6] = ad(Tn0[:,:,1])
            adi = ad(Tn0[:,:,1]*expξ(xi[:,1],q[i,1]))
            for j = 2:n
                A[(i-1)*6+1:6*i,(j-1)*6+1:6*j] = adi*ad(Tn0[:,:,j])
                adi = adi*ad(Tn0[:,:,j]*expξ(xi[:,j],q[i,j]))
            end
            A[(i-1)*6+1:6*i,n*6+1:6*(n+1)] = adi*ad(Tn0[:,:,n+1])
        end
        x = (A'A + λ*I)\A'y
        # Update the candidate nominal parameters
        for j = 1:n+1
            ξ = x[(j-1)*6+1:6j]
            Tn0cand[:,:,j] = Tn0[:,:,j]*expξ2(ξ,1)
        end

        et[iter+1],er[iter+1] = evalErrorLPOE(xi,Tn0cand,Ta,q)
        println("Error: ",round(et[iter+1], digits=5), " λ: ", λ, " Norm dx: ", round(norm(x), digits=5))

        if norm(x) < 1e-16
            return Tn0,xi, et, er
        end
        if et[iter+1] > et[iter]
            et[iter+1] = et[iter]
            er[iter+1] = er[iter]
            λ *= 10
        else
            λ/=10
            xi = deepcopy(xini)
            Tn0 = deepcopy(Tn0cand)
        end
    end
    return Tn0, xi, et, er
end

"""
    calibLPOEdual(xi,Tn0,q;maxiter=10, λ=1.0)

Performs dual arm calibration given joint twists xi and nominal transformations Tn0
Returns calibrated nominal transformations.
This function does not converge to the correct kinematic parameters since the problem is underspecified. It only converges to kinematic parameters that make the two arms agree.
"""
function calibLPOEdual(xin,Tn0in,q;maxiter=10, λ=1.0)

    xi      = deepcopy(xin)
    Tn0     = deepcopy(Tn0in)
    n       = size(xi,2)-1
    N       = size(q,1)
    y       = zeros(6N)
    A       = zeros(6N,6*(n+1),2)
    xini    = deepcopy(xin)
    Tn0cand = deepcopy(Tn0)
    et      = zeros(maxiter+1)
    er      = zeros(maxiter+1)
    et[1],er[1]   = evalErrorDual(xi,Tn0,q)
    println("Error: ",et[1], " λ: ", λ)

    for iter = 1:maxiter # do a few iterations of the calibration
        for i = 1:N
            T1 = fkineLPOE(Tn0[:,:,:,1],xi[:,:,1],q[i,:,1])
            T2 = fkineLPOE(Tn0[:,:,:,2],xi[:,:,2],q[i,:,2]) # Use the other arm as measurement
            y[(i-1)*6+1:6i] = twistcoords(logT(T2*trinv(T1)))
            for arm = 1:2
                A[(i-1)*6+1:6*i,1:6,arm] = ad(Tn0[:,:,1,arm])
                adi = ad(Tn0[:,:,1,arm]*expξ(xi[:,1,arm],q[i,1,arm]))
                for j = 2:n
                    A[(i-1)*6+1:6*i,(j-1)*6+1:6*j,arm] = adi*ad(Tn0[:,:,j,arm])
                    adi = adi*ad(Tn0[:,:,j,arm]*expξ(xi[:,j,arm],q[i,j,arm]))
                end
                A[(i-1)*6+1:6*i,n*6+1:6*(n+1),arm] = adi*ad(Tn0[:,:,n+1,arm])

            end
        end
        x = zeros(3)
        for arm = 1:2
            x = (A[:,:,arm]'A[:,:,arm] + λ*I)\A[:,:,arm]'y
            # Update the candidate nominal parameters
            s = arm == 1 ? 1 : -1 # Determine the sign of the correction
            for j = 1:n+1
                ξ = 0.5*s*x[(j-1)*6+1:6j] # Only apply half the update to each arm, verified to be essential
                Tn0cand[:,:,j,arm] = Tn0[:,:,j,arm]*expξ2(ξ,1)
            end

        end
        et[iter+1],er[iter+1] = evalErrorDual(xi,Tn0cand,q)
        println("Error: ",round(et[iter+1], digits=5), " λ: ", λ, " Norm dx: ", round(norm(x), digits=5))

        if norm(x) < 1e-14
            return Tn0,xi, et, er
        end
        if et[iter+1] > et[iter]
            et[iter+1] = et[iter]
            er[iter+1] = er[iter]
            λ *= 10
        else
            λ/=10
            xi = deepcopy(xini)
            Tn0 = deepcopy(Tn0cand)
        end

    end

    return Tn0, xi, et, er

end

function evalErrorLPOE(xin,Tn0,Ta,q)
    N = size(Ta,3)
    et = 0.0
    er = 0.0
    for i = 1:N
        Tfull = fkineLPOE(Tn0,xin,q[i,:])
        et += norm(Ta[1:3,4,i]-Tfull[1:3,4])
        er += norm(skewcoords(logR(Ta[1:3,1:3,i]*Tfull[1:3,1:3]')))
    end
    return et/N, er/N
end

function evalErrorDual(xin,Tn0,q)
    N = size(Tn0,3)
    et = 0.0
    er = 0.0
    for i = 1:N
        T1 = fkineLPOE(Tn0[:,:,:,1],xin[:,:,1],q[i,:,1])
        T2 = fkineLPOE(Tn0[:,:,:,2],xin[:,:,2],q[i,:,2])
        et += norm(T2[1:3,4]-T1[1:3,4])
        er += norm(skewcoords(logR(T2[1:3,1:3]*T1[1:3,1:3]')))
    end
    return et/N, er/N
end

"""
q, xin, Tn0, Tn0mod, Ta = simulateCalibration_LPOE(N)
Create simulated data for use with `calibLPOE`
"""
function simulateCalibration_LPOE(N)
    Random.seed!(1)
    n       = 6
    q       = 2π*rand(N,n)
    dh      = DH7600()
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, I4);
    xin = DH2twistsLPOE(Tn0)
    Tn0mod = deepcopy(Tn0)
    for i = 1:n+1
        Tn0mod[1:3,4,i] += 0.1randn(3)
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkineLPOE(Tn0,xin,q[i,:])
        AAA,BBB,T = jacobian(q[i,:],dh, I4);
        @assert T ≈ Ta[:,:,i]
    end
    return q, xin, Tn0, Tn0mod, Ta
end

"""
q, xin, Tn0, Tn0mod, Ta = simulateCalibration_LPOE_dual(N)

Create simulated data for use with `calibLPOE_dual`
"""
function simulateCalibration_LPOE_dual(N)
    Random.seed!(1)
    n       = 6
    qt      = 2π*rand(N,n)
    q1      = qt + 0.01*π/180*randn(N,n)
    q2      = qt + 0.01*π/180*randn(N,n)
    q       = cat(q1,q2, dims=3)
    dh      = DH7600()
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, I4);
    xin = DH2twistsLPOE(Tn0)
    xin = cat(xin,xin, dims=3)
    Tn0mod = cat(deepcopy(Tn0),deepcopy(Tn0), dims=4)
    for i = 1:n+1
        Tn0mod[1:3,4,i,1] += 0.1randn(3)
        Tn0mod[1:3,4,i,2] += 0.1randn(3)
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkineLPOE(Tn0,xin[:,:,1],qt[i,:])
    end

    return q, xin, Tn0, Tn0mod, Ta
end

# """
#     xin, et, er = calibPOE(Xin,Ta,q; maxiter=50, λ = 10000.0)
#
# Performs knematic calibration using the POE formulation of kinematics.
# - `xin` Twists
# - `et` Translational errors as function of iteration of algorithm
# - `er` Rotational errors as function of iteration of algorithm
# - `λ` Regularization parameter
#
# See `runtests.jl` for simulated usage. Some functions to generate simulation data for calibration purposes are provided in $(@__FILE__)
# """
# function calibPOE(Xin,Ta,q; maxiter=50, λ = 100.0)
#
#     xin     = copy(Xin)
#     n       = size(xin,2)-2
#     N       = size(q,1)
#     y       = zeros(6N)
#     A       = zeros(6N,6*(n+1))
#     xini    = similar(xin)
#     et      = zeros(maxiter+1)
#     er      = zeros(maxiter+1)
#     et[1],er[1]   = evalErrorPOE(xin,Ta,q)
#     println("Error: ",et[1], " λ: ", λ)
#     @assert size(Ta,3) == N
#
#     for iter = 1:maxiter # do a few iterations of the calibration
#
#         for i = 1:N
#             Tfull = fkinePOE(xin,q[i,:]')
#             # populate the matrices of the linear estimation problem
#             y[(i-1)*6+1:6i] = twistcoords(logT(Ta[:,:,i]*trinv(Tfull)))
#             A[(i-1)*6+1:6*i,:] = Ai(q[i,:],xin)
#
#         end
#         x = [A;λ*I]\[y;zeros(size(A,2))]
#         @show norm(x)
#         # Update the candidate nominal parameters
#         for j = 1:n+1
#             ξ = x[(j-1)*6+1:6j]
#             xini[:,j] = xin[:,j] + ξ
#             # xini[:,j] = conformize(xini[:,j])
#         end
#
#         et[iter+1],er[iter+1] = evalErrorPOE(xini, Ta,q)
#         println("Error: ",round(et[iter+1], digits=5), " λ: ", λ, " Norm dx: ", round(norm(x), digits=5))
#
#         if norm(x) < 1e-10
#             return xin, et, er
#         end
#         if et[iter+1] > et[iter]
#             et[iter+1] = et[iter]
#             λ *= 10
#         else
#             λ/=10
#             xin .= xini
#         end
#
#     end
#
#     return xin, et, er
#
# end

# function calibPOE_offsets_from_points(Xin,Q;maxiter=50, λ = 10000.0)
#     xin     = copy(Xin)
#     n       = size(xin,2)-2
#     Ndatasets = size(Q,1)
#     Npoints = sum(map(x-> size(x,1),Q))
#
#     y       = zeros(6Npoints)
#     A       = zeros(6Npoints,n-1)
#     xini    = similar(xin)
#     et      = zeros(maxiter+1)
#     er      = zeros(maxiter+1)
#     et[1]   = evalErrorPOE_offsets_from_points(xin,Q, zeros(n-1))
#     println("Error: ",et[1], " λ: ", λ)
#     δq      = zeros(n-1)
#
#     for iter = 1:maxiter # do a few iterations of the calibration
#
#         ii = 1
#
#         for p = 1:Ndatasets
#             q = Q[p]
#             N = size(q,1)
#             Tfull = cat(3,[fkinePOE(xin,q[i,:]'+[δq;0]) for i = 1:N]...)
#             Tm = squeeze(mean(Tfull,3),3)
#             orthonormal!(Tm)
#             for i = 1:N
#                 # populate the matrices of the linear estimation problem
#                 y[ii:ii+5] = twistcoords(logT(Tm*trinv(Tfull[:,:,i])))
#                 A[ii:ii+5,:] = xii(q[i,1:(n-1)]'+δq, xin[:,1:(n-1)])
#                 ii += 6
#             end
#         end
#
#         x = (A'A + λ*I)\A'y
#         # Update the candidate nominal parameters
#         δqi = δq + x
#
#
#         et[iter+1] = evalErrorPOE_offsets_from_points(xin, Q, δqi)
#         println("$iter Error: ",round(et[iter+1], digits=5), " λ: ", λ, " Norm dx: ", round(norm(x), digits=5))
#
#         if norm(x) < 1e-10
#             return δq, et
#         end
#         if et[iter+1] > et[iter]
#             et[iter+1] = et[iter]
#             λ *= 10
#         else
#             # λ = max(λ/10,1e-8)
#             δq = deepcopy(δqi)
#         end
#
#     end
#
#     return δq, et
#
# end

# function evalErrorPOE(xin,Ta,q)
#     xini = copy(xin)
#     N = size(Ta,3)
#     et = 0.0
#     er = 0.0
#     for i = 1:N
#         Tfull = fkinePOE(xin,q[i,:])
#         et += norm(Ta[1:3,4,i]-Tfull[1:3,4])
#         er += norm(skewcoords(logR(Ta[1:3,1:3,i]*Tfull[1:3,1:3]')))
#     end
#     return et/N, er/N
# end
#
# function evalErrorPOE_offsets_from_points(xi,Q,dq)
#     Ndatasets = size(Q,1)
#     J = Array{Float64}(Ndatasets)
#     for p = 1:Ndatasets
#         q = Q[p]
#         N = size(q,1)
#         T = Array{typeof(dq[1])}(4,4,N)
#         for i = 1:N
#             T[:,:,i] = fkinePOE(xi,q[i,:] + [dq;0])
#         end
#         xyz = squeeze(T[1:3,4,:],2)
#         mxyz = mean(xyz,2)
#         Jpos = sum(abs((xyz.-mxyz)./mxyz))/N
#
#         Rm = squeeze(mean(T[1:3,1:3,:],3),3)
#         orthonormal!(Rm)
#         Jori = 0.0
#         @inbounds for i = 1:N
#             Jori += Rangle(T[1:3,1:3,i],Rm)
#         end
#         Jori /= Rangle(Rm)
#         Jori /= N
#
#         J[p] = Jpos + Jori
#     end
#     return sum(J)
# end

# """
# q, xin, T0, xinmod, Ta = simulateCalibration_POE(N)
# Create simulated data for use with `calibPOE`
# """
# function simulateCalibration_POE(N)
#     Random.seed!(1)
#     n       = 6
#     q       = 2π*rand(N,n)
#     dh      = DH7600()
#     AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, I4);
#     xin     = DH2twistsPOE(dh)
#     xinmod   = deepcopy(xin)
#     for i = 1:n+1
#         xinmod[:,i]  += [0.001randn(3); 0.01π/180*randn(3)]
#         #conformize(xinmod[:,i])
#     end
#     Ta  = zeros(4,4,N)
#     for i = 1:N
#         Ta[:,:,i] = fkinePOE(xin,q[i,:])
#     end
#     return q, xin, T0, xinmod, Ta
# end

max_forks_repo_name: baggepinnen/Robotlib


text:
using Plots

using simulator.utilities
using simulator.matrixcore
using simulator.fresnelltools

mirror1_r = 0.6
mirror2_r = 0.6
distance = 3e-4

mirror1 = [
    1 - mirror1_r mirror1_r;
    mirror1_r     1 - mirror1_r
]
mirror2 = [
    1 - mirror2_r mirror2_r;
    mirror2_r     1 - mirror2_r
]

gap(ν) = FresnellSlab(1, 2*π*ν/c_0, distance, 0)

transmittance(ν) = CascadeScattering([mirror1, gap(ν), mirror2])[1, 1]
golden(ν) = (1 - mirror1_r) * (1 - mirror2_r) * ℯ^(-im * distance * 2*π*ν / c_0) / (1 - mirror1_r * mirror2_r * ℯ^(-2im * distance * 2*π*ν / c_0))

xs = 1e9:1e8:2e12
transmission = map(ν -> abs(transmittance(ν))^2, xs)
transmission_g = map(ν -> abs(golden(ν))^2, xs)

reflection = 1 .- transmission
reflection_g = 1 .- transmission_g

plot(xs, transmission, label="transmission")
plot!(xs, reflection, label="reflection")
plot!(xs, transmission_g, label="transmission_g")
plot!(xs, reflection_g, label="reflection_g")

ext: jl


text:
# Make jump variables and a dict mapping hold vars to jump vars
function _initialize_hold_variables(trans_model::JuMP.Model,
                                      inf_model::InfiniteOpt.InfiniteModel)
    # search inf_model for hold vars and make a jump var for one that is used
    for index in sort(collect(keys(inf_model.vars)))
        var = inf_model.vars[index]
        if isa(var, InfiniteOpt.HoldVariable)
            hvref = InfiniteOpt.HoldVariableRef(inf_model, index)
            if is_used(hvref)
                vref = JuMP.add_variable(trans_model,
                                         JuMP.ScalarVariable(var.info),
                                         JuMP.name(hvref))
                transcription_data(trans_model).hold_to_var[hvref] = vref
            end
        end
    end
    return
end

# Return a vector of arrays containing the supports
function _list_supports(prefs::Tuple)::Vector
    return [InfiniteOpt.supports(pref) for pref in prefs]
end

# Make an index mapping for parameter support combinations for an infinite variable
function _make_supports(prefs::Tuple)::Vector
    support_list = _list_supports(prefs)
    combos = Iterators.product(support_list...)
    return [combo for combo in Iterators.take(combos, length(combos))]
end

# Make jump variables and a dict mapping infinite/point vars to jump vars
function _initialize_infinite_variables(trans_model::JuMP.Model,
                                        inf_model::InfiniteOpt.InfiniteModel)
    # search inf_model for infinite vars and make a jump var for all of its supports
    for index in sort(collect(keys(inf_model.vars)))
        var = inf_model.vars[index]
        if isa(var, InfiniteOpt.InfiniteVariable)
            ivref = InfiniteOpt.InfiniteVariableRef(inf_model, index)
            if InfiniteOpt.is_used(ivref)
                prefs = InfiniteOpt.parameter_refs(ivref)
                supports = _make_supports(prefs)
                transcription_data(trans_model).infvar_to_supports[ivref] = supports
                vrefs = Vector{JuMP.VariableRef}(undef, length(supports))
                name = InfiniteOpt._root_name(ivref)
                for i in eachindex(vrefs)
                    # TODO Perhaps add different naming options...
                    vrefs[i] = JuMP.add_variable(trans_model,
                                                 JuMP.ScalarVariable(var.info),
                                                 string(name, "(support: ", i, ")"))
                end
                transcription_data(trans_model).infinite_to_vars[ivref] = vrefs
            end
        end
    end
    return
end

# Map the point variable reference to a transcribed variable based off of the support
function _update_point_mapping(trans_model::JuMP.Model,
                               pvref::InfiniteOpt.PointVariableRef,
                               ivref::InfiniteOpt.InfiniteVariableRef,
                               support::Tuple)
    supps = InfiniteOpt.variable_supports(trans_model, ivref)
    # search for the support and update mapping
    for i in eachindex(supps)
       if all(isapprox.(support, supps[i]))
           vrefs = transcription_variable(trans_model, ivref)
           transcription_data(trans_model).point_to_var[pvref] = vrefs[i]
           return
       end
    end
    # fallback that should not be needed
    error("Couldn't find variable to map $pvref to.")
end

# Override the info of the jump variable with the point variable's if any is provided
function _update_point_info(trans_model::JuMP.Model,
                            pvref::InfiniteOpt.PointVariableRef)
    vref = transcription_variable(trans_model, pvref)
    if JuMP.has_lower_bound(pvref) && !JuMP.has_lower_bound(vref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_lower_bound(vref, JuMP.lower_bound(pvref))
    end
    if JuMP.has_upper_bound(pvref) && !JuMP.has_upper_bound(vref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_upper_bound(vref, JuMP.upper_bound(pvref))
    end
    if JuMP.is_fixed(pvref) && !JuMP.is_fixed(vref)
        JuMP.fix(vref, JuMP.fix_value(pvref), force = true)
    end
    if JuMP.is_binary(pvref) && !JuMP.is_binary(vref)
        if JuMP.is_integer(vref)
            JuMP.unset_integer(vref)
        end
        JuMP.set_binary(vref)
    elseif JuMP.is_integer(pvref) && !JuMP.is_integer(vref)
        if JuMP.is_binary(vref)
            JuMP.unset_binary(vref)
        end
        JuMP.set_integer(vref)
    end
    if JuMP.start_value(pvref) != JuMP.start_value(vref)
        JuMP.set_start_value(vref, JuMP.start_value(pvref))
    end
    return
end

# Map point variables to the correct transcribed infinite variable
function _map_point_variables(trans_model::JuMP.Model,
                              inf_model::InfiniteOpt.InfiniteModel)
    # search inf_model for point vars and map them to jump vars if they are used
    for (index, var) in inf_model.vars
        if isa(var, InfiniteOpt.PointVariable)
            pvref = InfiniteOpt.PointVariableRef(inf_model, index)
            if InfiniteOpt.is_used(pvref)
                ivref = InfiniteOpt.infinite_variable_ref(pvref)
                support = InfiniteOpt.parameter_values(pvref)
                _update_point_mapping(trans_model, pvref, ivref, support)
                _update_point_info(trans_model, pvref)
            end
        end
    end
    return
end

# Return the tuple index of a parameter based off of group_id
function _parameter_tuple_index(pref::InfiniteOpt.ParameterRef,
                                prefs::Tuple)::Int
    group = InfiniteOpt.group_id(pref)
    groups = InfiniteOpt._group.(prefs)
    if !(group in groups)
        return -1
    else
        return findfirst(isequal(group), groups)
    end
end

# Return the support value corresponding to a parameter reference
# return NaN is the parameter is not contained in prefs
function _parameter_value(pref::InfiniteOpt.ParameterRef, support::Tuple,
                          prefs::Tuple)::Float64
    pref_index = _parameter_tuple_index(pref, prefs)
    if pref_index == -1
        return NaN
    elseif isa(prefs[pref_index], InfiniteOpt.ParameterRef)
        return support[pref_index]
    else
        for (k, v) in prefs[pref_index].data
            if v == pref
                return support[pref_index].data[k]
            end
        end
    end
    return NaN
end

## Helper function for mapping InfiniteOpt variables to jump variables
# FiniteVariableRef
function _map_to_variable(fvref::InfiniteOpt.FiniteVariableRef, support::Tuple,
                          prefs::Tuple, trans_model::JuMP.Model)::JuMP.VariableRef
    return transcription_variable(trans_model, fvref)
end

# InfiniteVariableRef
function _map_to_variable(ivref::InfiniteOpt.InfiniteVariableRef, support::Tuple,
                          prefs::Tuple, trans_model::JuMP.Model)::JuMP.VariableRef
    # reduce support to only include the relavent parameter id groups
    ivref_groups = InfiniteOpt._group.(InfiniteOpt.parameter_refs(ivref))
    support_groups = InfiniteOpt._group.(prefs)
    reduced_support = Tuple(support[findfirst(isequal(group), support_groups)]
                            for group in ivref_groups)
    # find the jump variable associated with the support
    supps = InfiniteOpt.variable_supports(trans_model, ivref)
    for i in eachindex(supps)
        if all(isapprox.(reduced_support, supps[i]))
            return transcription_variable(trans_model, ivref)[i]
        end
    end
    # this shouldn't be needed, but provided as a backup
    error("Couldn't find JuMP variable corresponding to $ivref.")
end

# ReducedInfiniteVariableRef
function _map_to_variable(rvref::InfiniteOpt.ReducedInfiniteVariableRef,
                          support::Tuple, prefs::Tuple,
                          trans_model::JuMP.Model)::JuMP.VariableRef
    # parse the reduced parameters and modify the support to include them
    ivref = InfiniteOpt.infinite_variable_ref(rvref)
    orig_groups = InfiniteOpt._group.(InfiniteOpt.parameter_refs(ivref))
    support_groups = InfiniteOpt._group.(prefs)
    support = Any[i for i in support]
    prefs = Any[pref for pref in prefs]
    for (index, value) in InfiniteOpt.eval_supports(rvref)
        if orig_groups[index] in support_groups
            support[findfirst(isequal(orig_groups[index]), support_groups)] = value
        else
            push!(support, value)
            push!(prefs, InfiniteOpt.parameter_refs(ivref)[index])
        end
    end
    prefs = Tuple(pref for pref in prefs)
    support = Tuple(i for i in support)
    # call the infinite variable function with the updated support and prefs
    return _map_to_variable(ivref, support, prefs, trans_model)
end

# ParameterRef
function _map_to_variable(pref::InfiniteOpt.ParameterRef, support::Tuple,
                          prefs::Tuple, trans_model::JuMP.Model)::Float64
    # find pref in prefs and return associated support value
    value = _parameter_value(pref, support, prefs)
    # this error is a backup that should not be needed
    value === NaN && error("Couldn't find support corresponding to $pref.")
    return value
end

# Return a Bool whether the support satisfies the bounds or not
function _support_in_bounds(support::Tuple, prefs::Tuple, bounds::Dict)::Bool
    for (pref, set) in bounds
        value = _parameter_value(pref, support, prefs)
        if value === NaN
            continue
        end
        if value < set.lower_bound || value > set.upper_bound
            return false
        end
    end
    return true
end

# Return a bool array indicating which supports are in bounds
function _supports_in_bounds(supports::Vector, prefs::Tuple, bounds::Dict)::Vector
    return [_support_in_bounds(supports[i], prefs, bounds)
            for i in eachindex(supports)]
end

## Convert jump scalar expressions with InfiniteOpt variables into transcribed relations
# PointVariableRef and HoldVariableRef --> return scalar jump object
function _make_transcription_function(vref::InfiniteOpt.FiniteVariableRef,
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple
    return transcription_variable(trans_model, vref), ()
end

# InfiniteVariableRef --> return tuple of expressions, prefs, and support mapping
function _make_transcription_function(vref::InfiniteOpt.InfiniteVariableRef,
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple
    # Update the supports to include only those in the bounds and return vars
    if length(bounds) != 0
        supports = InfiniteOpt.variable_supports(trans_model, vref)
        prefs = InfiniteOpt.parameter_refs(vref)
        old_support_indices = _supports_in_bounds(supports, prefs, bounds)
        new_supports = supports[old_support_indices]
        return transcription_variable(trans_model, vref)[old_support_indices],
               InfiniteOpt.parameter_refs(vref), new_supports
    # easy case just return the info
    else
        return transcription_variable(trans_model, vref),
               InfiniteOpt.parameter_refs(vref),
               InfiniteOpt.variable_supports(trans_model, vref)
    end
end

# ParameterRef --> return tuple of numbers, pref, and support mapping
function _make_transcription_function(pref::InfiniteOpt.ParameterRef,
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple
    supports = unique(InfiniteOpt.supports(pref))
    # truncate the supports if out of bounds
    if pref in keys(bounds)
        old_support_indices = bounds[pref].lower_bound .<= supports .<= bounds[pref].upper_bound
        new_supports = supports[old_support_indices]
        return new_supports, (pref, ), new_supports
    # easy case just return the supports
    else
        return supports, (pref, ), supports
    end
end

# MeasureRef --> return depends if finite or infinite
function _make_transcription_function(mref::InfiniteOpt.MeasureRef,
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple
    # expand the measure
    func = InfiniteOpt.measure_function(mref)
    data = InfiniteOpt.measure_data(mref)
    new_func = InfiniteOpt._possible_convert(InfiniteOpt.FiniteVariableRef,
                                         InfiniteOpt.expand_measure(func, data,
                                            trans_model, _update_point_mapping))
    # will either call finite variable function or general variable function
    return _make_transcription_function(new_func, trans_model, bounds)
end

# GenericAffExpr of FiniteVariableRefs --> return scalar jump object
function _make_transcription_function(expr::JuMP.GenericAffExpr{C,
                                               <:InfiniteOpt.FiniteVariableRef},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    # replace finite vars with jump vars
    pairs = [transcription_variable(trans_model, var) => coef
             for (var, coef) in expr.terms]
    return JuMP.GenericAffExpr(expr.constant,
                         JuMP._new_ordered_dict(JuMP.VariableRef, C, pairs)), ()
end

# GenericQuadExpr of FiniteVariableRefs --> return scalar jump object
function _make_transcription_function(expr::JuMP.GenericQuadExpr{C,
                                               <:InfiniteOpt.FiniteVariableRef},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    # replace finite vars with jump vars
    pairs = Vector{Pair{JuMP.UnorderedPair{JuMP.VariableRef}, C}}(undef,
                                                             length(expr.terms))
    counter = 1
    for k in keys(expr.terms)
        a = transcription_variable(trans_model, k.a)
        b = transcription_variable(trans_model, k.b)
        pairs[counter] = JuMP.UnorderedPair(a, b) => expr.terms[k]
        counter += 1
    end
    aff = _make_transcription_function(expr.aff, trans_model)[1]
    return JuMP.GenericQuadExpr(aff,
     JuMP._new_ordered_dict(JuMP.UnorderedPair{JuMP.VariableRef}, C, pairs)), ()
end

# GenericAffExpr of GeneralVariableRefs --> return tuple of numbers, pref,
# and support mapping
function _make_transcription_function(expr::JuMP.GenericAffExpr{C,
                                              <:InfiniteOpt.GeneralVariableRef},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    # check if there is only 1 var to dispatch to that transcription method
    if length(expr.terms) == 1
        var = first(keys(expr.terms))
        results = _make_transcription_function(var, trans_model, bounds)
        # results[2] contains prefs if is infinite or it is empty otherwise
        if length(results[2]) != 0
            exprs = expr.terms[var] * results[1] .+ expr.constant
            return exprs, results[2], results[3]
        else
            return expr.terms[var] * results[1] + expr.constant, ()
        end
    end
    # check to see if there are measures and expand them
    # TODO rework paradigm so expressions don't need to be checked for measures
    has_measure = false
    for var in keys(expr.terms)
        if var isa InfiniteOpt.MeasureRef
            has_measure = true
            break
        end
    end
    if has_measure
        expr = InfiniteOpt._expand_measures(expr, trans_model,
                                            _update_point_mapping)
    end
    # dispatch to quadratic method if the measures contained quadratic terms
    if isa(expr, JuMP.GenericQuadExpr)
        return _make_transcription_function(expr, trans_model, bounds)
    end
    # determine the common set of prefs and make all of the support combos
    prefs = InfiniteOpt._all_parameter_refs(expr)
    supports = _make_supports(prefs)
    if length(bounds) != 0
        old_support_indices = _supports_in_bounds(supports, prefs, bounds)
        supports = supports[old_support_indices]
    end
    exprs = [zero(JuMP.GenericAffExpr{C, JuMP.VariableRef})
             for i in eachindex(supports)]
    # make an expression for each support
    for i in eachindex(exprs)
        exprs[i].constant = expr.constant
        for (var, coef) in expr.terms
            # replace each variable with appropriate jump var
            new_var = _map_to_variable(var, supports[i], prefs,
                                       trans_model)
            if isa(new_var, JuMP.VariableRef)
                JuMP.add_to_expression!(exprs[i], coef, new_var)
            else
                JuMP.add_to_expression!(exprs[i], coef * new_var)
            end
        end
    end
    return exprs, prefs, supports
end

# GenericQuadExpr of GeneralVariableRefs
function _make_transcription_function(expr::JuMP.GenericQuadExpr{C,
                                              <:InfiniteOpt.GeneralVariableRef},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    # check if there is only 1 var to dispatch to that transcription method
    if length(expr.terms) == 0 && length(expr.aff.terms) == 1
      var = first(keys(expr.aff.terms))
      results = _make_transcription_function(var, trans_model, bounds)
      # results[2] contains prefs if is infinite or it is empty otherwise
      if length(results[2]) != 0
          exprs = expr.aff.terms[var] * results[1] .+ expr.aff.constant
          return exprs, results[2], results[3]
      else
          return expr.aff.terms[var] * results[1] + expr.aff.constant, ()
      end
    end
    # check to see if there are measures and expand them
    # TODO rework paradigm so expressions don't need to be checked for measures
    all_variables = InfiniteOpt._all_function_variables(expr)
    has_measure = false
    for var in all_variables
        if var isa InfiniteOpt.MeasureRef
            has_measure = true
            break
        end
    end
    if has_measure
        # TODO expand measures in place by deleting entry and appending to the end
        # TODO thus the measure search can be done by expand measures
        expr = InfiniteOpt._expand_measures(expr, trans_model,
                                            _update_point_mapping)
    end
    # determine the common set of prefs and make all of the support combos
    prefs = InfiniteOpt._all_parameter_refs(expr)
    supports = _make_supports(prefs)
    if length(bounds) != 0
        old_support_indices = _supports_in_bounds(supports, prefs, bounds)
        supports = supports[old_support_indices]
    end
    exprs = [zero(JuMP.GenericQuadExpr{C, JuMP.VariableRef})
             for i in eachindex(supports)]
    # make an expression for each support
    for i in eachindex(exprs)
        exprs[i].aff.constant = expr.aff.constant
        for (var, coef) in expr.aff.terms
            # replace each variable with appropriate jump var
            new_var = _map_to_variable(var, supports[i], prefs,
                                       trans_model)
            if isa(new_var, JuMP.VariableRef)
                JuMP.add_to_expression!(exprs[i], coef, new_var)
            else
                JuMP.add_to_expression!(exprs[i], coef * new_var)
            end
        end
        for (pair, coef) in expr.terms
            # replace each variable with appropriate jump var
            var_a = _map_to_variable(pair.a, supports[i], prefs,
                                     trans_model)
            var_b = _map_to_variable(pair.b, supports[i], prefs,
                                     trans_model)
            if isa(var_a, JuMP.VariableRef) && isa(var_b, JuMP.VariableRef)
                JuMP.add_to_expression!(exprs[i], coef, var_a, var_b)
            elseif isa(var_b, JuMP.VariableRef)
                JuMP.add_to_expression!(exprs[i], coef * var_a, var_b)
            elseif isa(var_a, JuMP.VariableRef)
                JuMP.add_to_expression!(exprs[i], coef * var_b, var_a)
            else
                JuMP.add_to_expression!(exprs[i], coef * var_a * var_b)
            end
        end
    end
    # TODO Check if there are still quadratic terms and convert to AffExpr if necessary
    return exprs, prefs, supports
end

# GenericAffExpr and GenericQuadExpr of MeasureFiniteVariableRefs --> return
# depends on whether finite
function _make_transcription_function(expr::Union{JuMP.GenericAffExpr{C,
                                        <:InfiniteOpt.MeasureFiniteVariableRef},
                                        JuMP.GenericQuadExpr{C,
                                        <:InfiniteOpt.MeasureFiniteVariableRef}},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    expr = InfiniteOpt._possible_convert(FiniteVariableRef,
                                         InfiniteOpt._expand_measures(expr,
                                         trans_model, _update_point_mapping))
    return _make_transcription_function(expr, trans_model, bounds)
end

# Empty jump variable expr (for constraints of form number <= number)
function _make_transcription_function(expr::JuMP.GenericAffExpr{C,
                                                              JuMP.VariableRef},
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple where {C}
    return expr, ()
end

# Fall back function for other jump objects
function _make_transcription_function(expr::JuMP.AbstractJuMPScalar,
                                      trans_model::JuMP.Model,
                                      bounds::Dict = Dict())::Tuple
    type = typeof(expr)
    error("Unsupported transcription of expression of type $type.")
end

## Construct the objective and error is contains non finite variables
function _set_objective(trans_model::JuMP.Model,
                        inf_model::InfiniteOpt.InfiniteModel)
    trans_obj, = _make_transcription_function(JuMP.objective_function(inf_model),
                                              trans_model)
    isa(trans_obj, Vector) && error("Objective is not finite, ensure all " *
                                    "infinite variables/parameters in measures " *
                                    "are evaluated completely.")
    JuMP.set_objective(trans_model, JuMP.objective_sense(inf_model), trans_obj)
    return
end

## Define helper functions for setting constraint mappings
# InfiniteConstraintRef
function _set_mapping(icref::InfiniteOpt.InfiniteConstraintRef,
                      crefs::Vector{<:JuMP.ConstraintRef})
    transcription_data(JuMP.owner_model(crefs[1])).infinite_to_constrs[icref] = crefs
    return
end

# MeasureConstraintRef (infinite)
function _set_mapping(mcref::InfiniteOpt.MeasureConstraintRef,
                      crefs::Vector{<:JuMP.ConstraintRef})
    transcription_data(JuMP.owner_model(crefs[1])).measure_to_constrs[mcref] = crefs
    return
end

# MeasureConstraintRef (finite)
function _set_mapping(mcref::InfiniteOpt.MeasureConstraintRef,
                      cref::JuMP.ConstraintRef)
    transcription_data(JuMP.owner_model(cref)).measure_to_constrs[mcref] = [cref]
    return
end

# FiniteConstraintRef
function _set_mapping(fcref::InfiniteOpt.FiniteConstraintRef,
                      cref::JuMP.ConstraintRef)
    transcription_data(JuMP.owner_model(cref)).finite_to_constr[fcref] = cref
    return
end

## Define helper functions for setting constraint supports
# InfiniteConstraintRef
function _set_supports(trans_model::JuMP.Model,
                       icref::InfiniteOpt.InfiniteConstraintRef,
                       supports::Vector)
    transcription_data(trans_model).infconstr_to_supports[icref] = supports
    return
end

# MeasureConstraintRef
function _set_supports(trans_model::JuMP.Model,
                       mcref::InfiniteOpt.MeasureConstraintRef,
                       supports::Vector)
    transcription_data(trans_model).measconstr_to_supports[mcref] = supports
    return
end

## Define helper functions for setting constraint parameter reference tuples
# InfiniteConstraintRef
function _set_parameter_refs(trans_model::JuMP.Model,
                             icref::InfiniteOpt.InfiniteConstraintRef,
                             prefs::Tuple)
    transcription_data(trans_model).infconstr_to_params[icref] = prefs
    return
end

# MeasureConstraintRef
function _set_parameter_refs(trans_model::JuMP.Model,
                             mcref::InfiniteOpt.MeasureConstraintRef,
                             prefs::Tuple)
    transcription_data(trans_model).measconstr_to_params[mcref] = prefs
    return
end

# Extract the root name of a constraint
function InfiniteOpt._root_name(cref::InfiniteOpt.GeneralConstraintRef)
    name = JuMP.name(cref)
    if length(name) == 0
        return "noname"
    else
        first_bracket = findfirst(isequal('['), name)
        if first_bracket == nothing
            return name
        else
            return name[1:first_bracket-1]
        end
    end
end

## leverage _make_transcription_function to transcribe the constraints
function _set_constraints(trans_model::JuMP.Model, inf_model::InfiniteOpt.InfiniteModel)
    # transform and add constraints that haven't already been added through add_variable
    ordered_indices = sort(collect(keys(inf_model.constrs)))
    for index in ordered_indices
        constr = inf_model.constrs[index]
        if !inf_model.constr_in_var_info[index]
            # extract the reference and transcribe the jump object function
            icref = InfiniteOpt._make_constraint_ref(inf_model, index)
            if isa(constr, BoundedScalarConstraint)
                results = _make_transcription_function(constr.func, trans_model,
                                                       constr.bounds.intervals)
            else
                results = _make_transcription_function(constr.func, trans_model)
            end
            # extract the name and create a transcribed constraint
            root_name = InfiniteOpt._root_name(icref)
            # results[2] contains prefs if is infinite or it is empty otherwise
            if length(results[2]) != 0
                crefs = Vector{JuMP.ConstraintRef}(undef, length(results[1]))
                for i in eachindex(crefs)
                    con = JuMP.build_constraint(error, results[1][i], constr.set)
                    # TODO Perhaps improve naming
                    name = string(root_name, "(Support: ", i, ")")
                    crefs[i] = JuMP.add_constraint(trans_model, con, name)
                end
                # update the mappings
                _set_mapping(icref, crefs)
                _set_parameter_refs(trans_model, icref, results[2])
                _set_supports(trans_model, icref, results[3])
            else
                # We have a finite constraint, just add normally and update mapping.
                con = JuMP.build_constraint(error, results[1], constr.set)
                cref = JuMP.add_constraint(trans_model, con, root_name)
                _set_mapping(icref, cref)
            end
        end
    end
    return
end

## Helper functions for mapping InfiniteOpt var info constrs to the transcribed ones
# FiniteVariableRef
function _map_info_constraints(fvref::InfiniteOpt.FiniteVariableRef,
                               trans_model::JuMP.Model)
    vref = transcription_variable(trans_model, fvref)
    # Check if variables have info constraints and map if they do
    if JuMP.has_lower_bound(fvref)
        _set_mapping(JuMP.LowerBoundRef(fvref), JuMP.LowerBoundRef(vref))
    end
    if JuMP.has_upper_bound(fvref)
        _set_mapping(JuMP.UpperBoundRef(fvref), JuMP.UpperBoundRef(vref))
    end
    if JuMP.is_fixed(fvref)
        _set_mapping(JuMP.FixRef(fvref), JuMP.FixRef(vref))
    end
    if JuMP.is_integer(fvref)
        _set_mapping(JuMP.IntegerRef(fvref), JuMP.IntegerRef(vref))
    end
    if JuMP.is_binary(fvref)
        _set_mapping(JuMP.BinaryRef(fvref), JuMP.BinaryRef(vref))
    end
    return
end

# InfiniteVariableRef
function _map_info_constraints(ivref::InfiniteOpt.InfiniteVariableRef,
                               trans_model::JuMP.Model)
    vrefs = transcription_variable(trans_model, ivref)
    # TODO Prealocate the size
    # Check if both variables have a constraint and map if they do
    if JuMP.has_lower_bound(ivref)
        crefs = JuMP.ConstraintRef[]
        supports = Tuple[]
        for i in eachindex(vrefs)
            if JuMP.has_lower_bound(vrefs[i])
                push!(crefs, JuMP.LowerBoundRef(vrefs[i]))
                push!(supports, InfiniteOpt.variable_supports(trans_model, ivref)[i])
            end
        end
        if length(crefs) != 0
            _set_mapping(JuMP.LowerBoundRef(ivref), crefs)
            _set_parameter_refs(trans_model, JuMP.LowerBoundRef(ivref),
                                InfiniteOpt.parameter_refs(ivref))
            _set_supports(trans_model, JuMP.LowerBoundRef(ivref), supports)
        end
    end
    if JuMP.has_upper_bound(ivref)
        crefs = JuMP.ConstraintRef[]
        supports = Tuple[]
        for i in eachindex(vrefs)
            if JuMP.has_upper_bound(vrefs[i])
                push!(crefs, JuMP.UpperBoundRef(vrefs[i]))
                push!(supports, InfiniteOpt.variable_supports(trans_model, ivref)[i])
            end
        end
        if length(crefs) != 0
            _set_mapping(JuMP.UpperBoundRef(ivref), crefs)
            _set_parameter_refs(trans_model, JuMP.UpperBoundRef(ivref),
                                InfiniteOpt.parameter_refs(ivref))
            _set_supports(trans_model, JuMP.UpperBoundRef(ivref), supports)
        end
    end
    if JuMP.is_fixed(ivref)
        crefs = JuMP.ConstraintRef[]
        supports = Tuple[]
        for i in eachindex(vrefs)
            if JuMP.is_fixed(vrefs[i])
                push!(crefs, JuMP.FixRef(vrefs[i]))
                push!(supports, InfiniteOpt.variable_supports(trans_model, ivref)[i])
            end
        end
        if length(crefs) != 0
            _set_mapping(JuMP.FixRef(ivref), crefs)
            _set_parameter_refs(trans_model, JuMP.FixRef(ivref),
                                InfiniteOpt.parameter_refs(ivref))
            _set_supports(trans_model, JuMP.FixRef(ivref), supports)
        end
    end
    if JuMP.is_integer(ivref)
        crefs = JuMP.ConstraintRef[]
        supports = Tuple[]
        for i in eachindex(vrefs)
            if JuMP.is_integer(vrefs[i])
                push!(crefs, JuMP.IntegerRef(vrefs[i]))
                push!(supports, InfiniteOpt.variable_supports(trans_model, ivref)[i])
            end
        end
        if length(crefs) != 0
            _set_mapping(JuMP.IntegerRef(ivref), crefs)
            _set_parameter_refs(trans_model, JuMP.IntegerRef(ivref),
                                InfiniteOpt.parameter_refs(ivref))
            _set_supports(trans_model, JuMP.IntegerRef(ivref), supports)
        end
    end
    if JuMP.is_binary(ivref)
        crefs = JuMP.ConstraintRef[]
        supports = Tuple[]
        for i in eachindex(vrefs)
            if JuMP.is_binary(vrefs[i])
                push!(crefs, JuMP.BinaryRef(vrefs[i]))
                push!(supports, InfiniteOpt.variable_supports(trans_model, ivref)[i])
            end
        end
        if length(crefs) != 0
            _set_mapping(JuMP.BinaryRef(ivref), crefs)
            _set_parameter_refs(trans_model, JuMP.BinaryRef(ivref),
                                InfiniteOpt.parameter_refs(ivref))
            _set_supports(trans_model, JuMP.BinaryRef(ivref), supports)
        end
    end
    return
end

## Map the variable info constraints between the two models
function _map_variable_info_constraints(trans_model::JuMP.Model,
                                        inf_model::InfiniteOpt.InfiniteModel)
    for index in keys(inf_model.vars)
        ivref = InfiniteOpt._make_variable_ref(inf_model, index)
        if InfiniteOpt.is_used(ivref)
            _map_info_constraints(ivref, trans_model)
        end
    end
    return
end

# Fill in a given trans_model with the reformulation of the inf_model
# This assumes that trans_model is empty
function _build_transcription_model!(trans_model::JuMP.Model,
                                     inf_model::InfiniteOpt.InfiniteModel)
    InfiniteOpt.fill_in_supports!(inf_model)
    _initialize_hold_variables(trans_model, inf_model)
    _initialize_infinite_variables(trans_model, inf_model)
    _map_point_variables(trans_model, inf_model)
    if JuMP.objective_sense(inf_model) != MOI.FEASIBILITY_SENSE
        _set_objective(trans_model, inf_model)
    end
    _map_variable_info_constraints(trans_model, inf_model)
    # TODO optimize performance --> bottlenecked with _make_transcription_function
    _set_constraints(trans_model, inf_model)
end

"""
    TranscriptionModel(model::InfiniteModel, [optimizer_constructor;
                       caching_mode::MOIU.CachingOptimizerMode = MOIU.AUTOMATIC,
                       bridge_constraints::Bool = true])

Return a `TranscriptionModel` of `model`. This transcribes all of the variables,
constraints, and objective.

**Example**
```julia-repl
julia> TranscriptionModel(model)
A JuMP Model
Feasibility problem with:
Variables: 130
`GenericAffExpr{Float64,VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 25 constraints
`GenericAffExpr{Float64,VariableRef}`-in-`MathOptInterface.GreaterThan{Float64}`: 100 constraint
`GenericAffExpr{Float64,VariableRef}`-in-`MathOptInterface.LessThan{Float64}`: 84 constraints
`VariableRef`-in-`MathOptInterface.EqualTo{Float64}`: 40 constraints
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 40 constraints
`VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 25 constraints
`VariableRef`-in-`MathOptInterface.Integer`: 40 constraints
`VariableRef`-in-`MathOptInterface.ZeroOne`: 40 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
function TranscriptionModel(inf_model::InfiniteOpt.InfiniteModel,
                            optimizer_constructor = nothing;
                            kwargs...)::JuMP.Model
    if isa(optimizer_constructor, Nothing)
        trans_model = TranscriptionModel(; kwargs...)
    else
        trans_model = TranscriptionModel(optimizer_constructor; kwargs...)
    end
    _build_transcription_model!(trans_model, inf_model)
    return trans_model
end


text:
struct PursuerEvador3DMissile <: AbstractMissile
    pursuer::PointMass3DMissile
    evador::PointMass3DMissile
end

function State(env::PursuerEvador3DMissile)
    @unpack pursuer, evador = env
    return function (x0_pursuer, x0_evador)
        ComponentArray(pursuer=x0_pursuer, evador=x0_evador)
    end
end

function Dynamics!(env::PursuerEvador3DMissile)
    @unpack pursuer, evador = env
    @Loggable function dynamics!(dx, x, params, t; u_pursuer, u_evador)
        @onlylog r = norm(x.pursuer.p - x.evador.p)
        @nested_log :pursuer Dynamics!(pursuer)(dx.pursuer, x.pursuer, nothing, t; u=u_pursuer)
        @nested_log :evador Dynamics!(evador)(dx.evador, x.evador, nothing, t; u=u_evador)
    end
end


text:
# ---
# title: 258. Add Digits
# id: problem258
# author: zhwang
# date: 2022-03-03
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/add-digits/description/>
# hidden: true
# ---
#
# Given a non-negative integer `num`, repeatedly add all its digits until the
# result has only one digit.
#
# **Example:**
#
#
#
#     Input: 38
#     Output: 2
#     Explanation: The process is like: 3 + 8 = 11, 1 + 1 = 2.
#                  Since 2 has only one digit, return it.
#
#
# **Follow up:**
# Could you do it without any loop/recursion in O(1) runtime?
#
#
## @lc code=start
using LeetCode

function add_digits(num::Int)::Int
    num <= 9 && return num
    return iszero(num % 9) ? 9 : num % 9
end
## @lc code=end


text:
"""
Tools to compute partition functions.
"""

"""
    partition_function(
        g::Array{<: Real}, χ::Array{<: Unitful.Energy, 1}, temp::Unitful.Temperature
    )
    partition_function(data::AtomicStage, temp::Unitful.Temperature)

Calculate partition function from a set of levels defined by arrays
of statistical weights and energies (or an AtomicStage struct),
for a given temperature.
"""
function partition_function(
    g::Array{<: Real}, χ::Array{<: Unitful.Energy, 1}, temp::Unitful.Temperature,
)
    return sum(g .* exp.(-χ ./ (k_B * temp)))
end

function partition_function(atom::AtomicStage, temp::Unitful.Temperature)
    return partition_function(atom.g, atom.χ, temp)
end

"""
    partition_function_interpolator(
        atom::AtomicStage, temperatures::Array{<: Unitful.Temperature},
    )
    partition_function_interpolator(
        element::String, stage, temperatures::Array{<: Unitful.Temperature};
        source="NIST"
    )

Returns an interpolator for the partition function of a given atomic stage.

# Arguments
- `atom`: an `AtomicStage` structure with the data for the given stage
- `temperatures`: array with temperatures used to build the interpolation table

# Returns
- `interpolator`: linear interpolator object that takes temperature to give
  partition function
"""
function partition_function_interpolator(
    atom::AtomicStage, temperatures::Array{<: Unitful.Temperature},
)
    pfunc = partition_function.(Ref(atom), temperatures)
    return LinearInterpolation(temperatures, pfunc, extrapolation_bc=Line())
end

function partition_function_interpolator(
    element::String, stage, temperatures::Array{<: Unitful.Temperature};
    source="NIST"
)
    if source == "NIST"
        atom = read_NIST(element, stage)
    else
        error("NotImplemented: Atomic data source $source not supported")
    end
    return partition_function_interpolator(atom, temperatures)
end


text:
using PointerStructs
using Documenter

DocMeta.setdocmeta!(PointerStructs, :DocTestSetup, :(using PointerStructs); recursive=true)

makedocs(;
    modules=[PointerStructs],
    authors="Julia Computing and contributors",
    repo="https://github.com/JuliaSIMD/PointerStructs.jl/blob/{commit}{path}#{line}",
    sitename="PointerStructs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSIMD.github.io/PointerStructs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSIMD/PointerStructs.jl",
)


text:
"Problem 14 in the Hock-Schittkowski suite"
function hs14()

  nlp = Model()

  @variable(nlp, x[i=1:2])
  setvalue(x[1], 2)
  setvalue(x[2], 2)

  @NLobjective(
    nlp,
    Min,
    (x[1] - 2)^2 + (x[2] - 1)^2
  )

  @NLconstraint(
    nlp,
    -x[1]^2/4 - x[2]^2 + 1 >= 0
  )

  @constraint(
    nlp,
    x[1] - 2 * x[2] + 1 == 0
  )

  return nlp
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/CH4cycle.jl")

    addch4cycle(m, use_permafrost)

    set_param!(m, :CH4Cycle, :e_globalCH4emissions, readpagedata(m,"test/validationdata/$valdir/e_globalCH4emissions.csv"))
    set_param!(m, :CH4Cycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/$valdir/rtl_g_landtemperature.csv"))
    set_param!(m, :CH4Cycle, :permtce_permafrostemissions, readpagedata(m,"test/validationdata/$valdir/perm_tot_ce_ch4.csv"))
    set_param!(m,:CH4Cycle,:y_year_0,2015.)

    p = load_parameters(m)
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    #running Model
    run(m)

    conc=m[:CH4Cycle,  :c_CH4concentration]
    conc_compare=readpagedata(m,"test/validationdata/$valdir/c_ch4concentration.csv")

    @test conc ≈ conc_compare rtol=1e-4
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
"""
    Environment(atmos::Atmosphere, wind::Wind, grav::Gravity)

Environment information composed of atmosphere, wind and gravity.

# Constructors
    Environment(atmos::Atmosphere, wind::Wind, grav::Gravity)
    Environment(pos::Position; atmos="ISA1978", wind="NoWind", grav="const")
"""
struct Environment
    atmos::Atmosphere
    wind::Wind
    grav::Gravity
end

"""
    Environment(pos::Position; atmos="ISA1978", wind="NoWind", grav="const")

Create an environment for the given position. Models can be chosen for:
    - atmos: atmosphere model.
        - ISA1978: default. International Standard Atmosphere 1978.
        - ISAF16: International Standard Atmosphere 1978 in F16 Stevens model.
    - wind: wind model.
        - NoWind: default.
    - grav: gravity model.
        - const: constant vertical gravity.
"""
function Environment(pos::Position; atmos="ISA1978", wind="NoWind", grav="const")
    # atm
    atmos = uppercase(atmos)
    if atmos == "ISA1978"
        atm_model = AtmosphereISA(pos)
    elseif atmos == "ISAF16"
        atm_model = AtmosphereF16(pos)
    else
        error("Unknown atmospheric model ($atmos).")
    end

    # wind
    wind = uppercase(wind)
    if wind == "NOWIND"
        wind_model = ConstantWind()
    else
        error("Unknown wind model ($wind).")
    end

    # Gravity
    grav = uppercase(grav)
    if grav == "CONST"
        grav_model = EarthConstantGravity()
    else
        error("Unknown gravity model ($grav).")
    end

    Environment(atm_model, wind_model, grav_model)
end

get_atmos(env::Environment) = env.atmos
get_temperature(env::Environment) = get_temperature(get_atmos(env))
get_pressure(env::Environment) = get_pressure(get_atmos(env))
get_density(env::Environment) = get_density(get_atmos(env))
get_sound_velocity(env::Environment) = get_sound_velocity(get_atmos(env))

get_wind(env::Environment) = env.wind
get_wind_dir_int_ver(env::Environment) = get_wind_dir_int_ver(get_wind(env))
get_wind_NED(env::Environment) = get_wind_NED(get_wind(env))
get_wind_body(env::Environment, att::Attitude) = get_wind_body(get_wind(env))
get_wind_direction(env::Environment) = get_wind_direction(get_wind(env))
get_wind_intensity(env::Environment) = get_wind_intensity(get_wind(env))
get_wind_vertical(env::Environment) = get_wind_vertical(get_wind(env))

get_gravity(env::Environment) = env.grav
get_gravity_horizon(env::Environment) = get_gravity_horizon(get_gravity(env))
get_gravity_accel(env::Environment) = get_gravity_accel(get_gravity(env))
get_gravity_body(env::Environment, att::Attitude) = get_gravity_body(get_gravity(env), att)

function calculate_environment(env::Environment, pos::Position)
    atmos = calculate_atmosphere(get_atmos(env), pos)
    wind = calculate_wind(get_wind(env), pos)
    grav = calculate_gravity(get_gravity(env), pos)
    Environment(atmos, wind, grav)
end


text:
@testset "search_ora" begin
    categories = ["Production Process"]
    cat_apps = Dict("Production Process" => ["Example A"])
    # test that all categories can be searched
    for c in categories
        @test search_ora(c) == cat_apps[c]
    end
    # test that uppercase arg is fine
    @test search_ora(uppercase(categories[1])) == cat_apps[categories[1]]

    # test that lowercase arg is fine
    @test search_ora(lowercase(categories[1])) == cat_apps[categories[1]]

    # test that throws ArgumentError when cat not exist
    @test_throws ArgumentError search_ora("This shouldn't exist")

    # test that works with no argument (lists categories)
    @test search_ora() == categories

end


text:
using Particles
using Test

@testset "Particles.jl" begin
    # Write your tests here.
end


text:
# Discrete univariate

convert(::Type{Binomial}, d::Bernoulli) = Binomial(1, d.p)

# Continuous univariate

convert(::Type{Gamma}, d::Exponential) = Gamma(1.0, d.θ)
convert(::Type{Gamma}, d::Erlang) = Gamma(d.α, d.θ)


text:
using Test

@testset "masked_relu.jl" begin

    @testset "MaskedReLU" begin
        mask = [-1, 1, 0, 0, 1, 1]
        @testset "Initialize without tightening algorithm" begin
            p = MaskedReLU(mask)
            @test p.tightening_algorithm === nothing
        end

        @testset "Initialize with tightening algorithm" begin
            p = MaskedReLU(mask, MIPVerify.mip)
            @test !(p.tightening_algorithm === nothing)
            @test p.tightening_algorithm == MIPVerify.mip
        end

        @testset "Base.show" begin
            p = MaskedReLU(mask)
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) ==
                  "MaskedReLU with expected input size (6,). (1 zeroed, 3 as-is, 2 rectified)."
        end
    end

end


text:
module DatasetsSynthetic

using Random
export make_blobs, make_circles, make_moons

uniform_sample_in_zero_maxval(p, maxval) = maxval .* (1 .- rand(p))
uniform_sample_in_minval_maxval(p, minval, maxval) = (maxval-minval) .* rand(p) .+ minval
uniform_sample_in_minval_maxval(n, p, minval, maxval) = (maxval-minval) .* rand(n, p) .+ minval
normal_sample(p, mu, var) = mu .+ sqrt(var) .* randn(p)

"""
Shuffles the rows of an Array `X` and the values of a vector `y` using a randomly
generated permutation. The same permutation is used to shuffle both `X` and `y`.
"""
function shuffle_Xy(X, y; random_seed=Random.GLOBAL_RNG)
    Random.seed!(random_seed)
    perm = randperm(length(y))
    return X[perm,:], y[perm]
end

"""
make_blobs(n::Int=100;
           p::Int=2,
           centers::Int=3,
           cluster_std=1.0,
           center_box=(-10.,10.),
           element_type=Float64,
           random_seed=Random.GLOBAL_RNG,
           return_centers=false)

Generates a dataset with `n` examples of dimension `p` and returns a vector containing
as integers the membership of the different points generated.

The data is roughly grouped around several `centers`  which are created using `cluster_std`.
The data lives inside `center_box` in the case it is randomly generated.

- If `centers` is an integer the centroids are created randomly.
- If `centers` is an Array containing points the centroids are picked from `centers`.
- If `return_centers=true` the centroids of the blobs are returned.
"""
function make_blobs(n::Int=100; p::Int=2,
                    shuffle::Bool=false, centers::Int=3, cluster_std::Real=1.0, center_box=(-10.,10.),
                    element_type=Float64, random_seed=Random.GLOBAL_RNG, return_centers=false, verbose=0)

    Random.seed!(random_seed)
    X = zeros(n, p)
    y = zeros(n)

    if typeof(centers) <: Int
        n_centers = centers
        center_sample = uniform_sample_in_minval_maxval(n, p, center_box[1], center_box[2])
    else
        n_centers = length(centers)
    end

    if typeof(cluster_std) <: AbstractFloat
        cluster_std = cluster_std * randn(n_centers)
    end

    # generates div(n, n_centers) examples assigned to each center
    n_per_center = [div(n, n_centers) for x  in 1:n_centers]

    # adds the reamainding examples to each center up to n
    n_per_center = fill(div(n,n_centers), n_centers)
    n_per_center[end] += rem(n,n_centers)

    # generates the actual vectors close to each center blob
    start_ind = 1
    for (i, (n_blob, std, center)) in enumerate(zip(n_per_center, cluster_std, centers))
        ind_center = start_ind:(start_ind + n_per_center[i]-1)
        X[ind_center,:] .= center' .+ std .* randn(element_type, (n_per_center[i], p));
        y[ind_center] .= i
        if verbose>0
            println("center $i with $(n_per_center[i]) points created")
        end
        start_ind += n_per_center[i]
    end

    if shuffle
       X, y = shuffle_Xy(X, y ; random_seed=random_seed)
    end

    if return_centers
        return X, y, centers
    else
        return X, y
    end
end

"""
make_circles(n::Int=100; shuffle::Bool=true, noise::Number=0., random_seed=Random.GLOBAL_RNG, factor::Number=0.8)

Generates a dataset with `n` bi-dimensional examples. Samples are created
from two circles. One of the circles inside the other. The `noise` scalar
can be used to add noise to the generation process. The scalar `factor`
corresponds to the radius of the smallest circle.
"""
function make_circles(n::Int=100; shuffle::Bool=true, noise::Number=0., random_seed=Random.GLOBAL_RNG, factor::Number=0.8)

    Random.seed!(random_seed)

    n_out = div(n, 2)
    n_in = n - n_out

    linrange_out = Array(LinRange(0, 2 * pi, n_out))
    linrange_in  = Array(LinRange(0, 2 * pi, n_in))

    outer_circ_x = cos.(linrange_out)
    outer_circ_y = sin.(linrange_out)
    inner_circ_x = cos.(linrange_in) .* factor
    inner_circ_y = sin.(linrange_in) .* factor

    X = hcat(vcat(outer_circ_x, inner_circ_x), vcat(outer_circ_y, inner_circ_y))
    y = vcat(ones(Int,n_out), 2*ones(Int,n_in))

    if shuffle
       X, y = shuffle_Xy(X, y ; random_seed=random_seed)
    end

    X .+= noise .* rand(n, 2)

    return X, y
end

"""
make_moons(n::Int=100; shuffle::Bool=true, noise::Number=0.,
           translation::Number=0.5, factor::Number=1.0, random_seed=Random.GLOBAL_RNG)

Generates `n` examples sampling from two moons. The `noise` can be changed to add
noise to the samples.
"""
function make_moons(n::Int=100; shuffle::Bool=true, noise::Number=0.,
                   translation::Number=0.5, factor::Number=1.0, random_seed=Random.GLOBAL_RNG)

    Random.seed!(random_seed)

    n_out = div(n, 2)
    n_in = n - n_out

    linrange_out = Array(LinRange(0, pi, n_out))
    linrange_in  = Array(LinRange(0, pi, n_in))

    outer_circ_x = cos.(linrange_out)
    outer_circ_y = sin.(linrange_out)
    inner_circ_x = 1 .- cos.(linrange_in) .* factor
    inner_circ_y = 1 .- sin.(linrange_in) .* factor .- translation

    X = hcat(vcat(outer_circ_x, inner_circ_x), vcat(outer_circ_y, inner_circ_y))
    y = vcat(ones(Int,n_out), 2*ones(Int,n_in))

    if shuffle
       X, y = shuffle_Xy(X, y ; random_seed=random_seed)
    end

    X .+= noise .* rand(n, 2)

    return X, y
end

end #module


text:
    const user_agents = [
        "chrome", "firefox", "internet_explorer", "opera", "safari"
    ]

    const windows_platform_tokens = [
        "Windows 95", "Windows 98", "Windows 98; Win 9x 4.90", "Windows CE",
        "Windows NT 4.0", "Windows NT 5.0", "Windows NT 5.01",
        "Windows NT 5.1", "Windows NT 5.2", "Windows NT 6.0", "Windows NT 6.1",
        "Windows NT 6.2"
    ]

    const linux_processors = ["i686", "x86_64"]

    const mac_processors = ["Intel", "PPC", "U; Intel", "U; PPC"]

    const langs = ["en-US", "sl-SI", "it-IT"]

    mac_processor()::String = rand(mac_processors)

    linux_processor()::String = rand(linux_processors)

    user_agent()::String = rand(user_agents)

    windows_platform_token()::String = rand(windows_platform_tokens)

    linux_platform_token()::String = string("X11; Linux", linux_processor())

    mac_platform_token()::String = string("Macintosh; ",mac_processor()," Mac OS X 10_",rand(5:8),"_",rand(0:9))

    function chrome()::String
        saf = string(rand(531:536),rand(0:2))
        platforms::Tuple = (linux_platform_token(),windows_platform_token(),mac_platform_token())

        return string("Mozilla/5.0 " ,"(",rand(platforms),") AppleWebKit/",saf," (KHTML, like Gecko) Chrome/",rand(13:70),".0.",rand(800:899),".0 Safari/",saf)
    end

    function firefox()::String
        ver::Tuple = (
                string("Gecko/",date_time_between()," Firefox/",rand(4:15),".0"),
                string("Gecko/",date_time_between()," Firefox/3.6.", rand(1:20)),
                string("Gecko/",date_time_between()," Firefox/3.8")
        )

        platforms::Tuple = (
            string("(",windows_platform_token(),"; ",rand(langs),"; rv:1.9.",rand(0:2),".20) ",rand(ver)),
            string("(",linux_platform_token(),"; rv:1.9.",rand(5:7),".20) ",rand(ver)),
            string("(",mac_platform_token(),"; rv:1.9.",rand(2:6),".20) ",rand(ver))
        )

        return string("Mozilla/5.0 " , rand(platforms))
    end

   function safari()::String
    saf::String = string(rand(531:535),".",rand(1:50),".",rand(1:7))
    ver::String=""

        if rand(0:1) == 0
            ver = string(rand(4:5),".",rand(0:1))
        else
            ver = string(rand(4:5),".0.",rand(1:5))
        end

        platforms::Tuple = (
        string("(Windows; U; ",windows_platform_token(),") AppleWebKit/",saf," (KHTML, like Gecko) Version/",ver," Safari/",saf),
        string("(",mac_platform_token()," rv:",rand(2:6),".0; ",rand(langs)," AppleWebKit/",saf," (KHTML, like Gecko) Version/",ver," Safari/",saf),
        string("(iPod; U; CPU iPhone OS ",rand(3:4),"_",rand(0:3)," like Mac OS X; ",rand(langs)," AppleWebKit/",saf," (KHTML, like Gecko) Version/",rand(3:4),".0.5 Mobile/8B",rand(111:119)," Safari/6",saf)
        );

        return string("Mozilla/5.0 ", rand(platforms))
    end

    function opera()::String

        platforms::Tuple = ( linux_platform_token(), windows_platform_token())

        return string("Opera/",rand(8:9),".",rand(10:99),".","(",rand(platforms),"; ",rand(langs)," Presto/2.9.",rand(160:190)," Version/",rand(10:12),".00")
    end

    internet_explorer()::String = string("Mozilla/5.0 (compatible; MSIE ",rand(5:9),".0; ",windows_platform_token(),"; Trident/",rand(3:5),".",rand(0:1))


text:
using ITensors
using ITensorInfiniteMPS

##############################################################################
# VUMPS parameters
#

maxdim = 30 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
outer_iters = 5 # Number of times to increase the bond dimension

model_params = (t=1.0, U=10.0, V=0.0)

##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#

N = 2 # Unit cell size

function electron_space_shift(q̃nf, q̃sz)
  return [
    QN(("Nf", 0 - q̃nf, -1), ("Sz", 0 - q̃sz)) => 1,
    QN(("Nf", 1 - q̃nf, -1), ("Sz", 1 - q̃sz)) => 1,
    QN(("Nf", 1 - q̃nf, -1), ("Sz", -1 - q̃sz)) => 1,
    QN(("Nf", 2 - q̃nf, -1), ("Sz", 0 - q̃sz)) => 1,
  ]
end

electron_space = fill(electron_space_shift(1, 0), N)
s = infsiteinds("Electron", N; space=electron_space)
initstate(n) = isodd(n) ? "↑" : "↓"
ψ = InfMPS(s, initstate)

model = Model"hubbard"()
@show model, model_params

# Form the Hamiltonian
H = InfiniteITensorSum(model, s; model_params...)

# Check translational invariance
println("\nCheck translational invariance of initial infinite MPS")
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

outputlevel = 1
vumps_kwargs = (tol=1e-8, maxiter=max_vumps_iters, outputlevel=outputlevel)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

# For now, to increase the bond dimension you must alternate
# between steps of VUMPS and subspace expansion (which outputs
# a new state that is equal to the original state but with
# a larger bond dimension)

println("\nRun VUMPS on initial product state, unit cell size $N")
ψ = vumps(H, ψ; vumps_kwargs...)

for _ in 1:outer_iters
  println("\nIncrease bond dimension")
  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
  println("Run VUMPS with new bond dimension")
  global ψ = vumps(H, ψ; vumps_kwargs...)
end

# Check translational invariance
println("\nCheck translational invariance of optimized infinite MPS")
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n)
  return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[]
end

function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2)
  n1, n2 = n1n2
  ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2]
  return (noprime(ϕ * h) * dag(ϕ))[]
end

function expect_two_site(ψ::MPS, h::ITensor, n1n2)
  n1, n2 = n1n2
  ψ = orthogonalize(ψ, n1)
  ϕ = ψ[n1] * ψ[n2]
  return (noprime(ϕ * h) * dag(ϕ))[]
end

Nup = [expect(ψ, "Nup", n) for n in 1:N]
Ndn = [expect(ψ, "Ndn", n) for n in 1:N]
Sz = [expect(ψ, "Sz", n) for n in 1:N]

bs = [(1, 2), (2, 3)]
energy_infinite = map(b -> expect_two_site(ψ, H[b], b), bs)

#
# Compare to DMRG
#

Nfinite = 100
sfinite = siteinds("Electron", Nfinite; conserve_qns=true)
Hfinite = MPO(model, sfinite; model_params...)
ψfinite = randomMPS(sfinite, initstate; linkdims=10)
println("\nQN sector of starting finite MPS")
@show flux(ψfinite)
sweeps = Sweeps(30)
maxdims =
  min.(maxdim, [2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 16, 16, 16, 16, 32, 32, 32, 32, 50])
@show maxdims
setmaxdim!(sweeps, maxdims...)
setcutoff!(sweeps, cutoff)
println("\nRun DMRG on $Nfinite sites")
energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
println("\nEnergy density")
@show energy_finite_total / Nfinite

nfinite = Nfinite ÷ 2 - 1
bsfinite = [(nfinite, nfinite + 1), (nfinite + 1, nfinite + 2)]
hfinite(b) = ITensor(model, sfinite[b[1]], sfinite[b[2]]; model_params...)
energy_finite = map(b -> expect_two_site(ψfinite, hfinite(b), b), bsfinite)

Nup_finite = ITensors.expect(ψfinite, "Nup")[nfinite:(nfinite + 1)]
Ndn_finite = ITensors.expect(ψfinite, "Ndn")[nfinite:(nfinite + 1)]
Sz_finite = ITensors.expect(ψfinite, "Sz")[nfinite:(nfinite + 1)]

energy_exact = reference(model, Observable("energy"); U=model_params.U / model_params.t)

println("\nResults from VUMPS")
@show energy_infinite
@show energy_exact
@show Nup
@show Ndn
@show Nup .+ Ndn
@show Sz

println("\nResults from DMRG")
@show energy_finite
@show Nup_finite
@show Ndn_finite
@show Nup_finite .+ Ndn_finite
@show Sz_finite

nothing


text:
function log_fma(x::Float64)
    if x > 0.0
        if isinf(x)
            Inf
        else
            xu = reinterpret(Uint64,x)
            k = int(xu >> 52) & 0x07ff
            if k == 0 # x is subnormal
                x *= 0x1p54 # normalise significand
                xu = reinterpret(Uint64,x)
                k = int(xu >> 52) & 0x07ff - 54
            end
            k -= 1022
            xu = (xu & 0x800f_ffff_ffff_ffff) | 0x3fe0_0000_0000_0000
            s = reinterpret(Float64,xu)
            if s < 0.7071067811865476
                s *= 2.0
                k -= 1
            end
            f = s-1.0
            fk = float64(k)

            ln2_hi = 6.93147180369123816490e-01
            ln2_lo = 1.90821492927058770002e-10

            lg1 = 6.666666666666735130e-01
            lg2 = 3.999999999940941908e-01
            lg3 = 2.857142874366239149e-01
            lg4 = 2.222219843214978396e-01
            lg5 = 1.818357216161805012e-01
            lg6 = 1.531383769920937332e-01
            lg7 = 1.479819860511658591e-01

            hfsq = 0.5*f*f
            s = f/(2.0+f)
            z = s*s
            w = z*z
            t1 = w*@horner_fma(w,lg2,lg4,lg6)
            t2 = @horner_fma(w,lg1,lg3,lg5,lg7)
            R = fma(t2,z,t1)
            fk*ln2_hi+(f-(hfsq-(s*(hfsq+R)+(fk*ln2_lo))))
        end
    elseif x == 0.0
        -Inf
    else
        throw(DomainError())
    end
end

function log_muladd(x::Float64)
    if x > 0.0
        if isinf(x)
            Inf
        else
            xu = reinterpret(Uint64,x)
            k = int(xu >> 52) & 0x07ff
            if k == 0 # x is subnormal
                x *= 0x1p54 # normalise significand
                xu = reinterpret(Uint64,x)
                k = int(xu >> 52) & 0x07ff - 54
            end
            k -= 1022
            xu = (xu & 0x800f_ffff_ffff_ffff) | 0x3fe0_0000_0000_0000
            s = reinterpret(Float64,xu)
            if s < 0.7071067811865476
                s *= 2.0
                k -= 1
            end
            f = s-1.0
            fk = float64(k)

            ln2_hi = 6.93147180369123816490e-01
            ln2_lo = 1.90821492927058770002e-10

            lg1 = 6.666666666666735130e-01
            lg2 = 3.999999999940941908e-01
            lg3 = 2.857142874366239149e-01
            lg4 = 2.222219843214978396e-01
            lg5 = 1.818357216161805012e-01
            lg6 = 1.531383769920937332e-01
            lg7 = 1.479819860511658591e-01

            hfsq = 0.5*f*f
            s = f/(2.0+f)
            z = s*s
            w = z*z
            t1 = w*@horner_muladd(w,lg2,lg4,lg6)
            t2 = @horner_muladd(w,lg1,lg3,lg5,lg7)
            R = muladd(t2,z,t1)
            fk*ln2_hi+(f-(hfsq-(s*(hfsq+R)+(fk*ln2_lo))))
        end
    elseif x == 0.0
        -Inf
    else
        throw(DomainError())
    end
end


text:
function _JStar(
    center,
    radius,
    color,
    linewidth,
    npoints,
    ratio,
    orientation,
    action,
    reversepath,
)
    sethue(color)
    setline(linewidth)
    star(center, radius, npoints, ratio, orientation, action, reversepath = reversepath)
    return center
end

"""
    1. JStar(center::Point, radius; kwargs...)
    2. JStar(xcenter, ycenter, radius; kwargs...)
        - same as 1. with `center = Point(xcenter, ycenter)`

Draw a star centered at a position.
Return the center of the star.

# Keywords for all
- `color` color of the outline or fill of the star (default: "black")
- `linewidth` linewidth of the outline (default: 1)
- `action` defines whether the rectangle should be outlined (`:stroke`) or filled (`:fill`)
- `npoints` number of points the star has (default: 5)
- `ratio` height of the smaller radius relative to the larger radius (default: 0.5)
- `orientation` orientation of the star given by an angle (default: 0)
- `reversepath` if true it reverses the path and therefore creates a hole (default: true)
"""
JStar(
    center::Point,
    radius;
    color = "black",
    linewidth = 1,
    npoints = 5,
    ratio = 0.5,
    orientation = 0,
    action = :stroke,
    reversepath = false,
) =
    (
        args...;
        center = center,
        radius = radius,
        color = color,
        linewidth = linewidth,
        action = action,
        npoints = npoints,
    ) -> _JStar(
        center,
        radius,
        color,
        linewidth,
        npoints,
        ratio,
        orientation,
        action,
        reversepath,
    )

JStar(xcenter, ycenter, radius; kwargs...) =
    JStar(Point(xcenter, ycenter), radius; kwargs...)


text:
import Oceananigans.TimeSteppers: calculate_tendencies!

using Oceananigans.Utils: work_layout
using Oceananigans: fields

using KernelAbstractions: @index, @kernel, Event, MultiEvent

using Oceananigans.Architectures: device

using Oceananigans.BoundaryConditions

"""
    calculate_tendencies!(model::ShallowWaterModel)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(model::ShallowWaterModel)

    # Note:
    #
    # "tendencies" is a NamedTuple of OffsetArrays corresponding to the tendency data for use
    # in GPU computations.
    #
    # "model.timestepper.Gⁿ" is a NamedTuple of Fields, whose data also corresponds to
    # tendency data.

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_interior_tendency_contributions!(model.timestepper.Gⁿ,
                                               model.architecture,
                                               model.grid,
                                               model.gravitational_acceleration,
                                               model.advection,
                                               model.coriolis,
                                               model.closure,
                                               model.bathymetry,
                                               model.solution,
                                               model.tracers,
                                               model.diffusivity_fields,
                                               model.forcing,
                                               model.clock,
                                               model.formulation)

    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_boundary_tendency_contributions!(model.timestepper.Gⁿ,
                                               model.architecture,
                                               model.solution,
                                               model.tracers,
                                               model.clock,
                                               fields(model))

    return nothing
end

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_tendency_contributions!(tendencies,
                                                    arch,
                                                    grid,
                                                    gravitational_acceleration,
                                                    advection,
                                                    coriolis,
                                                    closure,
                                                    bathymetry,
                                                    solution,
                                                    tracers,
                                                    diffusivities,
                                                    forcings,
                                                    clock,
                                                    formulation)

    workgroup, worksize = work_layout(grid, :xyz)

    calculate_Guh_kernel! = calculate_Guh!(device(arch), workgroup, worksize)
    calculate_Gvh_kernel! = calculate_Gvh!(device(arch), workgroup, worksize)
    calculate_Gh_kernel!  =  calculate_Gh!(device(arch), workgroup, worksize)
    calculate_Gc_kernel!  =  calculate_Gc!(device(arch), workgroup, worksize)

    barrier = Event(device(arch))

    args_vel = (grid, gravitational_acceleration, advection.momentum, coriolis, closure,
                      bathymetry, solution, tracers, diffusivities, forcings, clock, formulation)
    args_h   = (grid, gravitational_acceleration, advection.mass, coriolis, closure,
                      solution, tracers, diffusivities, forcings, clock, formulation)

    Guh_event = calculate_Guh_kernel!(tendencies[1], args_vel...; dependencies = barrier)
    Gvh_event = calculate_Gvh_kernel!(tendencies[2], args_vel...; dependencies = barrier)
    Gh_event  =  calculate_Gh_kernel!(tendencies[3], args_h...;   dependencies = barrier)

    events = [Guh_event, Gvh_event, Gh_event]

    for (tracer_index, tracer_name) in enumerate(propertynames(tracers))
        @inbounds c_tendency = tendencies[tracer_index+3]
        @inbounds forcing = forcings[tracer_index+3]
        @inbounds c_advection = advection[tracer_name]

        Gc_event = calculate_Gc_kernel!(c_tendency, grid, Val(tracer_index), c_advection, closure, solution,
                                        tracers, diffusivities, forcing, clock, formulation, dependencies=barrier)

        push!(events, Gc_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Tendency calculators for the transports and height: uh, vh, h
#####

""" Calculate the right-hand-side of the uh-transport equation. """
@kernel function calculate_Guh!(Guh,
                                grid,
                                gravitational_acceleration,
                                advection,
                                coriolis,
                                closure,
                                bathymetry,
                                solution,
                                tracers,
                                diffusivities,
                                forcings,
                                clock,
                                formulation)

    i, j, k = @index(Global, NTuple)

    @inbounds Guh[i, j, k] = uh_solution_tendency(i, j, k, grid, gravitational_acceleration, advection, coriolis, closure,
                                                    bathymetry, solution, tracers, diffusivities, forcings, clock, formulation)
end

""" Calculate the right-hand-side of the vh-transport equation. """
@kernel function calculate_Gvh!(Gvh,
                                grid,
                                gravitational_acceleration,
                                advection,
                                coriolis,
                                closure,
                                bathymetry,
                                solution,
                                tracers,
                                diffusivities,
                                forcings,
                                clock,
                                formulation)

    i, j, k = @index(Global, NTuple)

    @inbounds Gvh[i, j, k] = vh_solution_tendency(i, j, k, grid, gravitational_acceleration, advection, coriolis, closure,
                                                    bathymetry, solution, tracers, diffusivities, forcings, clock, formulation)
end

""" Calculate the right-hand-side of the height equation. """
@kernel function calculate_Gh!(Gh,
                               grid,
                               gravitational_acceleration,
                               advection,
                               coriolis,
                               closure,
                               solution,
                               tracers,
                               diffusivities,
                               forcings,
                               clock,
                               formulation)

    i, j, k = @index(Global, NTuple)

    @inbounds Gh[i, j, k] = h_solution_tendency(i, j, k, grid, gravitational_acceleration, advection, coriolis, closure,
                                                solution, tracers, diffusivities, forcings, clock, formulation)
end

#####
##### Tracer(s)
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
@kernel function calculate_Gc!(Gc,
                               grid,
                               tracer_index,
                               advection,
                               closure,
                               solution,
                               tracers,
                               diffusivities,
                               forcing,
                               clock,
                               formulation)

    i, j, k = @index(Global, NTuple)

    @inbounds Gc[i, j, k] = tracer_tendency(i, j, k, grid, tracer_index, advection, closure, solution, tracers,
                                            diffusivities, forcing, clock, formulation)
end

#####
##### Boundary contributions to tendencies due to user-prescribed fluxes
#####

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_tendency_contributions!(Gⁿ, arch, solution, tracers, clock, model_fields)

    barrier = Event(device(arch))

    prognostic_fields = merge(solution, tracers)

    events = []

    # Solution fields and tracer fields
    for i in 1:length(Gⁿ)
        x_bcs_event = apply_x_bcs!(Gⁿ[i], prognostic_fields[i], arch, barrier, clock, model_fields)
        y_bcs_event = apply_y_bcs!(Gⁿ[i], prognostic_fields[i], arch, barrier, clock, model_fields)
        push!(events, x_bcs_event, y_bcs_event)
    end

    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end


text:
norm(a::SVector{N,T}) where{N,T} = sqrt(sum(a.^2))

const GEO_FAC = 4.0 / 3.0 * pi
const KERNELCONST = 16.0 / pi
@inline ramp(x) = max(0, x);
@inline function kernel_cubic(x::T) where {T}
    return T(KERNELCONST * (ramp(1.0 - x)^3 - 4.0 * ramp(0.5 - x)^3))
end


text:
  # Clear
  const p = 0
  const q = 0
  const r = 0
  const θ7 = 10
  const A1s = 0
  const B1s = 0
  const ρ = 1.225
  const G = 280
  const m = 280
  const Ix = 26
  const Iy = 113
  const Iz = 89
  const Ixy = -1.8
  const Iyz = 0.1
  const Izx = 4.3
  const g = 9.8
  const x = 1.64
  const y = 0
  const z = 0.42

  const R = 2.54
  const Ω = 76.8 # rad/s
  # Ω = 734rpm
  const Ib = 7.806
  const Mb = 5.588
  const a00 = 7.19
  const alpha00 = -0.75
  const b = 0.178
  const n = 2
  const eroot = 0
  const deltw = 0
  const theta1 = -10
  const sigma = 0.045
  const alpha00 = -0.75
  const X_mr = x-1.65
  const Y_mr = 0-y
  const Z_mr = z-1.24
  # 尾桨原始参数
  const theta1_T = 0
  const n = 2
  const ΩT = 470
  const sigmaT = 0.11
  # 4488r/min
  const RT = 0.38
  const a00T = 6.05
  const alpha00T = 0
  const bT = 0.068
  const K_QT = 1
  const X_T = x-4.65
  const Y_T = -0.12-y
  const Z_T = z-0.473
  #平尾原始参数
  const K_MH = 1
  const K_QH = 1
  # v_01 = v_0
  const alpha_sf = 0
  const phi_H = -4
  const A_H = 0.3
  const X_H = x-3.81
  const Y_H = 0
  const Z_H = z-0.4073
  #垂尾原始参数
  #机身原始参数
  const C_DF = 0.2
  const C_LF = 0.02
  const C_SF = 0.06
  const C_MxF = -0.001
  const C_MyF = 0.001
  const C_MzF = -0.01
  const A_F = 0.72
  const l_F = 4.908
  const X_F = 0
  const Y_F = 0
  const Z_F = 0

  Trzb = [0 Z_mr -Y_mr; -Z_mr 0 X_mr; Y_mr -X_mr 0]
  u = [20.0,0.0,0.0]+Trzb*[p,q,r]
  Vx = u[1]
  Vy = u[2]
  Vz = u[3]
  q1 = 1/2*ρ*π*R^2*(Ω*R)^2
  μ = sqrt(Vx*Vx+Vy*Vy)/Ω/R
  lamda00 = -Vz/Ω/R
  theta0 = θ7+7
  if abs(Vx)<=1e-3&&abs(Vy)<=1e-3
      alphas = π/2
  else
      alphas = atan(Vz/sqrt(Vx*Vx+Vy*Vy))
  end

  lamda_b = (ρ*a00*b*R^4)/Ib
  vO0 = 0.000001
  vOK = 1
  vOT = 0.000001
  eM = 0.92*a00*sigma*(theta0/57.3*(1/3+μ^2/2)+
                          theta1/57.3*(1+μ^2)/4-
                            (vOT-lamda00)/2-μ*B1s/57.3/2)       # CT旋翼拉力系数
  fM = 4*sqrt(μ^2+(vOT-lamda00)^2)
  FvOT = vOT-eM/fM
  i = 0
  while abs(FvOT)>0.000001
      if i != 0
          if FvOT>0
              vOK = vOT
          else
              vO0 = vOT
          end
      end
      i = 1
      aM = 0.92*a00*sigma*(theta0/57.3*
                            (1/3+μ^2/2)+theta1/57.3*
                              (1+μ^2)/4-(vO0-lamda00)/2-μ*B1s/57.3/2)
      bM = 4*sqrt(μ^2+(vO0-lamda00)^2)
      cM = 0.92*a00*sigma*(theta0/57.3*
                            (1/3+μ^2/2)+theta1/57.3*
                              (1+μ^2)/4-(vOK-lamda00)/2-μ*B1s/57.3/2)
      dM = 4*sqrt(μ^2+(vOK-lamda00)^2)
      FvO0 = vO0-aM/bM
      FvOK = vOK-cM/dM
      vOT = vOK-FvOK*(vOK-vO0)/(FvOK- FvO0)
      eM = 0.92*a00*sigma*(theta0/57.3*
                            (1/3+μ^2/2)+theta1/57.3*
                              (1+μ^2)/4-(vOT-lamda00)/2-μ*B1s/57.3/2)
      fM = 4*sqrt(μ^2+(vOT-lamda00)^2)
      FvOT = vOT-eM/fM
  end

  v_0 = vOT
  v_01 = v_0*Ω*R
  C_T = eM
  beta0 = 3/8*lamda_b*C_T/a00/sigma-1.5*9.8*R/((Ω*R)^2)
  beta00 = beta0*57.3
  a1 = (2*μ/(1-0.5*μ*μ)*
        (4/3*theta0/57.3+theta1/57.3+μ*alphas-v_0/Ω/R)-
          (1+1.5*μ*μ)*B1s/57.3/(1-0.5*μ*μ))
  a11 = a1*57.3
  b1 = 1/(1+0.5*μ*μ)*(4/3*μ*beta0+v_0/Ω/R)+A1s/57.3
  b11 = b1*57.3
  xt = R
  xs = 0.1*R
  rstep = (xt-xs)/20
  r1 = Array{Float64}(21)
  for j in 1:21
    r1[j] = xs+(j-1)*rstep
  end
  # r1 = xs:rstep:xt
  ψ = linspace(0,2*π,12)
  Ts1 = 0
  Hs1 = 0
  Ss1 = 0
  M1 = 0
  P1 = 0
  Cd1 = Array{Float64}(11,20)
  DQ1 = Array{Float64}(11,20)
  DSs1 = Array{Float64}(11,20)
  Ss11 = Array{Float64}(11)
  Hs11 = Array{Float64}(11)
  for i = 1:11
      Ts = 0
      Hs = 0
      Ss = 0
      M = 0
      P = 0
      sp = sin(ψ[i])
      cp = cos(ψ[i])
      beta = beta0-a1*cp-b1*sp
      beta_d = Ω*(a1*sp-b1*cp)
      beta_dd = Ω*Ω*(a1*cp+b1*sp)
      for j = 1:20
          rad = r1[j]
          theta = θ7+(rad/R-0.7)*theta1-A1s*cp-B1s*sp
          theta = theta/57.3
          ut = Ω*R*(rad/R+μ*sp)
          up = -(Ω*R*(-lamda00-μ*beta*cp)-v_0*Ω*R*(1+rad/R*cp)-rad*beta_d)
          alpha = atan(up/ut)
          alpha11 = theta-alpha-alpha00/57.3
  #         mach = sqrt(up^2+ut^2)/340
  #         CLCD0 = CLCD(alpha11,mach)
  #         CL0 = CLCD0 [1]
  #         CD0 = CLCD0[2]
          CL = a00*alpha11
          Cd = 0.008-0.003*CL+0.01*CL*CL
          Cd1[i,j] = Cd
          q2 = 0.5*ρ*(up^2+ut^2)
          Dy = q2*CL*b
          Dx = q2*Cd*b
          DT = Dy*cos(alpha)-Dx*sin(alpha)
          DQ = Dx*cos(alpha)+Dy*sin(alpha)
          DQ1[i,j] = DQ
          DTs = DT*cos(beta)
          DHs = DQ*sin(ψ[i])-DT*sin(beta)*cos(ψ[i])
          DSs = -DQ*cos(ψ[i])-DT*sin(beta)*sin(ψ[i])
          DSs1[i,j] = DSs
          DM = DQ*rad
          DP = DQ*rad*Ω
          Ts = Ts+DTs*rstep
          Hs = Hs+DHs*rstep
          Ss = Ss+DSs*rstep
          M = M+DM*rstep
          P = P+DP*rstep
      end
      Ts1 = Ts1+Ts
      Hs1 = Hs1+Hs
      Ss1 = Ss1+Ss
      Ss11[i] = Ss
      Hs11[i] = Hs
      M1 = M1+M
      P1 = P1+P
  end
  Ts = Ts1*n/11
  Hs = Hs1*n/11
  Ss = Ss1*n/11
  M = M1*n/11
  P = P1*n/11
  Mx = 0.5*n*eroot*Mb*Ω^2*b1
  My = 0.5*n*eroot*Mb*Ω^2*a1
  Trzb1 = [cos(deltw) -sin(deltw) 0; 0 0 1; sin(deltw) cos(deltw) 0]
  #旋翼侧滑角
  betaw = asin(Vy/sqrt(Vx*Vx+Vy*Vy))
  Ts = Ts
  Hs = Hs*cos(betaw)+Ss*sin(betaw)
  Ss = -Hs*sin(betaw)+Ss*cos(betaw)
  # v_0 = v_0*Ω*R
  beta0 = beta0*57.3
  a1 = a1*57.3
  b1 = b1*57.3
  # C_T = C_T0
  F = Trzb1*[-Hs,-Ts,-Ss]
  M = Trzb1*[Mx,M,My]+-1*Trzb*F
  FM = [F,M]
  FX = FM
  Cd1 = Cd1
  DQ1 = DQ1
  DSs1 = DSs1
  Ss11 = Ss11
  Hs11 = Hs11

max_stars_repo_path: src/[For Validation]matlab_rotor.jl

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_path: src/[For Validation]matlab_rotor.jl

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_path: src/[For Validation]matlab_rotor.jl

max_forks_repo_name: helicopterliu/RotorTrim


text:
#!/usr/bin/env julia

using LightXML
using OrderedCollections

function pdb_standard_bonds(residues, pdb_names)
    pdb_names_root = root(parse_file(pdb_names))
    glossary = OrderedDict()
    alternatives = OrderedDict()
    for element in pdb_names_root["Residue"]
        attributes = attributes_dict(element)
        atoms_dict = haskey(attributes, "type") ? copy(glossary[attributes["type"]]) : OrderedDict()
        for item in element["Atom"]
            list = join(values(attributes_dict(item)), "|")
            atoms_dict[attribute(item, "name")] = "\\b($(list))\\b"
        end
        glossary[attributes["name"]] = atoms_dict
        alt_names = []
        while haskey(attributes, "alt$(length(alt_names) + 1)")
            push!(alt_names, attributes["alt$(length(alt_names) + 1)"])
        end
        alternatives[attributes["name"]] = alt_names
    end

    residues_root = root(parse_file(residues))
    residue_dict = OrderedDict{String, Vector{Vector{String}}}()
    for element in residues_root["Residue"]
        attributes = attributes_dict(element)
        name = attributes["name"]
        residue_dict[name] = [map(x->get(glossary[name], x, "\\b($x)\\b"),
                              values(attributes_dict(item)))
                              for item in element["Bond"]]
        for alt in alternatives[name]
            residue_dict[alt] = residue_dict[name]
        end
    end
    return residue_dict
end

const ELEMENTS = LittleDict(
    "H" => 1.008,
    "C" => 12.011,
    "N" => 14.007,
    "O" => 15.999,
    "P" => 30.973762,
    "S" => 32.06,
)

sanitized(str) = replace(replace(replace(str, "-" => "_"), "'" => "p"), "*" => "a")

function pdb_aliases_xml(pdb_bonds)
    regex_set = OrderedSet(regex for list in values(pdb_bonds) for bond in list for regex in bond)

    xdoc = XMLDocument()
    xroot = create_root(xdoc, "Residues")
    elements = new_child(xroot, "Elements")
    for (element, mass) in ELEMENTS
        set_attributes(
            new_child(elements, "Element"),
            Dict("name" => element, "mass" => mass),
        )
    end

    regular_expressions = new_child(xroot, "RegularExpressions")
    regex_id = Dict()
    for (id, regex) in enumerate(regex_set)
        regex_entry = new_child(regular_expressions, "Regex")
        set_attributes(regex_entry, Dict("id" => id, "code" => sanitized(regex)))
        regex_id[regex] = id
    end

    for (resname, list) in pdb_bonds
        residue = new_child(xroot, "Residue")
        set_attributes(residue, Dict("name" => resname))
        for bond in list
            ids = Dict("id$i" => regex_id[regex] for (i, regex) in enumerate(bond))
            set_attributes(new_child(residue, "Bond"), ids)
        end
    end
    return xdoc
end

url = "raw.githubusercontent.com/openmm/openmm/master/wrappers/python/simtk/openmm/app/data/"
files = ["residues.xml", "pdbNames.xml"]
pipes = []
for file in files
    push!(pipes, download(url*file))
end

pdb_bonds = pdb_standard_bonds(broadcast(joinpath, @__DIR__, pipes)...)
pdb_aliases_content = pdb_aliases_xml(pdb_bonds)

io = open(joinpath(@__DIR__, "..", "src", "data", "pdb_aliases.xml"), "w")
write(io, string(pdb_aliases_content))
close(io)


text:
"""
A module which provides types, and transformations
for provided and user-provided orbital frames, e.g.
Earth-Centered-Inertial, Heliocentric-Inertial,
Synodic.

# Extended help

**Exports**

$(EXPORTS)

**Imports**

$(IMPORTS)
"""
module OrbitalFrames

# Frames
export OrbitalFrame, Inertial, Rotating
export BodycentricInertial, BarycentricInertial
export BodycentricRotating, BarycentricRotating
export isinertial, isrotating, isbodycentric, isbarycentric
export @frame

# Transforms
export AstrodynamicsTransform
export Transform

using Unitful
using DocStringExtensions
using CoordinateTransformations

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("Frames.jl")
include("Transforms.jl")

end # module


text:
module Funds

    using Dates
    import ..YearMonths
    import ..HTTP

    struct FundData
        fund_type::Union{String, Nothing}
        cnpj::String
        date::Date
        total_value::Float64
        quota_value::Float64
        net_worth::Float64
        deposits::Float64
        withdraws::Float64
        number_quota_holders::Int64
    end

    function get_month_quotes(ym::YearMonths.YearMonth)::Vector{FundData}
        @assert ym >= YearMonths.YearMonth(2017, 1)

        raw = _get_data(ym)

        out = Vector{FundData}()

        header = split(String(raw[1]), ";")

        for i in 2:length(raw)

            line = raw[i]

            if line == ""
                continue
            end

            s = split(String(line), ";")

            if ym <= YearMonths.YearMonth(2020,3)

                @assert length(s) == 8 "bad CVM data, this line should have 8 fields: $s"
                fund_type = nothing
                cnpj = String(s[1])
                date = Date(String(s[2]))
                total_value = parse(Float64, String(s[3]))
                quota_value = parse(Float64, String(s[4]))
                net_worth = parse(Float64, String(s[5]))
                deposits = parse(Float64, String(s[6]))
                withdraws = parse(Float64, String(s[7]))
                number_quota_holders = parse(Int64, String(s[8]))

            else
                # after 2020-03
                @assert length(s) == 9 "bad CVM data, this line should have 9 fields: $s"
                fund_type = String(s[1])
                cnpj = String(s[2])
                date = Date(String(s[3]))
                total_value = parse(Float64, String(s[4]))
                quota_value = parse(Float64, String(s[5]))
                net_worth = parse(Float64, String(s[6]))
                deposits = parse(Float64, String(s[7]))
                withdraws = parse(Float64, String(s[8]))
                number_quota_holders = parse(Int64, String(s[9]))

            end
            d = FundData(
                fund_type, cnpj, date,
                total_value, quota_value, net_worth,
                deposits, withdraws, number_quota_holders
            )
            push!(out, d)

        end

        return(out)
    end

    function _get_url(ym::YearMonths.YearMonth)::String
        str_year = "$(YearMonths.year(ym))"
        str_month = lpad("$(YearMonths.month(ym))", 2, "0")
        url = "http://dados.cvm.gov.br/dados/FI/DOC/INF_DIARIO/DADOS/inf_diario_fi_$str_year$str_month.csv"
        return url
    end

    function _get_data(ym::YearMonths.YearMonth)
        url = _get_url(ym)
        r = HTTP.get(url, status_exception = false)
        if r.status == 200
            return split(String(r.body), "\r\n")
        elseif r.status == 404
            throw(BoundsError("CVM fund data not found for $ym"))
        else
            error("error retrieving CVM fund data: $r")
        end
    end

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/subsolver_config/config/clp.jl
# Configuration adjustment subroutines for Clp.
#############################################################################

function set_default_config!(::Val{:clp}, ext::ExtensionType, m::T) where T <: MOI.AbstractOptimizer

    MOI.set(m, MOI.RawParameter("PrimalTolerance"),       1E-7)
    MOI.set(m, MOI.RawParameter("DualTolerance"),         1E-7)
    MOI.set(m, MOI.RawParameter("DualObjectiveLimit"),    1e308)
    MOI.set(m, MOI.RawParameter("MaximumIterations"),     2147483647)
    MOI.set(m, MOI.RawParameter("PresolveType"),          0)
    MOI.set(m, MOI.RawParameter("SolveType"),             5)
    MOI.set(m, MOI.RawParameter("InfeasibleReturn"),      1)
    MOI.set(m, MOI.RawParameter("Scaling"),               3)
    MOI.set(m, MOI.RawParameter("Perturbation"),          100)

    return nothing
end