text:
using LowRankModels, Random
import StatsBase: sample, Weights, mean
import LinearAlgebra: norm

# test MNL Ordinal loss

## generate data
Random.seed!(1);
m,n,k = 100,100,2;
kfit = k+1
nlevels = 5; # number of levels
d = nlevels-1 # embedding dimension
D = n*d;
# coordinates of covariates
X_real = randn(m,k)
# directions of observations
Y_real = randn(k,n)
# measurement thresholds
T_real = k*randn(d,n) # notice x^T y has variance k; so making this bigger makes the problem easier
for j=1:n # this scheme doesn't work to ensure uniform sampling
	T_real[:,j] = sort(T_real[:,j])
end

signedsums = Array{Float64}(undef, d, nlevels)
for i=1:d
    for j=1:nlevels
        signedsums[i,j] = i<j ? 1 : -1
    end
end

XY = X_real*Y_real;
XYplusT = zeros(Float64, (m,D))
A = zeros(Int, (m, n))
for i=1:m
	for j=1:n
		u = XY[i,j] .+ T_real[:,j]
		XYplusT[i,(j-1)*d .+ (1:d)] = u
		diffs = u'*signedsums
		wv = Weights(Float64[exp(-diffs[l]) for l in 1:nlevels])
		l = sample(wv)
		A[i,j] = l
	end
end
# loss is insensitive to shifts; regularizer should pick this shift
XYplusT = XYplusT .- mean(XYplusT, dims=2)

# and the model
losses = fill(MultinomialOrdinalLoss(nlevels),n)
rx, ry = lastentry1(QuadReg(.01)), MNLOrdinalReg(QuadReg(.01)) #lastentry_unpenalized(QuadReg(10));
Yord = randn(kfit,D)
yidxs = get_yidxs(losses)
for j=1:n
	prox!(ry, view(Yord,:,yidxs[j]), 1)
end
glrm = GLRM(A,losses,rx,ry,kfit, scale=false, offset=false, X=randn(kfit,m), Y=Yord);

# fit w/o initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y#[:,1:d:D];
@show ch.objective
println("After fitting, parameters differ from true parameters by $(norm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs.(A_imputed - A)) / prod(size(A))*100)% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)% of entries would be wrong.)\n")

# initialize
init_svd!(glrm)
XYh = glrm.X' * glrm.Y
println("After initialization with the svd, parameters differ from true parameters by $(norm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After initialization with the svd, $(sum(A_imputed .!= A) / prod(size(A))*100)% of imputed entries are wrong")
println("After initialization with the svd, imputed entries are off by $(sum(abs,(A_imputed - A)) / prod(size(A))*100)% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)% of entries would be wrong.)\n")

# fit w/ initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y#[:,1:d:D];
@show ch.objective
println("After fitting, parameters differ from true parameters by $(norm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs,(A_imputed - A)) / prod(size(A))*100)% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)% of entries would be wrong.)\n")

# test scaling
mul!(glrm.ry, 3)


text:
#####################################################
# Surya
# NREL
# November 2021
# Call SIIP PSY Tabular Data Parser
#####################################################################################
const PSY = PowerSystems
#####################################################################################
# Main Function
#####################################################################################
function parse_tabular_data(csv_dir::String,base_MVA::Float64,rt_flag::Bool;ts_pointers_file::Union{Nothing, String} = nothing)

    dir_name = @__DIR__
    user_descriptors_file = joinpath(dir_name,"Descriptors","user_descriptors.yaml")
    generator_mapping_file = joinpath(dir_name,"Descriptors","generator_mapping.yaml")

    if (ts_pointers_file === nothing)
        @warn "Time series pointers file type wasn't passed. Using timeseries_pointers.csv"
        ts_pointers_file = "CSV"
    end

    if (ts_pointers_file == "CSV")
        @info "Using timeseries_pointers.csv as time series metadata file...."
        ts_pointers_file= "timeseries_pointers.csv"
    elseif (ts_pointers_file == "JSON")
        @info "Using timeseries_pointers.json as time series metadata file...."
        ts_pointers_file= "timeseries_pointers.json"
    else
        error("Unrecognized time series pointers file type")
    end
    timeseries_pointers_file = joinpath(csv_dir, ts_pointers_file)

    rawsys = PSY.PowerSystemTableData(
        csv_dir,
        base_MVA,
        user_descriptors_file,
        timeseries_metadata_file = timeseries_pointers_file,
        generator_mapping_file = generator_mapping_file,
    );

    sys_DA = PSY.System(rawsys; time_series_resolution = Dates.Hour(1));

    if (rt_flag)
        sys_RT = PSY.System(rawsys; time_series_resolution = Dates.Minute(5));
        @info "Successfully generated both DA and RT PSY Systems."
        return sys_DA, sys_RT
    else
        @info "Successfully generated DA PSY System."
        return sys_DA
    end
end

ext: jl


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Convert a satellite state vector from an Earth-Centered, Earth-Fixed (ECEF)
#   reference frame to a Earth-Centered Inertial (ECI) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export svECEFtoECI

"""
    svECEFtoECI(sv::SatelliteStateVector, ECEF, ECI, JD_UTC [, eop_data])

Convert the satellite state vector `sv` from the Earth-Centered, Earth-Fixed
(ECEF) reference frame `ECEF` to the Earth-Centered Inertial (ECI) reference
frame at the Julian day `JD_UTC` [UTC]. The `eop_data` may be required depending
on the selection of the input and output reference system. For more information,
see the documentation of the function `rECEFtoECI`.

!!! info

    It is assumed that the input velocity and acceleration in `sv` are obtained
    by an observer on the ECEF frame. Thus, the output will contain the velocity
    and acceleration as measured by an observer on the ECI frame.

"""
function svECEFtoECI(sv::SatelliteStateVector, T_ECEF::Val{:ITRF},
                     T_ECI::T_ECIs, JD_UTC::Number, eop_data::EOPData_IAU1980)

    # First, convert from the ITRF to PEF.
    sv_pef = svECEFtoECEF(sv, Val(:ITRF), Val(:PEF), JD_UTC, eop_data)

    # Now, convert from PEF to ECI.
    return svECEFtoECI(sv_pef, Val(:PEF), T_ECI, JD_UTC, eop_data)
end

function svECEFtoECI(sv::SatelliteStateVector, T_ECEF::Val{:ITRF},
                     T_ECI::T_ECIs_IAU_2006, JD_UTC::Number,
                     eop_data::EOPData_IAU2000A)

    # First, convert from the ITRF to TIRS.
    sv_tirs = svECEFtoECEF(sv, Val(:ITRF), Val(:TIRS), JD_UTC, eop_data)

    # Now, convert from TIRS to ECI.
    return svECEFtoECI(sv_tirs, Val(:TIRS), T_ECI, JD_UTC, eop_data)
end

function svECEFtoECI(sv::SatelliteStateVector,
                     T_ECEF::Union{Val{:PEF},Val{:TIRS}},
                     T_ECI::Union{T_ECIs, T_ECIs_IAU_2006}, JD_UTC::Number,
                     eop_data::Union{Nothing, EOPData_IAU1980, EOPData_IAU2000A} = nothing)

    # Get the matrix that converts the ECEF to the ECI.
    if eop_data == nothing
        D = rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC)
    else
        D = rECEFtoECI(DCM, T_ECEF, T_ECI, JD_UTC, eop_data)
    end

    # Since the ECI and ECEF frames have a relative velocity between them, then
    # we must account from it when converting the velocity and acceleration. The
    # angular velocity between those frames is computed using `we` and corrected
    # by the length of day (LOD) parameter of the EOP data, if available.
    ω  = we * ( 1 - (eop_data != nothing ? eop_data.LOD(JD_UTC)/86400 : 0 ) )
    vω = SVector{3}(0,0,ω)

    # Compute the position in the ECI frame.
    r_eci = D*sv.r

    # Compute the velocity in the ECI frame.
    vω_x_r = vω × sv.r
    v_eci = D*(sv.v + vω_x_r )

    # Compute the acceleration in the ECI frame.
    a_eci = D*(sv.a + vω × vω_x_r + 2vω × sv.v)

    return satsv(sv.t, r_eci, v_eci, a_eci)
end


text:
export STCG

using LinearAlgebra
using CUTEst
using NLPModels, LinearOperators, Krylov, SolverTools, SolverBenchmark,JuMP, Ipopt

function STCG(nlp;x :: AbstractVector=copy(nlp.meta.x0),
    atol :: Real=√eps(eltype(x)), rtol :: Real=√eps(eltype(x)),max_iter = 1_000_000, max_time = 30.0)
    t₀ = time()
    Δt = time() - t₀

    f(x) = obj(nlp,x)
    ∇f(x) = grad(nlp,x)
    gx = ∇f(x)
    H(x) = hess(nlp,x)
    Δ = 1.0
    η = 1e-3
    iter = 0
    ϵ = atol + rtol * norm(gx)

    status =:unknown

    p= x .* 0
    z = zero(∇f(x))
    while norm(∇f(x)) > ϵ

        r = ∇f(x)
        d = -r
        n = length(∇f(x))
        fx = f(x)
        gx = ∇f(x)
        Hx = H(x)

        if norm(r) < ϵ
            p = z
        end

        for j = 1:n

            if dot(d, Hx*d) ≤ 0
                model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
                @variable(model, t)
                @objective(model, Min, fx + dot(gx, z + t*d) + dot((z + t*d), Hx*(z + t*d))/2)
                for i = 1:length(z)
                    @constraint(model, (z[i] + t*d[i])^2 == Δ^2)
                end
                @constraint(model, t ≥ 0)
                optimize!(model)
                τ = value.(t)
                p = z + τ*d
            end
            α = dot(r, r)/dot(d, Hx*d)
            znew = z + α*d

            #end 1st test

            if norm((znew)) ≥ Δ
                model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
                @variable(model, t)
                for i = 1:length(z)
                    @objective(model, Min, z[i] + t*d[i])
                    @constraint(model, (z[i] + t*d[i])^2 == Δ^2)
                end
                @constraint(model, t ≥ 0)
                optimize!(model)
                τ = value.(t)
                p = z + τ*d
            end
            rnew = r + α*Hx*d

            #end 2nd test

            if norm((rnew)) < ϵ
                p = znew
            end

            βnew = dot(rnew, rnew)/dot(r, r)
            dnew = -rnew + βnew*d

            #end 3rd test
            if Δt > max_time
                status = :max_time
                break
            end

            if iter > max_iter
                status = :max_iter
                break
            end
            z, r, d = znew, rnew, dnew
        end
        Δt = time() - t₀

        ared = f(x) - f(x + p)
        pred = f(x) - (fx + dot(gx, p) + dot(p, Hx*p)/2)

        ρ = ared/pred

        if norm(ρ) < 0.25
            Δ = 0.25*Δ
        elseif norm(ρ) > 0.75 && norm(ρ) == Δ
            Δ = min(2*Δ, ϵ)
        else
             Δ = Δ
        end

        if norm(ρ) > η
            x = x + p
        else
            x = x
        end

        Δt = time() - t₀
        iter += 1

        if iter > max_iter
            status = :max_iter
            break
        end

        if Δt > max_time
            status = :max_time
            break
        end

    end

    if norm(∇f(x)) <= ϵ
        status = :first_order
    elseif (Δt > max_time || iter > max_iter)
        if Δt > max_time
            status = :max_time
        elseif iter > max_iter
            status = :max_iter
        end
    end
    return GenericExecutionStats(status,nlp,objective=f(x),solution=x,dual_feas=norm(∇f(x)),iter = iter, elapsed_time = Δt)
end


text:
# Detect Elemental integer size
function ElIntType()
    using64 = Cint[0]
    err = ccall((:ElUsing64BitInt, libEl), Cuint, (Ptr{Cint},), using64)
    err == 0 || throw(ElError(err))
    return using64[1] == 1 ? Int64 : Int32
end
const ElInt = ElIntType()

function ElCommType()
    sameSizeAsInt = Cint[0]
    err = ccall((:ElMPICommSameSizeAsInteger, libEl), Cuint, (Ptr{Cint},),
      sameSizeAsInt)
    err == 0 || throw(ElError(err))
    return sameSizeAsInt[1] == 1 ? Cint : Ptr{Void}
end
const ElComm = ElCommType()

function ElGroupType()
    sameSizeAsInt = Cint[0]
    err = ccall((:ElMPIGroupSameSizeAsInteger, libEl), Cuint, (Ptr{Cint},),
      sameSizeAsInt)
    err == 0 || throw(ElError(err))
    return sameSizeAsInt[1] == 1 ? Cint : Ptr{Void}
end
const ElGroup = ElGroupType()

# Detect Elemental Bool type
function ElBoolType()
    # NOTE: Returning Uint8 when C claims that sizeof(bool) is 1 byte leads
    #       to improperly passed structs to Elemental's C interface. This is
    #       worth investigating and might be an alignment issue.

    # warn("Hardcoding ElBool to Cint")
    boolsize = Ref(zero(Cuint))
    err = ccall((:ElSizeOfCBool, libEl), Cuint, (Ref{Cuint},), boolsize)
    err == 0 || throw(ElError(err))
    return boolsize[] == 1 ? UInt8 : Cint
end
const ElBool = ElBoolType()

function ElBool(value::Bool)
    if value
      return ElBool(1)
    else
      return ElBool(0)
    end
end

using Base.LinAlg: BlasFloat, BlasReal, BlasComplex

typealias ElElementType Union{ElInt,Float32,Float64,Complex64,Complex128}

typealias ElFloatType Union{Float32,Float64} # TODO: Maybe just use BlasReal here

abstract ElementalMatrix{T} <: AbstractMatrix{T}
eltype{T}(A::ElementalMatrix{T}) = T

# Error is handled in error.jl as an Exception

typealias SortType Cint
const UNSORTED   = Cint(0)
const DESCENDING = Cint(1)
const ASCENDING  = Cint(2)

typealias Dist Cint
const MC		= Cint(0)
const MD		= Cint(1)
const MR		= Cint(2)
const VC		= Cint(3)
const VR		= Cint(4)
const STAR		= Cint(5)
const CIRC		= Cint(6)

typealias Orientation Cint
const NORMAL 	= Cint(0)
const TRANSPOSE = Cint(1)
const ADJOINT 	= Cint(2)

typealias UpperOrLower Cint
const LOWER = Cint(0)
const UPPER = Cint(1)

# Get MPIWorldComm
function CommWorldValue()
    r = Ref{ElComm}(0)
    ccall((:ElMPICommWorld, libEl), Cuint, (Ref{ElComm},), r)
    return r[]
end
const CommWorld = CommWorldValue()


text:
##############

@views function DamageAssembly( f::Fields2D, Th_BC::Thermal_BC, dom::ModelDomain, mat::MaterialParameters, dt::Float64 )
    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    Num       = reshape(1:ncx*ncy,ncx,ncy)
    iC        = Num
    iW        =  ones(size(Num)); iW[2:end-0,:] = Num[1:end-1,:]
    iE        =  ones(size(Num)); iE[1:end-1,:] = Num[2:end-0,:]
    iS        =  ones(size(Num)); iS[:,2:end-0] = Num[:,1:end-1]
    iN        =  ones(size(Num)); iN[:,1:end-1] = Num[:,2:end-0]
    kW        = ones(size(Num));#f.GDam .* f.lDam
    kE        = ones(size(Num));#f.GDam .* f.lDam
    kS        = ones(size(Num));#f.GDam .* f.lDam
    kN        = ones(size(Num));#f.GDam .* f.lDam
    BCW       = zeros(size(Num)); BCW[1,:]   .= Th_BC.type_W
    BCE       = zeros(size(Num)); BCE[end,:] .= Th_BC.type_E
    BCS       = zeros(size(Num)); BCS[:,1]   .= Th_BC.type_S
    BCN       = zeros(size(Num)); BCN[:,end] .= Th_BC.type_S
    # Pinned
    pC        = Th_BC.Pinned
    pW        = zeros(size(Num)); pW[2:end-0,:] = pC[1:end-1,:]
    pE        = zeros(size(Num)); pE[1:end-1,:] = pC[2:end-0,:]
    pS        = zeros(size(Num)); pS[:,2:end-0] = pC[:,1:end-1]
    pN        = zeros(size(Num)); pN[:,1:end-1] = pC[:,2:end-0]
    # Periodic stencil
    if Th_BC.type_W==0 && Th_BC.type_E==0
        println("Periodic Poisson stencil (Damage)")
        iW[  1,:] = Num[end,:]
        iE[end,:] = Num[  1,:]
    end
    # Matrix coeffcients
    fW     = 1.0 .+ 1.0.*(BCW.==1) .- 1.0.*(BCW.==2)
    fE     = 1.0 .+ 1.0.*(BCE.==1) .- 1.0.*(BCE.==2)
    fS     = 1.0 .+ 1.0.*(BCS.==1) .- 1.0.*(BCS.==2)
    fN     = 1.0 .+ 1.0.*(BCN.==1) .- 1.0.*(BCN.==2)
    cC     =  fN.*kN./dy^2 .+ fS.*kS./dy^2 .+ fE.*kE./dx^2 .+ fW.*kW./dx^2 .+ ( 2*f.We ./ (mat.gDam * mat.lDam) .+  1.0 ./mat.lDam^2 .+ mat.eDam ./ (dt.*mat.gDam * mat.lDam))
    cW     = -kW./dx^2; cW[BCW.!=0] .= 0.0
    cE     = -kE./dx^2; cE[BCE.!=0] .= 0.0
    cS     = -kS./dy^2; cS[BCS.!=0] .= 0.0
    cN     = -kN./dy^2; cN[BCN.!=0] .= 0.0
    # Pinning
    cC[pC.==1] .= 1.0; cW[pC.==1] .= 0.0; cE[pC.==1] .= 0.0; cS[pC.==1] .= 0.0; cN[pC.==1] .= 0.0
    cW[pW.==1] .= 0.0
    cE[pE.==1] .= 0.0
    cS[pS.==1] .= 0.0
    cN[pN.==1] .= 0.0
    # Sparse matrix
    I = [iC[:]; iC[:]; iC[:]; iC[:]; iC[:]]
    J = [iC[:]; iW[:]; iE[:]; iS[:]; iN[:]]
    V = [cC[:]; cW[:]; cE[:]; cS[:]; cN[:]]
    K = sparse( I, J, V)
    droptol!(K, 1e-13)
    # display(UnicodePlots.spy(K))
    return K
end
export DamageAssembly

##############

function DamageResidual!( f::Fields2D, Th_BC::Thermal_BC, dom::ModelDomain, mat::MaterialParameters, dt::Float64 )
    """ Damage residual """
    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    T = f.phiDam
    # Field with ghost values
    Tex = zeros(ncx+2, ncy+2)
    Tex[2:end-1,2:end-1] .= T
    Tex[  1,2:end-1]     .= (Th_BC.type_W==0) * T[end,:] .+ (Th_BC.type_W==2) * T[  1,:] .+ (Th_BC.type_W==1) * (2.0*Th_BC.Dir_W[:] - T[  1,:])
    Tex[end,2:end-1]     .= (Th_BC.type_E==0) * T[  1,:] .+ (Th_BC.type_E==2) * T[end,:] .+ (Th_BC.type_E==1) * (2.0*Th_BC.Dir_E[:] - T[end,:])
    Tex[2:end-1,  1]     .=                                 (Th_BC.type_S==2) * T[:,  1] .+ (Th_BC.type_S==1) * (2.0*Th_BC.Dir_S[:] - T[:,  1])
    Tex[2:end-1,end]     .=                                 (Th_BC.type_N==2) * T[:,end] .+ (Th_BC.type_N==1) * (2.0*Th_BC.Dir_N[:] - T[:,end])
    # Compute residual
    @tturbo for i = 1:ncx
        for j = 1:ncy
            pinned = Th_BC.Pinned[i,j]==1
            qW = (Tex[i+1,j+1] - Tex[i+0,j+1])/dx
            qE = (Tex[i+2,j+1] - Tex[i+1,j+1])/dx
            qS = (Tex[i+1,j+1] - Tex[i+1,j+0])/dy
            qN = (Tex[i+1,j+2] - Tex[i+1,j+1])/dy
            f.FDam[i,j] =  (pinned==0) * ( f.phiDam[i,j] * ( 2*f.We[i,j] / (mat.gDam*mat.lDam) +  1.0/mat.lDam^2 ) - 2.0*f.We[i,j] / (mat.gDam*mat.lDam) + mat.eDam*(f.phiDam[i,j]-f.phiDam0[i,j]) / (dt*mat.gDam*mat.lDam) - ( (qE-qW)/dx + (qN-qS)/dy ) ) + (pinned==1) * 0.0

            # f.FDam[i,j] =  1.0/(f.GDam[i,j]*f.lDam[i,j]) * (-2f.We[i,j] + ( f.GDam[i,j]/f.lDam[i,j] + 2f.We[i,j] )*f.phiDam[i,j]) - ( (qE-qW)/dx + (qN-qS)/dy )
        end
    end
    # println(norm(2f.We))
    # println(norm(f.FDam))
    # error("dd")
end
export DamageResidual!

##############

@views function InitialDamageAssembly( f::Fields2D, Th_BC::Thermal_BC, dom::ModelDomain, mat::MaterialParameters )
    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    Num       = reshape(1:ncx*ncy,ncx,ncy)
    iC        = Num
    iW        =  ones(size(Num)); iW[2:end-0,:] = Num[1:end-1,:]
    iE        =  ones(size(Num)); iE[1:end-1,:] = Num[2:end-0,:]
    iS        =  ones(size(Num)); iS[:,2:end-0] = Num[:,1:end-1]
    iN        =  ones(size(Num)); iN[:,1:end-1] = Num[:,2:end-0]
    kW        = ones(size(Num));
    kE        = ones(size(Num));
    kS        = ones(size(Num));
    kN        = ones(size(Num));
    BCW       = zeros(size(Num)); BCW[1,:]   .= Th_BC.type_W
    BCE       = zeros(size(Num)); BCE[end,:] .= Th_BC.type_E
    BCS       = zeros(size(Num)); BCS[:,1]   .= Th_BC.type_S
    BCN       = zeros(size(Num)); BCN[:,end] .= Th_BC.type_S
    # Pinned
    pC        = Th_BC.Pinned
    pW        = zeros(size(Num)); pW[2:end-0,:] = pC[1:end-1,:]
    pE        = zeros(size(Num)); pE[1:end-1,:] = pC[2:end-0,:]
    pS        = zeros(size(Num)); pS[:,2:end-0] = pC[:,1:end-1]
    pN        = zeros(size(Num)); pN[:,1:end-1] = pC[:,2:end-0]
    # Periodic stencil
    if Th_BC.type_W==0 && Th_BC.type_E==0
        println("Periodic Poisson stencil (Damage)")
        iW[  1,:] = Num[end,:]
        iE[end,:] = Num[  1,:]
    end
    # Matrix coeffcients
    fW     = 1.0 .+ 1.0.*(BCW.==1) .- 1.0.*(BCW.==2)
    fE     = 1.0 .+ 1.0.*(BCE.==1) .- 1.0.*(BCE.==2)
    fS     = 1.0 .+ 1.0.*(BCS.==1) .- 1.0.*(BCS.==2)
    fN     = 1.0 .+ 1.0.*(BCN.==1) .- 1.0.*(BCN.==2)
    cC     =  fN.*kN./dy^2 .+ fS.*kS./dy^2 .+ fE.*kE./dx^2 .+ fW.*kW./dx^2 .+ 1.0 / mat.lDam^2
    cW     = -kW./dx^2; cW[BCW.!=0] .= 0.0
    cE     = -kE./dx^2; cE[BCE.!=0] .= 0.0
    cS     = -kS./dy^2; cS[BCS.!=0] .= 0.0
    cN     = -kN./dy^2; cN[BCN.!=0] .= 0.0
    # Pinning
    cC[pC.==1] .= 1.0; cW[pC.==1] .= 0.0; cE[pC.==1] .= 0.0; cS[pC.==1] .= 0.0; cN[pC.==1] .= 0.0
    cW[pW.==1] .= 0.0
    cE[pE.==1] .= 0.0
    cS[pS.==1] .= 0.0
    cN[pN.==1] .= 0.0
    # Sparse matrix
    I = [iC[:]; iC[:]; iC[:]; iC[:]; iC[:]]
    J = [iC[:]; iW[:]; iE[:]; iS[:]; iN[:]]
    V = [cC[:]; cW[:]; cE[:]; cS[:]; cN[:]]
    K = sparse( I, J, V)
    droptol!(K, 1e-13)
    # display(UnicodePlots.spy(K))
    return K
end
export InitialDamageAssembly

##############

function InitialDamageResidual!( f::Fields2D, Th_BC::Thermal_BC, dom::ModelDomain, mat::MaterialParameters )
    """ Damage residual """
    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    T = f.phiDam
    # Field with ghost values
    Tex = zeros(ncx+2, ncy+2)
    Tex[2:end-1,2:end-1] .= T
    Tex[  1,2:end-1]     .= (Th_BC.type_W==0) * T[end,:] .+ (Th_BC.type_W==2) * T[  1,:] .+ (Th_BC.type_W==1) * (2.0*Th_BC.Dir_W[:] - T[  1,:])
    Tex[end,2:end-1]     .= (Th_BC.type_E==0) * T[  1,:] .+ (Th_BC.type_E==2) * T[end,:] .+ (Th_BC.type_E==1) * (2.0*Th_BC.Dir_E[:] - T[end,:])
    Tex[2:end-1,  1]     .=                                 (Th_BC.type_S==2) * T[:,  1] .+ (Th_BC.type_S==1) * (2.0*Th_BC.Dir_S[:] - T[:,  1])
    Tex[2:end-1,end]     .=                                 (Th_BC.type_N==2) * T[:,end] .+ (Th_BC.type_N==1) * (2.0*Th_BC.Dir_N[:] - T[:,end])
    # Compute residual
    @tturbo for i = 1:ncx
        for j = 1:ncy
            pinned = Th_BC.Pinned[i,j]==1
            qW = (Tex[i+1,j+1] - Tex[i+0,j+1])/dx
            qE = (Tex[i+2,j+1] - Tex[i+1,j+1])/dx
            qS = (Tex[i+1,j+1] - Tex[i+1,j+0])/dy
            qN = (Tex[i+1,j+2] - Tex[i+1,j+1])/dy
            f.FDam[i,j] = (pinned==0) * ( f.phiDam[i,j] / mat.lDam[i,j]^2 - ( (qE-qW)/dx + (qN-qS)/dy ) ) + (pinned==1) * 0.0
        end
    end
end
export InitialDamageResidual!


text:
module Contingency

# package code goes here

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the operations with quaternions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/quaternion.jl
# =========================

# Functions: +
# ------------

@testset "Operations with quaternions: +"  begin
    q1 = Quaternion(1, 0, 0, 0)
    q2 = Quaternion(0, cosd(60), 0, sind(60))
    q3 = q1 + q2
    @test q3.q0 == 1
    @test q3.q1 == cosd(60)
    @test q3.q2 == 0
    @test q3.q3 == sind(60)
    @test eltype(q3) === promote_type(eltype(q1), eltype(q2))

    q1 = Quaternion{Float32}(1, 0, 0, 0)
    q2 = Quaternion{Float32}(0, cosd(60), 0, sind(60))
    q3 = q1 + q2
    @test q3.q0 == 1
    @test q3.q1 == cosd(60.0f0)
    @test q3.q2 == 0
    @test q3.q3 == sind(60.0f0)
    @test eltype(q3) === Float32

    q3 = q1 + I
    @test q3.q0 == 2
    @test q3.q1 == 0
    @test q3.q2 == 0
    @test q3.q3 == 0
    @test eltype(q3) === promote_type(eltype(q1), eltype(I))

    q4 = I + q1
    @test q3 === q4

    q3 = q1 + 1.0I
    @test q3.q0 == 2
    @test q3.q1 == 0
    @test q3.q2 == 0
    @test q3.q3 == 0
    @test eltype(q3) === promote_type(eltype(q1), eltype(1.0))
end

# Functions: -
# ------------

@testset "Operations with quaternions: -"  begin
    q1 = Quaternion(1, 2, 3, 4)
    q2 = -q1
    @test q2.q0 == -1
    @test q2.q1 == -2
    @test q2.q2 == -3
    @test q2.q3 == -4
    @test eltype(q2) === eltype(q1)

    q1 = Quaternion(1, 0, 0, 0)
    q2 = Quaternion(0, cosd(60), 0, sind(60))
    q3 = q1 - q2
    @test q3.q0 == 1
    @test q3.q1 == -cosd(60)
    @test q3.q2 == 0
    @test q3.q3 == -sind(60)
    @test eltype(q3) === promote_type(eltype(q1), eltype(q2))

    q1 = Quaternion{Float32}(1, 0, 0, 0)
    q2 = Quaternion{Float32}(0, cosd(60), 0, sind(60))
    q3 = q1 - q2
    @test q3.q0 == 1
    @test q3.q1 == -cosd(60.0f0)
    @test q3.q2 == 0
    @test q3.q3 == -sind(60.0f0)
    @test eltype(q3) === Float32

    q3 = q1 - I
    @test q3.q0 == 0
    @test q3.q1 == 0
    @test q3.q2 == 0
    @test q3.q3 == 0
    @test eltype(q3) === promote_type(eltype(q1), eltype(I))

    q4 = I - q1
    @test q4 == -q3

    q3 = q1 - 1.0I
    @test q3.q0 == 0
    @test q3.q1 == 0
    @test q3.q2 == 0
    @test q3.q3 == 0
    @test eltype(q3) === promote_type(eltype(q1), eltype(1.0))
end

# Functions: *
# ------------

@testset "Operations with quaternions: *" begin
    # Scalar multiplication.
    q1 = Quaternion(1, 0, 0, 0)

    q2 = 2q1
    @test q2.q0 == 2
    @test q2.q1 == 0
    @test q2.q2 == 0
    @test q2.q3 == 0
    @test eltype(q2) === eltype(q1)

    q3 = q1 * 2
    @test q3 === q2

    q2 = 2.0q1
    @test q2.q0 == 2
    @test q2.q1 == 0
    @test q2.q2 == 0
    @test q2.q3 == 0
    @test eltype(q2) === eltype(2.0)

    q3 = q1 * 2.0
    @test q3 === q2

    # Quaternion multiplication (Hamilton product).
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q3 = q1 * q2

    @test q3.q0 == q1.q0 * q2.q0 - q1.q1 * q2.q1 - q1.q2 * q2.q2 - q1.q3 * q2.q3
    @test q3.q1 == q1.q0 * q2.q1 + q1.q1 * q2.q0 + q1.q2 * q2.q3 - q1.q3 * q2.q2
    @test q3.q2 == q1.q0 * q2.q2 - q1.q1 * q2.q3 + q1.q2 * q2.q0 + q1.q3 * q2.q1
    @test q3.q3 == q1.q0 * q2.q3 + q1.q1 * q2.q2 - q1.q2 * q2.q1 + q1.q3 * q2.q0
    @test eltype(q3) === Float64

    q3 = I * q1
    @test q3 === q1

    q3 = q1 * I
    @test q3 === q1

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Int}(1, 1, 1, 1)
    q3 = q1 * q2

    @test q3.q0 == q1.q0 * q2.q0 - q1.q1 * q2.q1 - q1.q2 * q2.q2 - q1.q3 * q2.q3
    @test q3.q1 == q1.q0 * q2.q1 + q1.q1 * q2.q0 + q1.q2 * q2.q3 - q1.q3 * q2.q2
    @test q3.q2 == q1.q0 * q2.q2 - q1.q1 * q2.q3 + q1.q2 * q2.q0 + q1.q3 * q2.q1
    @test q3.q3 == q1.q0 * q2.q3 + q1.q1 * q2.q2 - q1.q2 * q2.q1 + q1.q3 * q2.q0
    @test eltype(q3) === promote_type(Int, Float32)

    # Quaternion multiplication with a vector.
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    v1 = randn(3)

    q2 = q1 * v1
    q3 = q1 * Quaternion(0, v1)
    @test q2 === q3

    q2 = v1 * q1
    q3 = Quaternion(0, v1) * q1
    @test q2 === q3

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    v1 = [1, 2, 3]

    q2 = q1 * v1
    q3 = q1 * Quaternion(0, v1)
    @test q2 === q3

    q2 = v1 * q1
    q3 = Quaternion(0, v1) * q1
    @test q2 === q3
end

# Functions: /
# ------------

@testset "Operations with quaternions: /" begin
    # Scalar division.
    q1 = Quaternion(2, 0, 0, 0)
    q2 = q1 / 2
    @test q2.q0 == 1
    @test q2.q1 == 0
    @test q2.q2 == 0
    @test q2.q3 == 0

    q1 = Quaternion(2, 1, 0, 1)
    q2 = 2 / q1
    @test q2.q0 == +2 / 3
    @test q2.q1 == -1 / 3
    @test q2.q2 == 0
    @test q2.q3 == -1 / 3

    q1 = Quaternion(1, 0, 1, 0)
    q2 = 2.0q1
    @test q2.q0 == 2
    @test q2.q1 == 0
    @test q2.q2 == 2
    @test q2.q3 == 0
    @test eltype(q2) === eltype(2.0)

    # Quaternion division (Hamilton product).
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q3 = q1 / q2
    norm_q2² = norm(q2) * norm(q2)

    @test q3.q0 ≈ (+ q1.q0 * q2.q0 + q1.q1 * q2.q1 + q1.q2 * q2.q2 + q1.q3 * q2.q3)/norm_q2²
    @test q3.q1 ≈ (- q1.q0 * q2.q1 + q1.q1 * q2.q0 - q1.q2 * q2.q3 + q1.q3 * q2.q2)/norm_q2²
    @test q3.q2 ≈ (- q1.q0 * q2.q2 + q1.q1 * q2.q3 + q1.q2 * q2.q0 - q1.q3 * q2.q1)/norm_q2²
    @test q3.q3 ≈ (- q1.q0 * q2.q3 - q1.q1 * q2.q2 + q1.q2 * q2.q1 + q1.q3 * q2.q0)/norm_q2²
    @test eltype(q3) === Float64

    q3 = I / q1
    @test q3 ≈ inv(q1)

    q3 = q1 / I
    @test q3 === q1

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Float32}(1, 1, 1, 1)
    q3 = q1 / q2
    norm_q2² = norm(q2) * norm(q2)

    @test q3.q0 ≈ (+ q1.q0 * q2.q0 + q1.q1 * q2.q1 + q1.q2 * q2.q2 + q1.q3 * q2.q3)/norm_q2²
    @test q3.q1 ≈ (- q1.q0 * q2.q1 + q1.q1 * q2.q0 - q1.q2 * q2.q3 + q1.q3 * q2.q2)/norm_q2²
    @test q3.q2 ≈ (- q1.q0 * q2.q2 + q1.q1 * q2.q3 + q1.q2 * q2.q0 - q1.q3 * q2.q1)/norm_q2²
    @test q3.q3 ≈ (- q1.q0 * q2.q3 - q1.q1 * q2.q2 + q1.q2 * q2.q1 + q1.q3 * q2.q0)/norm_q2²
    @test eltype(q3) === promote_type(Int, Float32)
end

# Functions: \
# ------------

@testset "Operations with quaternions: \\" begin
    # Quaternion division (Hamilton product).
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q3 = q1 \ q2
    q4 = inv(q1) * q2
    @test q3 ≈ q4

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q2 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q3 = q1 \ q2
    q4 = inv(q1) * q2
    @test q3 ≈ q4
    @test eltype(q3) == Float32

    # Quaternion multiplication with a vector.
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    v1 = randn(3)
    q2 = v1 \ q1
    q3 = inv(Quaternion(0, v1)) * q1
    @test q2 ≈ q3
    @test eltype(q2) == Float64

    q2 = q1 \ v1
    q3 = inv(q1) * Quaternion(0, v1)
    @test q2 ≈ q3
    @test eltype(q2) == Float64

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    v1 = randn(Float32, 3)
    q2 = v1 \ q1
    q3 = inv(Quaternion(0, v1)) * q1
    @test q2 ≈ q3
    @test eltype(q2) == Float32

    q2 = q1 \ v1
    q3 = inv(q1) * Quaternion(0, v1)
    @test q2 ≈ q3
    @test eltype(q2) == Float32

    # Operations with `I`.
    q1 = Quaternion{Float64}(randn(), randn(), randn(), randn())
    q2 = I \ q1
    @test q2 == q1

    q2 = q1 \ I
    @test q2 == inv(q1)

    q1 = Quaternion{Float32}(randn(), randn(), randn(), randn())
    q2 = I \ q1
    @test q2 == q1

    q2 = q1 \ I
    @test q2 == inv(q1)
end


text:
module FieldTracer

# Hongyang Zhou, hyzhou@umich.edu

using Meshes, Requires, MuladdMacro

export trace
export Euler, RK4

abstract type Algorithm end

struct Euler <: Algorithm end
struct RK4 <: Algorithm end

include("structured2d.jl")
include("structured3d.jl")
include("unstructured2d.jl")
include("utility/seed.jl")

"""
	 trace(fieldx, fieldy, startx, starty, gridx, gridy; alg=RK4(), kwargs...)

Stream tracing on structured mesh with field in 2D array and grid in range.
"""

"""
	 trace(fieldx, fieldy, fieldz, startx, starty, startz, gridx, gridy, gridz;
       alg=RK4(), kwargs...)

    trace(fieldx, fieldy, fieldz, startx, starty, startz, grid::CartesianGrid;
		 alg=RK4(), maxstep=20000, ds=0.01, gridType="ndgrid", direction="both")

Stream tracing on structured mesh with field in 3D array and grid in range.
"""
function trace(args...; alg::Algorithm=RK4(), kwargs...)

   if length(args) ≤ 6 # 2D
      if alg isa RK4
         trace2d_rk4(args...; kwargs...)
      elseif alg isa Euler
         trace2d_euler(args...; kwargs...)
      end
   else # 3D
      if alg isa RK4
         trace3d_rk4(args...; kwargs...)
      elseif alg isa Euler
         trace3d_euler(args...; kwargs...)
      end
   end

end

function __init__()
   @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
      include("utility/pyplot.jl")
  end
end

end


text:
"""
    InflationTrimmedMeanWeighted <: InflationFunction

Función de inflación para computar la media truncada ponderada

## Utilización
    function (inflfn::InflationTrimmedMeanWeighted)(base::VarCPIBase{T}) where T

Define cómo opera InflationTrimmedMeanWeighted sobre un objeto de tipo VarCPIBase.
"""
Base.@kwdef struct InflationTrimmedMeanWeighted <: InflationFunction
    l1::Float32
    l2::Float32
end

# Métodos para crear funciones de inflación a partir de enteros
"""
    InflationTrimmedMeanWeighted(l1::Real,l2::Real)
Nos permite utilizar valores que no necesariamente son Float32, como enteros o Float64.

# Ejemplo:
```julia-repl
julia> mtfn = InflationTrimmedMeanWeighted(25,75.5)
(::InflationTrimmedMeanWeighted) (generic function with 5 methods)
```
"""
InflationTrimmedMeanWeighted(l1::Real,l2::Real) = InflationTrimmedMeanWeighted(l1 = Float32(l1), l2 = Float32(l2))

"""
    measure_name(inflfn::InflationTrimmedMeanWeighted)

Nos indica qué medida se utiliza para una instancia de una función de inflación.

# Ejemplo:
```julia-repl
julia> mtfn = InflationTrimmedMeanWeighted(15.5,75.5)
julia> measure_name(mtfn)
"Media Truncada Ponderada (15.5 , 75.5)"
```
"""
function measure_name(inflfn::InflationTrimmedMeanWeighted)
    l1 = string(round(inflfn.l1, digits=2))
    l2 = string(round(inflfn.l2, digits=2))
    "Media Truncada Ponderada (" * l1 * ", " * l2 * ")"
end

# Extendemos `params`, que devuelve los parámetros de la medida de inflación
CPIDataBase.params(inflfn::InflationTrimmedMeanWeighted) = (inflfn.l1, inflfn.l2)

# Operación de InflationTrimmedMeanWeighted sobre VarCPIBase para obtener el
# resumen intermensual de esta metodología
function (inflfn::InflationTrimmedMeanWeighted)(base::VarCPIBase{T}) where T
    l1 = inflfn.l1
    l2 = inflfn.l2
    # l1 = min(inflfn.l1, inflfn.l2)
    # l2 = max(inflfn.l1, inflfn.l2)
    outVec   = Vector{T}(undef, periods(base))
    # para cada t: creamos parejas de variaciones con pesos,
    # ordenamos de acuerdo a variaciones, truncamos
    # renormalizamos para que los pesos sumen 1
    # sumamos el producto de variaciones con pesos

    # Número de gastos básicos
    g = size(base.v, 2)

    # Reservar la memoria para cómputos de media truncada
    sort_ids = [zeros(Int, g) for _ in 1:Threads.nthreads()]
    w_sorted_list = [zeros(T, g) for _ in 1:Threads.nthreads()]
    w_sorted_renorm_list = [zeros(T, g) for _ in 1:Threads.nthreads()]

    Threads.@threads for i in 1:periods(base)

        # Obtener los vectores de cada hilo
        j = Threads.threadid()
        sort_idx = sort_ids[j]
        w_sorted_acum = w_sorted_list[j]
        w_sorted_renorm = w_sorted_renorm_list[j]

        # Obtener índices de orden en sort_idx
        v_month = @view base.v[i, :]
        sortperm!(sort_idx, v_month)

        # Acumular las ponderaciones en w_sorted_acum
        w_sorted = @view base.w[sort_idx]
        cumsum!(w_sorted_acum, w_sorted)

        # Poner a cero las ponderaciones fuera de límites
        @inbounds for x in 1:g
            if w_sorted_acum[x] < l1 || w_sorted_acum[x] > l2
                w_sorted_acum[x] = 0
            else
                w_sorted_acum[x] = 1
            end
        end

        # Renormalizar las ponderaciones dentro de los límites
        w_sorted_renorm .= (w_sorted .* w_sorted_acum)
        w_sorted_renorm ./= sum(w_sorted_renorm)

        # Computar promedio ponderado de variaciones dentro de límites
        @inbounds outVec[i] = sum((@view v_month[sort_idx]) .* w_sorted_renorm)
    end

    outVec
end
# Método para recibir argumentos en forma de tupla
InflationTrimmedMeanWeighted(factors::Tuple{Real, Real}) = InflationTrimmedMeanWeighted(
    convert(Float32, factors[1]),
    convert(Float32, factors[2])
)

# Método para recibir argumentos en forma de vector
function InflationTrimmedMeanWeighted(factor_vec::Vector{<:Real})
    length(factor_vec) != 2 && return @error "Dimensión incorrecta del vector"
    InflationTrimmedMeanWeighted(
        convert(Float32, factor_vec[1]),
        convert(Float32, factor_vec[2])
    )
end


text:
module MiniTB

using TBComponents
using Parameters
#
# A toy tax-benefit system with outlines of the components
# a real model would need: models of people (and households)
# a parameter system, a holder for results, and some calculations
# using those things.
# Used in test building budget constraints.
# There's also some experiments of mine with constructors
# and copying strucs.
#

export calculate, DEFAULT_PERSON, modifiedcopy, TBParameters, Person, getnet
export modifiedcopy, DEFAULT_PARAMS, ZERO_PARAMS
export Gender, Male, Female, DEFAULT_WAGE, DEFAULT_HOURS
export NetType, NetIncome, TotalTaxes, BenefitsOnly
export calculatetax, calculatebenefit1, calculatebenefit2, calculate_internal

@enum NetType NetIncome TotalTaxes BenefitsOnly
@enum Gender Male Female

# experiment with types
const NullableFloat = Union{Missing,Float64}
const NullableInt = Union{Missing,Integer}
const NullableArray = Union{Missing,Array{Float64}}

mutable struct Person
   pid::BigInt
   wage::Float64
   hours::Float64
   age::Integer
   sex::Gender
end

mutable struct Household
   hid :: Integer
   rent::Float64
   people::Vector{Person}
end

const DEFAULT_HOURS = 30
const DEFAULT_WAGE = 5.0

const DEFAULT_PERSON = Person( BigInt(1), DEFAULT_HOURS*DEFAULT_WAGE, DEFAULT_HOURS, 40, Female)
const DEFAULT_HOUSEHOLD = Household( 1, 200.0, [DEFAULT_PERSON])

function modifiedcopy(
   copyFrom::Person;
   wage::NullableFloat = missing,
   hours::NullableFloat,
   age::NullableInt = missing,
)::Person

   Person(
      copyFrom.pid,
      wage !== missing ? wage : copyFrom.wage,
      hours !== missing ? hours : copyFrom.hours,
      age !== missing ? age : copyFrom.age,
      copyFrom.sex,
   )
end

mutable struct TBParameters
   it_allow::Float64
   it_rate::RateBands
   it_band::RateBands

   benefit1::Float64
   benefit2::Float64
   ben2_min_hours::Float64
   ben2_taper::Float64
   ben2_u_limit::Float64
   basic_income::Float64

   # attempt a constructor with named TBParameters
   function TBParameters(
      ;
      it_allow::Float64,
      it_rate::Array{Float64},
      it_band::Array{Float64},

      benefit1::Float64,
      benefit2::Float64,
      ben2_min_hours::Float64,
      ben2_taper::Float64,
      ben2_u_limit::Float64,
      basic_income::Float64
   )
      new(
         it_allow,
         it_rate,
         it_band,
         benefit1,
         benefit2,
         ben2_min_hours,
         ben2_taper,
         ben2_u_limit,
         basic_income
      )
   end
end

#
# Just a test of an idea
# e.g newpars = modifiedcopy( DEFAULT_PARAMS, it_allow=3_000 )
#
function modifiedcopy(
   copyFrom::TBParameters;
   it_allow::NullableFloat = missing,
   it_rate::NullableArray = missing,
   it_band::NullableArray = missing,

   benefit1::NullableFloat = missing,
   benefit2::NullableFloat = missing,
   ben2_min_hours::NullableFloat = missing,
   ben2_taper::NullableFloat = missing,
   ben2_u_limit::NullableFloat = missing,
   basic_income::NullableFloat = missing,

)::TBParameters

   x = it_allow !== missing ? it_allow : copyFrom.it_allow
   TBParameters(
      it_allow = it_allow !== missing ? it_allow : copyFrom.it_allow,
      it_rate = it_rate !== missing ? it_rate : copyFrom.it_rate,
      it_band = it_band !== missing ? it_band : copyFrom.it_band,

      benefit1 = benefit1 !== missing ? benefit1 : copyFrom.benefit1,
      benefit2 = benefit2 !== missing ? benefit2 : copyFrom.benefit2,
      ben2_min_hours = ben2_min_hours !== missing ? ben2_min_hours : copyFrom.ben2_min_hours,
      ben2_taper = ben2_taper !== missing ? ben2_taper : copyFrom.ben2_taper,
      ben2_u_limit = ben2_u_limit !== missing ? ben2_u_limit : copyFrom.ben2_u_limit,
      basic_income = basic_income !== missing ? basic_income : copyFrom.basic_income
   )
end

const DEFAULT_PARAMS = TBParameters(
   it_allow = weeklyise(12_500),
   it_rate = [0.20, 0.4],
   it_band = [weeklyise(50_000), 9999999999999999999.99],
   benefit1 = 73.00,
   benefit2 = 101.0, # weeklyise( 1_960.0+ 545+2_780.0),
   ben2_min_hours = 30.0,
   ben2_taper = 0.41,
   ben2_u_limit = 123.00,
   basic_income = 0.0
)

const ZERO_PARAMS = TBParameters(
   it_allow = 0.0,
   it_rate = [0.0,0.0],
   it_band = [weeklyise(50_000),99999999999999999999.99],
   benefit1 = 0.0,
   benefit2 = 0.0,
   ben2_min_hours = 0.0,
   ben2_taper = 0.0,
   ben2_u_limit = 0.0,
   basic_income = 0.0
)

const Results = Dict{Symbol,Any}

## need to include taxcalcs higher up

function calculatetax(pers::Person, params::TBParameters)::Float64
   taxable = max(0.0, pers.wage - params.it_allow)
   tc::TaxResult = calctaxdue(
      taxable    = taxable,
      rates      = params.it_rate,
      thresholds = params.it_band,
   )
   return tc.due
end

function calculatebenefit1(pers::Person, params::TBParameters)::Float64
   ben = params.benefit1 - pers.wage
   return max(0.0, ben )
end

function calculatebenefit2(pers::Person, params::TBParameters)::Float64
   b = pers.hours >= params.ben2_min_hours ? params.benefit2 : 0.0
   if pers.wage > params.ben2_u_limit
      b = max(0.0, b - (params.ben2_taper * (pers.wage - params.ben2_u_limit)))
   end
   return b
end

INCR = 0.001

function calculate_internal(pers::Person, params::TBParameters)::Results
   res = Results()
   res[:tax] = calculatetax(pers, params)
   res[:benefit1] = calculatebenefit1(pers, params)
   res[:benefit2] = calculatebenefit2(pers, params)
   res[:basic_income] = params.basic_income

   res[:netincome] = pers.wage + res[:benefit1] + res[:benefit2]+res[:basic_income] - res[:tax]
   return res
end

function calculate( pers::Person, params::TBParameters)::Results
   res1 = calculate_internal( pers, params )
   pers.wage += INCR
   res2 = calculate_internal( pers, params )
   m_metr = (res2[:netincome]-res1[:netincome])/INCR
   res1[:metr] = round( 1.0-m_metr, digits=4)
   tax_credit = res1[:netincome] - (m_metr*pers.wage)
   res1[:tax_credit] = round( tax_credit, digits=4)
   pers.wage -= INCR
   res1
end

end


text:
__precompile__()

module GggenomeSearch

export
	getResponse,
	generateGggenomeUrl,
	extractTopResultFromResponse,
	showGggenomeTopHit

using Requests

function getResponse(url::String)
    res = get(url)
	return(String(res.data))
end

function generateGggenomeUrl(seq::AbstractString, db::String, format::String)
	return(@sprintf "https://gggenome.dbcls.jp/%s/%s%s" db seq format)
end

function extractTopResultFromResponse(res_str::String)
	topResult = "No hit"
	for ln in split(chomp(res_str), "\n")
		if ln[1] == '#'
			continue
		else
			topResult = split(chomp(ln), "\t")[1]
			break
		end
	end
	return(topResult)
end

function showGggenomeTopHit(seq::AbstractString, db::String)
	url = generateGggenomeUrl(seq, db, ".txt")
	res_str = getResponse(url)
	return(extractTopResultFromResponse(res_str))
end

end # module


text:
## Primitive Object

abstract type PrimitiveObject <: Object end

export Empty
struct Empty <: PrimitiveObject end

export Sphere
struct Sphere <: PrimitiveObject
    center :: Array{Float64,1}
    radius :: Float64

    function Sphere(center,radius)
        if length(center) ≠ 3
            error("center of sphere must be an element of ℝ³")
        elseif radius == 0
            Empty()
        else
            new(center,abs(radius))
        end
    end
    function Sphere(point₁::RealVector,point₂::RealVector,point₃::RealVector,point₄::RealVector)
        center=Circumcenter(point₁,point₂,point₃,point₄)
        radius=Circumradius(point₁,point₂,point₃,point₄)
        return Sphere(center,radius)
    end
end

export Cylinder
struct Cylinder <: PrimitiveObject
    end1 :: Array{Float64,1}
    end2 :: Array{Float64,1}
    radius :: Float64

    function Cylinder(end1,end2,radius)
        if length(end1) ≠ 3
            error("end of cylinder must be an element of ℝ³")
        elseif length(end2) ≠ 3
            error("end of cylinder must be an element of ℝ³")
        elseif radius == 0
            Empty()
        elseif norm(end2-end1) == 0
            Empty()
        else
            new(end1,end2,abs(radius))
        end
    end
end

export Cone
struct Cone <: PrimitiveObject
    end1 :: Array{Float64,1}
    end2 :: Array{Float64,1}
    radius :: Float64

    function Cone(end1,end2,radius)
        if length(end1) ≠ 3
            error("end of cone must be an element of ℝ³")
        elseif length(end2) ≠ 3
            error("end of cone must be an element of ℝ³")
        elseif radius == 0
            Empty()
        elseif norm(end2-end1) == 0
            Empty()
        else
            new(end1,end2,abs(radius))
        end
    end
end

export Box
struct Box <: PrimitiveObject
    vertex1 :: Array{Float64,1}
    vertex2 :: Array{Float64,1}

    function Box(vertex1,vertex2)
        if length(vertex1) ≠ 3
            error("vertex of box must be an element of ℝ³")
        elseif length(vertex2) ≠ 3
            error("vertex of box must be an element of ℝ³")
        elseif norm(vertex2-vertex1) == 0
            Empty()
        else
            new(vertex1,vertex2)
        end
    end
end

export Disc
struct Disc <: PrimitiveObject
    center :: Array{Float64,1}
    normal :: Array{Float64,1}
    radius :: Float64

    function Disc(center,normal,radius)
        if length(center) ≠ 3
            error("center of disc must be an element of ℝ³")
        elseif length(normal) ≠ 3
            error("normal of disc must be an element of ℝ³")
        elseif norm(normal) == 0
            error("normal vector must be non-zero")
        elseif radius == 0
            Empty()
        else
            new(center,normal,abs(radius))
        end
    end
    function Disc(p₁::RealVector, p₂::RealVector, p₃::RealVector)
        center=Circumcenter(p₁,p₂,p₃)
        normal=NormalVector(p₁,p₂,p₃)
        R=norm(p₁-center)
        return Disc(center,normal,R)
    end
end

export Torus
struct Torus <: PrimitiveObject
    radius1 :: Float64
    radius2 :: Float64

    function Torus(radius1,radius2)
        if radius2 == 0
            Empty()
        elseif radius1 == 0
            Sphere([0,0,0],abs(radius2))
        else
            new(abs(radius1),abs(radius2))
        end
    end
    function Torus(center::RealVector, normal::RealVector, R::Float64, r::Float64)
        e₁=normalize(normal)
        e₂=OrthogonalVector(e₁)
        e₃=cross(e₁,e₂)
        A=hcat(e₃,e₁,e₂)
        return AffineTransform(Torus(R,r),A,center)
    end
    function Torus(p₁::RealVector, p₂::RealVector, p₃::RealVector, r::Float64; ε=10^(-4))
        α₁,α₂,α₃=Angles(p₁,p₂,p₃)
        if α₁>π-ε
            return Line(p₂,p₃,r)
        elseif α₂>π-ε
            return Line(p₃,p₁,r)
        elseif α₃>π-ε
            return Line(p₁,p₂,r)
        else
            center=Circumcenter(p₁,p₂,p₃)
            normal=NormalVector(p₁,p₂,p₃)
            R=(norm(p₁-center)+norm(p₂-center)+norm(p₃-center))/3
            return Torus(center,normal,R,r)
        end
    end

end

export Polygon
struct Polygon <: PrimitiveObject
    vertices :: Array{Array{Float64,1},1}

    function Polygon(vertices)
        if !all(e->e==3,length.(vertices))
            error("vertices of polygon must be an element of ℝ³")
        elseif rank(hcat(vertices...)-repeat(+(vertices...)/length(vertices),1,length(vertices)),atol=10^(-10)) ≠ 2
            Empty()
        else
            new(vertices)
        end
    end
    function Polygon(vertices::RealVector...)
        return Polygon([vertices...])
    end
end

function povray_script(sphere :: Sphere; preindent = 0)
    return "sphere{"*povray_script(sphere.center)*","*povray_script(sphere.radius)*"}"
end
function povray_script(cylinder :: Cylinder; preindent = 0)
    return "cylinder{"*povray_script(cylinder.end1)*","*povray_script(cylinder.end2)*","*povray_script(cylinder.radius)*"}"
end
function povray_script(cone :: Cone; preindent = 0)
    return "cone{"*povray_script(cone.end1)*","*povray_script(cone.radius)*","*povray_script(cone.end2)*",0}"
end
function povray_script(box :: Box; preindent = 0)
    return "box{"*povray_script(box.vertex1)*","*povray_script(box.vertex2)*"}"
end
function povray_script(disc :: Disc; preindent = 0)
    return "disc{"*povray_script(disc.center)*","*povray_script(disc.normal)*","*povray_script(disc.radius)*"}"
end
function povray_script(torus :: Torus; preindent = 0)
    return "torus{"*povray_script(torus.radius1)*","*povray_script(torus.radius2)*"}"
end
function povray_script(polygon :: Polygon; preindent = 0)
    n = length(polygon.vertices)
    return "polygon{" * povray_script(n) * ", " * povray_script(polygon.vertices) * "}"
end

import Base: isempty

function isempty(object::Object)
    return repr(object) == "Empty()"
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
function ChainRulesCore.rrule(
    ::typeof(linsolve),
    A::AbstractMatrix,
    b::AbstractVector,
    x₀,
    algorithm,
    a₀,
    a₁
)
    (x, info) = KrylovKit.linsolve(A, b, x₀, algorithm, a₀, a₁)
    project_A = ProjectTo(A)

    function linsolve_pullback(X̄)
        x̄ = unthunk(X̄[1])
        ∂self = NoTangent()
        ∂x₀ = ZeroTangent()
        ∂algorithm = NoTangent()
        ∂b, reverse_info = KrylovKit.linsolve(
            A', x̄, (zero(a₀) * zero(a₁)) * x̄, algorithm, conj(a₀), conj(a₁)
        )
        if info.converged > 0 && reverse_info.converged == 0
            @warn "The cotangent linear problem did not converge, whereas the primal linear problem did."
        end
        if A isa StridedMatrix
            ∂A = InplaceableThunk(
                Ā -> mul!(Ā, ∂b, x', -conj(a₁), true),
                @thunk(-conj(a₁) * ∂b * x'),
            )
        else
            ∂A = @thunk(project_A(-conj(a₁) * ∂b * x'))
        end
        ∂a₀ = @thunk(-dot(x, ∂b))
        if a₀ == zero(a₀) && a₁ == one(a₁)
            ∂a₁ = @thunk(-dot(b, ∂b))
        else
            ∂a₁ = @thunk(-dot((b - a₀ * x)/a₁, ∂b))
        end
        return ∂self, ∂A, ∂b, ∂x₀, ∂algorithm, ∂a₀, ∂a₁
    end
    return (x, info), linsolve_pullback
end

function ChainRulesCore.rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(linsolve),
    A::AbstractMatrix,
    b::AbstractVector,
    x₀,
    algorithm,
    a₀,
    a₁
)
    return rrule(linsolve, A, b, x₀, algorithm, a₀, a₁)
end

function ChainRulesCore.rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(linsolve),
    f,
    b,
    x₀,
    algorithm,
    a₀,
    a₁
)
    x, info = KrylovKit.linsolve(f, b, x₀, algorithm, a₀, a₁)

    # f defines a linear map => pullback defines action of the adjoint
    (y, f_pullback) = rrule_via_ad(config, f, x)
    fᴴ(xᴴ) = f_pullback(xᴴ)[2]
    # TODO can we avoid computing f_pullback if algorithm isa Union{CG,MINRES}?

    function linsolve_pullback(X̄)
        x̄ = unthunk(X̄[1])
        ∂self = NoTangent()
        ∂x₀ = ZeroTangent()
        ∂algorithm = NoTangent()
        (∂b, reverse_info) = KrylovKit.linsolve(
            fᴴ, x̄, (zero(a₀) * zero(a₁)) * x̄, algorithm, conj(a₀), conj(a₁)
        )
        if reverse_info.converged == 0
            @warn "Linear problem for reverse rule did not converge." reverse_info
        end
        ∂f = @thunk(f_pullback(-conj(a₁)*∂b)[1])
        ∂a₀ = @thunk(-dot(x, ∂b))
        # ∂a₁ = @thunk(-dot(f(x), ∂b))
        if a₀ == zero(a₀) && a₁ == one(a₁)
            ∂a₁ = @thunk(-dot(b, ∂b))
        else
            ∂a₁ = @thunk(-dot((b - a₀ * x)/a₁, ∂b))
        end
        return ∂self, ∂f, ∂b, ∂x₀, ∂algorithm, ∂a₀, ∂a₁
    end
    return (x, info), linsolve_pullback
end

# frule - currently untested

function ChainRulesCore.frule((_, ΔA, Δb, Δx₀, _, Δa₀, Δa₁)::Tuple, ::typeof(linsolve),
    A::AbstractMatrix, b::AbstractVector, x₀, algorithm, a₀, a₁)

    (x, info) = KrylovKit.linsolve(A, b, x₀, algorithm, a₀, a₁)

    if Δb isa ChainRulesCore.AbstractZero
        rhs = zero(b)
    else
        rhs = (1 - Δa₁)*Δb
    end
    if !iszero(Δa₀)
        rhs = axpy!(-Δa₀, x, rhs)
    end
    if !iszero(ΔA)
        rhs = mul!(rhs, ΔA, x, -a₁, true)
    end
    (Δx, forward_info) = KrylovKit.linsolve(A, rhs, zero(rhs), algorithm, a₀, a₁)
    if info.converged > 0 && forward_info.converged == 0
        @warn "The tangent linear problem did not converge, whereas the primal linear problem did."
    end
    return (x, info), (Δx, NoTangent())
end

function ChainRulesCore.frule(config::RuleConfig{>:HasForwardsMode}, tangents, ::typeof(linsolve),
    A::AbstractMatrix, b::AbstractVector, x₀, algorithm, a₀, a₁)
    frule(tangents, linsolve, A, b, x₀, algorithm, a₀, a₁)
end

function ChainRulesCore.frule(config::RuleConfig{>:HasForwardsMode}, (_, Δf, Δb, Δx₀, _, Δa₀, Δa₁),
    ::typeof(linsolve),
    f, b, x₀, algorithm, a₀, a₁)

    (x, info) = KrylovKit.linsolve(f, b, x₀, algorithm, a₀, a₁)

    if Δb isa AbstractZero
        rhs = false * b
    else
        rhs = (1 - Δa₁) * Δb
    end
    if !iszero(Δa₀)
        rhs = axpy!(-Δa₀, x, rhs)
    end
    if !(Δf isa AbstractZero)
        rhs = axpy!(-a₁, frule_via_ad(config, (Δf, ZeroTangent()), f, x), rhs)
    end
    (Δx, forward_info) = KrylovKit.linsolve(f, rhs, false*rhs, algorithm, a₀, a₁)
    if info.converged > 0 && forward_info.converged == 0
        @warn "The tangent linear problem did not converge, whereas the primal linear problem did."
    end
    return (x, info), (Δx, NoTangent())
end


text:
using Glob.glob

const FORTRAN_FUNCTIONS =
    [:dnrm2, :dasum, :ddot, :idamax, :dgemm, :dgemv, :dger,
     :dtrsm, :dtrmv, :dpotrf, :dpotrs, :dpotri, :dtrtri]

"""
    find_obj([Makefile])

Parse the original Csdp `Makefile` and return all the object names
(without extension `.o`) for `libsdp`.
Remark: You cannot use the Makefile directly as it produces only a
static library.
"""
function find_obj(makefile_path=Makefile)
     makefile = readstring(makefile_path)
    m = match(r"libsdp\.a\:(.+)", makefile)
    m != nothing || error("Could not find `libsdp.a` target in '$makefile_path'")
    objs = matchall(r"\w+\.o", m.captures[1])
    objs = String[splitext(o)[1] for o in [objs; basename(patchf)]]
end

"""
    patch_int()

Modify the files `lib/*.c` and the headers in `include/` to use `integer`
instead of `int`; the `printf` format `%d` is adapted, accordingly.
Notify that this happens *inplace*, i.e. without backup.
"""
function patch_int(; verbose::Bool = false)
    let patchsrc = "$srcdir/$(basename(patchf))"
        isfile(patchsrc) || cp(patchf, patchsrc)
    end
    if JULIA_LAPACK
        info("Patching int --> integer")
        cfiles = [glob("*.c", srcdir); [joinpath(srcdir, "..", "include", "$d.h")
                                        for d in ["declarations",
                                                  "blockmat",
                                                  "parameters"]]]
        for cfile in cfiles
            if verbose; println(cfile); end
            content = readstring(cfile)
            for (re,subst) in
                [(r"int ([^(]+);", s"integer \1;"),
                 (r"int ", s"integer "),
                 (r"int \*", s"integer *"),
                 (r"integer mycompare", s"int mycompare"),
                 (r"\(int\)", s"(integer)"),
                 (r"%d", s"%ld"),
                 (r"%2d", s"%2ld"),
                 (r" int ", s" integer ")]
                content = replace(content, re, subst)
            end
            open(cfile, "w") do io
                print(io, content)
            end
        end
    end
end

function compile_objs(JULIA_LAPACK=JULIA_LAPACK)
    if JULIA_LAPACK
        lapack = Libdl.dlpath(LinAlg.LAPACK.liblapack)
        lflag = replace(splitext(basename(lapack))[1], r"^lib", "")
        libs = ["-L$(dirname(lapack))", "-l$lflag"]
        info(libs)
        if endswith(LinAlg.LAPACK.liblapack, "64_")
            push!(cflags, "-march=x86-64", "-m64", "-Dinteger=long")
            for f in FORTRAN_FUNCTIONS
                let ext=string(BLAS.@blasfunc "")
                    push!(cflags, "-D$(f)_=$(f)_$ext")
                end
            end
            info(cflags)
        end
    else
        libs = ["-l$l" for l in ["blas", "lapack"]]
        @static if is_windows() unshift!(libs, "-L$libdir") end
    end

    for o in find_obj()
        info("CC $o.c")
        @static if is_unix()  push!(cflags, "-fPIC") end
        run(`$CC $cflags -o $builddir/$o.o -c $srcdir/$o.c`)
    end
    objs = ["$builddir/$o.o" for o in find_obj()]
    cmd = `gcc -shared -o $dlpath $objs $libs`
    try
        run(cmd)
        info("LINK --> $dlpath")
    catch
        println(cmd)
    end
end


text:
# This file was generated, do not modify it. # hide
coerce!(X, Finite => Count);


text:
@testset "Test generated structs" begin
    descriptor_file = joinpath(@__DIR__, "..", "src", "descriptors", "structs.json")
    existing_dir = joinpath(@__DIR__, "..", "src", "generated")
    output_dir = "tmp-test-generated-structs"
    if isdir(output_dir)
        rm(output_dir; recursive=true)
    end
    mkdir(output_dir)

    IS.generate_structs(descriptor_file, output_dir; print_results=false)

    matched = true
    try
        run(`diff $output_dir $existing_dir`)
    catch(err)
        @error "Generated structs do not match the descriptor file."
        matched = false
    finally
        rm(output_dir; recursive=true)
    end

    @test matched
end


text:
"""
    slab(data::AbstractData, h::Integer)

A "pancake" view into an underlying
data layout `data` at location `h`.
"""
function slab end
#=
# TODO: this could cause problems when it fails...
#slab(x::Number, inds...) = x
=#
slab(x, inds...) = x
slab(tup::Tuple, inds...) = map(x -> slab(x, inds...), tup)

"""
    column(data::AbstractData, i::Integer)

A contiguous "column" view into an underlying
data layout `data` at nodal point index `i`.
"""
function column end

column(x, inds...) = x
column(tup::Tuple, inds...) = map(x -> column(x, inds...), tup)


text:
function dden_complete(ygrid, W, eta, Z, mus, sig2; i, j)
  K = size(W, 2)
  L0 = size(eta[0], 3)
  L1 = size(eta[1], 3)
  L = Dict(0 => L0, 1 => L1)

  dden = [begin
            si = sqrt(sig2[i])
            dd = 0.0
            for k in 1:K
              ddk = 0.0
              z = Z[j, k]
              for ell in 1:L[z]
                mu_zl = mus[z][ell]
                eta_zijl = eta[z][i, j, ell]
                ddk += eta_zijl * pdf(Normal(mu_zl, si), y)
              end
              dd += W[i, k] * ddk
            end
            dd
          end for y in ygrid]

  return dden
end

function dden_complete(ygrid, W, eta, Z, mus, sig2)
  I, J = size(eta[0])[1:end-1]

  dden = [dden_complete(ygrid, W, eta, Z, mus, sig2, i=i, j=j)
          for i in 1:I, j in 1:J]

  return dden
end


text:
function gausslegendre(n::Integer)
    # GAUSSLEGENDRE(n) COMPUTE THE GAUSS-LEGENDRE NODES AND WEIGHTS IN O(n) time.

  if n <= 0
        Float64[], Float64[]
    elseif n == 1
        [0.0], [2.0]
    elseif n == 2
        [-1 / sqrt(3), 1 / sqrt(3)], [1.0, 1.0]
    elseif n == 3
        [-sqrt(3 / 5), 0.0, sqrt(3 / 5)], [5 / 9, 8 / 9, 5 / 9]
    elseif n == 4
        a = 2 / 7 * sqrt(6 / 5)
        ([-sqrt(3 / 7 + a), -sqrt(3/7-a), sqrt(3/7-a), sqrt(3/7+a)],
         [(18 - sqrt(30)) / 36, (18 + sqrt(30)) / 36,
          (18 + sqrt(30)) / 36, (18 - sqrt(30)) / 36])
    elseif n == 5
        b = 2 * sqrt(10 / 7)
        ([-sqrt(5 + b) / 3, -sqrt(5 - b) / 3, 0.0,
          sqrt(5 - b) / 3, sqrt(5 + b) / 3],
         [(322 - 13 * sqrt(70)) / 900, (322 + 13 * sqrt(70)) / 900, 128 / 225,
          (322 + 13 * sqrt(70)) / 900, (322 - 13 * sqrt(70)) / 900])
    elseif n <= 60
        # NEWTON'S METHOD WITH THREE-TERM RECURRENCE:
        rec(n)
    else
        # USE ASYMPTOTIC EXPANSIONS:
        asy(n)
    end
end

function asy(n)
    # COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING ASYMPTOTIC EXPANSIONS.
    # COMPLEXITY O(n).

    # Nodes and weights:
    m = (n + 1) >> 1
    a = besselZeroRoots(m)
    rmul!(a, 1 / (n + 0.5))
    x = legpts_nodes(n, a)
    w = legpts_weights(n, m, a)
    # Use symmetry to get the others:
    resize!(x, n)
    resize!(w, n)
    @inbounds for i in 1:m
        x[n + 1 - i] = x[i]
        w[n + 1 - i] = w[i]
    end
    @inbounds for i in 1:m
        x[i] = -x[i]
    end
    @inbounds mod(n, 2) == 1 && (x[m] = 0.0)
    x, w
end

function legpts_nodes(n, a)
    # ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE NODES.
    vn = 1 / (n + 0.5)
    m = length(a)
    nodes = cot.(a)
    vn² = vn * vn
    vn⁴ = vn² * vn²
    @inbounds if n <= 255
        vn⁶ = vn⁴ * vn²
        for i in 1:m
            u = nodes[i]
            u² = u^2
            ai = a[i]
            ai² = ai * ai
            ai³ = ai² * ai
            ai⁵ = ai² * ai³
            node = ai + (u - 1 / ai) / 8 * vn²
            v1 = (6 * (1 + u²) / ai + 25 / ai³ - u * muladd(31, u², 33)) / 384
            v2 = u * @evalpoly(u², 2595 / 15360, 6350 / 15360, 3779 / 15360)
            v3 = (1 + u²) * (-muladd(31 / 1024, u², 11 / 1024) / ai +
                             u / 512 / ai² + -25 / 3072 / ai³)
            v4 = (v2 - 1073 / 5120 / ai⁵ + v3)
            node = muladd(v1, vn⁴, node)
            node = muladd(v4, vn⁶, node)
            nodes[i] = node
        end
    elseif n <= 3950
        for i in 1:m
            u = nodes[i]
            u² = u^2
            ai = a[i]
            ai² = ai * ai
            ai³ = ai² * ai
            node = ai + (u - 1 / ai) / 8 * vn²
            v1 = (6 * (1 + u²) / ai + 25 / ai³ - u * muladd(31, u², 33)) / 384
            node = muladd(v1, vn⁴, node)
            nodes[i] = node
        end
    else
        for i in 1:m
            u = nodes[i]
            ai = a[i]
            node = ai + (u - 1 / ai) / 8 * vn²
            nodes[i] = node
        end
    end
    @inbounds for jj = 1:m
        nodes[jj] = cos(nodes[jj])
    end
    nodes
end

function legpts_weights(n, m, a)
    # ASYMPTOTIC EXPANSION FOR THE GAUSS-LEGENDRE WEIGHTS.
    vn = 1 / (n + 0.5)
    vn² = vn^2
    weights = Array{Float64}(undef, m)
    if n <= 850000
        @inbounds for i in 1:m
            weights[i] = cot(a[i])
        end
    end
    # Split out the part that can be vectorized by llvm
    @inbounds if n <= 170
        for i in 1:m
            u = weights[i]
            u² = u^2
            ai = a[i]
            ai⁻¹ = 1 / ai
            ai² = ai^2
            ai⁻² = 1 / ai²
            ua = u * ai
            W1 = muladd(ua - 1, ai⁻², 1.0) / 8
            W2 = @evalpoly(ai⁻², @evalpoly(u², -27.0, -84.0, -56.0),
                           muladd(-3.0, muladd(u², -2.0, 1.0), 6 * ua),
                           muladd(ua, -31.0, 81.0)) / 384
            W3 = @evalpoly(ai⁻¹, @evalpoly(u², 153 / 1024, 295 / 256, 187 / 96,
                                           151 / 160),
                           @evalpoly(u², -65 / 1024, -119 / 768, -35 / 384) * u,
                           @evalpoly(u², 5 / 512, 15 / 512, 7 / 384),
                           muladd(u², 1 / 512, -13 / 1536) * u,
                           muladd(u², -7 / 384, + 53 / 3072),
                           3749 / 15360 * u, -1125 / 1024)
            weights[i] = @evalpoly(vn², 1 / vn² + W1, W2, W3)
        end
    elseif n <= 1500
        for i in 1:m
            u = weights[i]
            u² = u^2
            ai = a[i]
            ai² = ai^2
            ai⁻² = 1 / ai²
            ua = u * ai
            W1 = muladd(ua - 1, ai⁻², 1.0) / 8
            W2 = @evalpoly(ai⁻², @evalpoly(u², -27.0, -84.0, -56.0),
                           muladd(-3.0, muladd(u², -2.0, 1.0), 6 * ua),
                           muladd(ua, -31.0, 81.0)) / 384
            weights[i] = muladd(vn², W2, 1 / vn² + W1)
        end
    elseif n <= 850000
        for i in 1:m
            u = weights[i]
            u² = u^2
            ai = a[i]
            ai² = ai^2
            ai⁻² = 1 / ai²
            ua = u * ai
            W1 = muladd(ua - 1, ai⁻², 1.0) / 8
            weights[i] = 1 / vn² + W1
        end
    else
        for i in 1:m
            weights[i] = 1 / vn²
        end
    end
    bJ1 = besselJ1(m)
    @inbounds for i in 1:m
        weights[i] = 2 / (bJ1[i] * (a[i] / sin(a[i])) * weights[i])
    end
    return weights
end

function rec(n)
    # COMPUTE GAUSS-LEGENDRE NODES AND WEIGHTS USING NEWTON'S METHOD.
    # THREE-TERM RECURENCE IS USED FOR EVALUATION. COMPLEXITY O(n^2).
    # Initial guesses:
    x0 = asy(n)[1]
    x = x0[1:n ÷ 2 + 1]
    # Perform Newton to find zeros of Legendre polynomial:
    PP1, PP2 = innerRec(n, x)
    @inbounds @simd for i in 1:length(x)
        x[i] -= PP1[i] / PP2[i]
    end
        # One more Newton for derivatives:
    PP1, PP2 = innerRec(n, x)
    @inbounds @simd for i in 1:length(x)
        x[i] -= PP1[i] / PP2[i]
    end

    # Use symmetry to get the other Legendre nodes and weights:
    m = length(x)
    resize!(x, n)
    resize!(PP2, n)
    @inbounds for i in 1:m-1
        x[n + 1 - i] = -x[i]
        PP2[n + 1 - i] = -PP2[i]
    end
    @inbounds for i in 1:n
        PP2[i] = 2 / ((1 - x[i]^2) * PP2[i]^2)
    end
    x,PP2
end

function innerRec(n, x)
    # EVALUATE LEGENDRE AND ITS DERIVATIVE USING THREE-TERM RECURRENCE RELATION.
    N = size(x, 1)
    myPm1 = Array{Float64}(undef, N)
    myPPm1 = Array{Float64}(undef, N)
    @inbounds for j = 1:N
        xj = x[j]
        Pm2 = 1.0
        Pm1 = xj
        PPm1 = 1.0
        PPm2 = 0.0
        for k = 1:(n - 1)
            Pm2, Pm1 = Pm1, muladd((2 * k + 1) * Pm1, xj, - k * Pm2) / (k + 1)
            PPm2, PPm1 = PPm1, ((2 * k + 1) * muladd(xj, PPm1, Pm2) -
                                k * PPm2) / (k + 1)
        end
        myPm1[j] = Pm1
        myPPm1[j] = PPm1
    end
    return myPm1, myPPm1
end

const besselZeros_20 = [2.4048255576957728, 5.5200781102863106,
                        8.6537279129110122, 11.791534439014281,
                        14.930917708487785, 18.071063967910922,
                        21.211636629879258, 24.352471530749302,
                        27.493479132040254, 30.634606468431975,
                        33.775820213573568, 36.917098353664044,
                        40.058425764628239, 43.199791713176730,
                        46.341188371661814, 49.482609897397817,
                        52.624051841114996, 55.765510755019979,
                        58.906983926080942, 62.048469190227170]

function besselZeroRoots(m)
    # BESSEL0ROOTS ROOTS OF BESSELJ(0,x). USE ASYMPTOTICS.
    # Use McMahon's expansion for the remainder (NIST, 10.21.19):
    jk = Array{Float64}(undef, m)
    p = (1071187749376 / 315, 0.0, -401743168 / 105, 0.0, 120928 / 15,
         0.0, -124 / 3, 0.0, 1.0, 0.0)
    # First 20 are precomputed:
    @inbounds for jj = 1:min(m, 20)
        jk[jj] = besselZeros_20[jj]
    end
    @inbounds for jj = 21:min(m, 47)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5], p[3])
    end
    @inbounds for jj = 48:min(m, 344)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5])
    end
    @inbounds for jj = 345:min(m,13191)
        ak = π * (jj - .25)
        ak82 = (.125 / ak)^2
        jk[jj] = ak + .125 / ak * muladd(ak82, p[7], 1.0)
    end
    @inbounds for jj = 13192:m
        ak = π * (jj - .25)
        jk[jj] = ak + .125 / ak
    end
    return jk
end

const besselJ1_10 = [0.2695141239419169, 0.1157801385822037,
                     0.07368635113640822, 0.05403757319811628,
                     0.04266142901724309, 0.03524210349099610,
                     0.03002107010305467, 0.02614739149530809,
                     0.02315912182469139, 0.02078382912226786]

function besselJ1(m)
    # BESSELJ1 EVALUATE BESSELJ(1,x)^2 AT ROOTS OF BESSELJ(0,x).
    # USE ASYMPTOTICS. Use Taylor series of (NIST, 10.17.3) and McMahon's
    # expansion (NIST, 10.21.19):
    Jk2 = Array{Float64}(undef, m)
    c = (-171497088497 / 15206400, 461797 / 1152, -172913 / 8064,
         151 / 80, -7 / 24, 0.0, 2.0)
    # First 10 are precomputed:
    @inbounds for jj = 1:min(m, 10)
        Jk2[jj] = besselJ1_10[jj]
    end
    @inbounds for jj = 11:min(m, 15)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3],
                                                  c[2], c[1]), ak2^2, c[7])
    end
    @inbounds for jj = 16:min(m, 21)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3], c[2]),
                                        ak2^2, c[7])
    end
    @inbounds for jj = 22:min(m,55)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3]),
                                        ak2^2, c[7])
    end
    @inbounds for jj = 56:min(m,279)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(muladd(ak2, c[4], c[5]), ak2^2, c[7])
    end
    @inbounds for jj = 280:min(m,2279)
        ak = π * (jj - .25)
        ak2 = (1 / ak)^2
        Jk2[jj] = 1 / (π * ak) * muladd(ak2^2, c[5], c[7])
    end
    @inbounds for jj = 2280:m
        ak = π * (jj - .25)
        Jk2[jj] = 1 / (π * ak) * c[7]
    end
    return Jk2
end


text:
module Juice

using Reexport

@reexport using LogicCircuits
@reexport using ProbabilisticCircuits
@reexport using DiscriminativeCircuits

end # module


text:
using SDC
import Plots
using LaTeXStrings
using PrettyTables
"""
Currently a demonstration of SDC solving
"""
function main()

    # 5-node quadrature
    n = 5
    # IMEX Euler, half explicit, half implicit (dalhquist equation)
    fE(t, u) = 0.5*u
    fI(t, u) = 0.5*u
    # Solve implicit euler analytically
    implicit_solver(α, t, b) = b/(1+α)

    # Create SDC timestepper with provided parameters (should change code
    # structure to avoid this)
    my_SDC = make_SDC(n, fE, fI, implicit_solver)

    # Analytic solution
    u(t) = exp(t)
    # Start time
    t = 0.0
    # Initial value
    u_n = u(t)
    # Timesteps
    Δtimes = [1.0, 1.0, 1.0, 1.0, 1.0]
    # Non-start times, derived from timesteps
    times = cumsum(Δtimes) .+ t

    # Perform timestepping
    u_actual = []
    u_approx = []
    for i in eachindex(Δtimes)
        t += Δtimes[i]
        u_n = my_SDC(u_n, t, Δtimes[i])
        push!(u_actual, u(t))
        push!(u_approx, u_n)
    end

    # Calculate errors
    u_diff_abs = abs.(u_actual - u_approx)
    u_diff_rel = [u_diff_abs[i]/u_actual[i] for i in eachindex(u_diff_abs)]

    # Produce table of results
    table_data = hcat(times, u_actual, u_approx, u_diff_abs, u_diff_rel)
    header = ["t", "u(t)", "est u(t)", "Abs Err", "Rel Err"]
    pretty_table(table_data; header=header)

    # Produce graph of results
    p1 = Plots.plot(times, u_diff_abs, yaxis=:log, title = "Absolute Error of e^t, t_0=0, variable timesteps",
                   legend=false)
    Plots.xlabel!(p1, "Time")
    Plots.ylabel!(p1, "Absolute Error")
    p2 = Plots.plot(times, u_diff_rel, yaxis=:log, title = "Relative Error of e^t, t_0=0, variable timesteps",
                   legend=false)
    Plots.xlabel!(p2, "Time")
    Plots.ylabel!(p2, "Relative Error")
    plt = Plots.plot(p1, p2, layout=(2,1))
    Plots.savefig(plt, "./dahlquist_test.png")
    return plt
end

display(main())


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

include("ordinal.jl")
export ordinal

include("rad2sec.jl")
export rad2sec

include("sec2rad.jl")
export sec2rad


text:
#=
using DelimiteFiles
=#

# NOTE: These simulation setting are identical to that of test-sims-6-5 EXCEPT
# that W_true is different. The W_true here has low abundance (< 10% each) for
# 2 pairs (i.e. 4) features. Thus, they should be harder to recover, and
# uncertainty in R_i should be greater.

# Simulation Name
simname = basename("$(@__DIR__)")
println(simname); flush(stdout)

# NOTE: write to scratchdir
function outdir_suffix(dataseed, mcmcseed, phi)
  return "dataseed$(dataseed)-mcmcseed$(mcmcseed)-phi$(phi)"
end

# NOTE: modify
settings = let
  Z = Matrix{Bool}(readdlm(joinpath(@__DIR__, "Z.txt"), Int, comments=true))
  W = readdlm(joinpath(@__DIR__, "W.txt"), comments=true)
  Nfac = 2000  # for Nfac=2000, wasn't able to recover feature 1.
  N = [1, 1] *  Nfac
  # NOTE: these are actually quite important.
  # Need pthin to be at least 4, but batchprop can be small, if N is
  # sufficiently large.
  pthin = 4
  batchprop = 0.05
  [Dict(:simname => simname,
        :repfam_dist_scale => phi,  # doesn't work for kmeans init when phi=1.0
        :N => N,
        :Z => Z,
        :W => W,
        :thin_samps => 2,
        :nburn => 500,
        :nsamps => 1000,
        :Lmcmc => Dict(0 => 2, 1 => 2),
        :Kmcmc => 10,
        :pthin => pthin,
        :batchprop => batchprop,
        :dataseed => dataseed,
        :mcmcseed => mcmcseed,
        :outdir_suffix => outdir_suffix(dataseed, mcmcseed, phi))
   for phi in [0.0, 1.0, 10.]  # pthin{2,10} didn't work for Nfac=2000
   for dataseed in 1:3  # batchprop{.05} didn't work for Nfac=2000
   for mcmcseed in 1:3]  # alpha{10} didn't work for Nfac=2000
end


text:
module SwaggerMarkdown

import JSONSchema
import YAML
import JSON
import UUIDs

const ROOT = dirname(dirname(@__FILE__))

include("OpenAPI.jl")
include("Utils.jl")

export OpenAPI, build, @swagger, validate_spec

const TMP = joinpath(ROOT, "tmp")

macro swagger(doc)
    if !isdir(TMP)
        mkdir(TMP)
    end
    return :( io = open(joinpath(SwaggerMarkdown.TMP, "$(string(UUIDs.uuid4())).yml"), "w"); write(io, $doc); close(io) )
end

"""
    build(spec::Union{OpenAPI, Dict{String, Any}})::Dict{String, Any}

Build a `Dict{String, Any}` as OpenAPI definitions.

# Arguments
Required:
- `spec::Union{OpenAPI, Dict{String, Any}}` : An `OpenAPI` or a `Dict{String, Any}}` with some pre-defined specs.
"""
function build(spec::Union{OpenAPI, Dict{String, Any}}; doc_parser::Function=parse_spec)::Dict{String, Any}
    if spec isa OpenAPI
        if isempty(spec.paths)
            spec.paths = doc_parser()
        end
        openApi = spec
        spec = Dict{String, Any}()
        version_name = openApi.version == "2.0" ? "swagger" : "openapi"
        spec[version_name] = openApi.version
        spec["info"] = openApi.info
        spec["paths"] = openApi.paths
        for (key, val) in openApi.optional_fields
            spec[key] = val
        end
    else # spec is Dict{String, Any} type
        if !haskey(spec, "paths")
            spec["paths"] = doc_parser()
        end
    end
    validate_spec(spec)
    return spec
end

"""
    build(file::String)::Dict{String, Any}

Build a `Dict{String, Any}()` as OpenAPI definitions from a `json` or `yml` file.

# Arguments
Required:
- `file::String` : A `json` or `yml` file used to build the OpenAPI definitions.
"""
function build(file::String; doc_parser::Function=parse_spec)::Dict{String, Any}
    if endswith(file, ".json")
        return build(JSON.parsefile(file), doc_parser=doc_parser)
    end
    if endswith(file, ".yml")
        return build(YAML.load_file(file, dicttype=Dict{String, Any}), doc_parser=doc_parser)
    end
end

end # module


text:
#=
given a symmetric sparsity pattern with some sparse entries unknown, find integer values
(within some bounds) for the unknowns such that the matrix is PSD-completable,
maximimizing the minimum eigenvalue

max y :
S + X - y I ∈ SparsePSD⋆
Xᵢⱼ ∈ ℤ, -M ≤ Xᵢⱼ ≤ M, ∀ (i,j) unknown

where:
- sparse parameter matrix S contains is the known entries
- sparse variable matrix X contains the unknown entries
- M is a bound on the absolute value of all unknown entries
=#

struct CompletablePSD <: ExampleInstance
    d::Int # side dimension
    sparsity::Float64 # fraction of elements in the pattern
    frac_known::Float64 # fraction of sparse entries known
    use_nat::Bool # use sparse PSD cone, else PSD cone formulation
end

function build(inst::CompletablePSD)
    d = inst.d
    @assert d >= 2

    # generate underlying positive definite matrix S0
    S0 = sprand(Bool, d, d, inst.sparsity)
    S0 = Matrix{Float64}(S0 * S0' + I)
    pattern = tril!(sprand(Bool, d, d, inst.sparsity) + I)
    (rows, cols, _) = findnz(pattern)
    S0vals = [S0[r, c] for (r, c) in zip(rows, cols)]
    num_sparse = length(S0vals)

    # select which entries of S0 will be known/unknown
    known_pattern = sprand(Bool, num_sparse, inst.frac_known)
    known_pattern[1] = known_pattern[end] = true
    num_unknown = num_sparse - sum(known_pattern)
    @assert num_unknown >= 1
    M = maximum(abs, S0vals[.!known_pattern])

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, y)
    JuMP.@objective(model, Max, y)

    JuMP.@variable(model, x[1:num_unknown], Int)
    JuMP.@constraint(model, -M .<= x .<= M)

    aff = zeros(JuMP.AffExpr, num_sparse)
    k = 1
    scal = (inst.use_nat ? rt2 : 1.0)
    for i in eachindex(known_pattern)
        if known_pattern[i]
            aff[i] = S0vals[i]
        else
            aff[i] = x[k]
            k += 1
        end
        if rows[i] == cols[i]
            aff[i] = aff[i] - y
        else
            aff[i] = scal * aff[i]
        end
    end
    @assert k == num_unknown + 1

    if inst.use_nat
        KT = Hypatia.PosSemidefTriSparseCone{Hypatia.Cones.PSDSparseDense, Float64, Float64}
        K = KT(d, rows, cols, true)
        JuMP.@constraint(model, aff in K)
    else
        mat = Matrix(sparse(rows, cols, aff))
        for i in 1:d, j in 1:i
            if iszero(mat[i, j])
                mat[i, j] = JuMP.@variable(model)
            end
        end
        LinearAlgebra.copytri!(mat, 'L')
        JuMP.@constraint(model, Symmetric(mat) in JuMP.PSDCone())
    end

    # save for use in tests
    model.ext[:x] = x
    model.ext[:S0_eigmin] = eigmin(Hermitian(S0))

    return model
end

function test_extra(inst::CompletablePSD, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    x = JuMP.value.(model.ext[:x])
    @test x ≈ round.(Int, x) atol = tol rtol = tol
    # check eigmin(S + X) >= eigmin(S0)
    S0_eigmin = model.ext[:S0_eigmin]
    @test JuMP.objective_value(model) >= S0_eigmin - tol
    return
end


text:
using L1Homotopy
using Base.Test

# write your own tests here
@testset "Homotopy" begin
    N = 100
    A = randn(N,N)
    x_in = Float64.(rand(N) .> 0.9)
    x_out = L1_homotopy(A, A*x_in, 0.1)
    # Solution is close but approximate on the support
    @test x_in ≈ x_out
    # Solution should be exactly equal to zero off the support
    @test (x_out .== 0) == (x_in .== 0)

end

@testset "Random matrix special cases" begin
    A = [-1.04 -0.70 -0.82  0.58 -0.01;
          1.04  1.29 -0.57 -0.02 -1.02;
         -0.29 -0.02  0.35 -0.22  0.98;
         -1.56  1.13  0.49 -0.23 -0.26;
         -0.72  1.39 -0.28  0.75  0.07]

    x_in = [0.0, 0.0, 1.0, 1.0, 0.0]

    x_out = L1_homotopy(A, A*x_in, 0.001)

    @test x_in ≈ x_out
    # This test doesn't yet pass, because the maintenance of the active set
    # doesn't handle the case where γ⁺ ≈ γ⁻ to numerical precision.
    @test_broken (x_out .== 0) == (x_in .== 0)
end


text:
function eigsolve_solver(; kwargs...)
  function solver(H, t, psi0; kws...)
    howmany = 1
    which = get(kwargs, :solver_which_eigenvalue, :SR)
    solver_kwargs = (;
      ishermitian=get(kwargs, :ishermitian, true),
      tol=get(kwargs, :solver_tol, 1E-14),
      krylovdim=get(kwargs, :solver_krylovdim, 3),
      maxiter=get(kwargs, :solver_maxiter, 1),
      verbosity=get(kwargs, :solver_verbosity, 0),
    )
    vals, vecs, info = eigsolve(H, psi0, howmany, which; solver_kwargs...)
    psi = vecs[1]
    return psi, info
  end
  return solver
end

function dmrg(H, psi0::MPS; kwargs...)
  t = Inf # DMRG is TDVP with an infinite timestep and no reverse step
  reverse_step = false
  psi = tdvp(eigsolve_solver(; kwargs...), H, t, psi0; reverse_step, kwargs...)
  return psi
end


text:
module EventUtils

export run_day, make_my_day, download_input

using Dates: day, today
using InteractiveUtils: edit
using Test: @test, @testset
using UnPack: @unpack

import ..Advent2020
import HTTP

include("download_data.jl")
include("make_or_run.jl")

end


text:
# NewPrefix(prefix_symbol, power_of_ten)
# Be careful not to define anything that conflicts with units
const Yotta = Prefix("Y",24)
const Zetta = Prefix("Z",21)
const Exa = Prefix("E",18)
const Peta = Prefix("P",15)
const Tera = Prefix("T",12)
const Giga = Prefix("G",9)
const Mega = Prefix("M",6)
const Kilo = Prefix("k",3)
const Hecto = Prefix("h",2)
const Deka = Prefix("da",1)
const Deci = Prefix("d",-1)
             Prefix("",0)
const Centi = Prefix("c",-2)
const Milli = Prefix("m",-3)
const Micro = Prefix("µ",-6)
const Nano = Prefix("n",-9)
const Pico = Prefix("p",-12)
const Femto = Prefix("f",-15)
const Atto = Prefix("a",-18)
const Zepto = Prefix("z",-21)
const Yocto = Prefix("y",-24)

export Yotta, Zetta, Exa, Peta, Tera, Giga, Mega, Kilo, Hecto, Deka, Deci, Centi, Milli
export Micro, Nano, Pico, Femto, Atto, Zepto, Yocto


text:
module ILKBackend

using LinearAlgebra

MatrixType = Matrix{Float64}

function matrix(r,c)
    return zeros(Float64,r,c)
end

function pose()
    return MatrixType(I,4,4)
end

function setPosition(pose, x,y,z)
    pose[1,4] = x
    pose[2,4] = y
    pose[3,4] = z
end

function setRotation(pose, xx, xy, xz, yx, yy, yz, zx, zy, zz)
    pose[1,1] = xx
    pose[1,2] = xy
    pose[1,3] = xz
    pose[2,1] = yx
    pose[2,2] = yy
    pose[2,3] = yz
    pose[3,1] = zx
    pose[3,2] = zy
    pose[3,3] = zz
end

function positionView(pose)
    return view(pose, 1:3, 4)
end

function rotationView(pose)
    return view(pose, 1:3, 1:3)
end

function zaxisView(pose)
    return view(pose, 1:3, 3)
end

function rot_z__a_x_b(angle)
    mx = zeros(4,4)
    s = sin(angle)
    c = cos(angle)

    mx[1,1] = c
    mx[1,2] = -s
    mx[2,1] = s
    mx[2,2] = c
    mx[3,3] = 1
    mx[4,4] = 1
    return mx
end

function rot_z__b_x_a(angle)
    mx = zeros(4,4)
    s = sin(angle)
    c = cos(angle)

    mx[1,1] = c
    mx[1,2] = s
    mx[2,1] = -s
    mx[2,2] = c
    mx[3,3] = 1
    mx[4,4] = 1
    return mx
end

function tr_z__a_x_b(length)
    mx = pose()
    mx[3,4] = length
    return mx
end

function tr_z__b_x_a(length)
    mx = pose()
    mx[3,4] = -length
    return mx
end

function linearCoords(_6dvector)
    return view(_6dvector, 4:6, 1)
end

function angularCoords(_6dvector)
    return view(_6dvector, 1:3, 1)
end

function geometricJacobianColumn_revolute(poi, jointOrigin, jointAxis, column)
    angularCoords(column)[:] = jointAxis
     linearCoords(column)[:] = LinearAlgebra.cross(jointAxis, poi - jointOrigin )
end

function geometricJacobianColumn_prismatic(jointAxis, column)
    angularCoords(column)[:] = zeros(3,1)
     linearCoords(column)[:] = jointAxis
end

mutable struct ik_pos_dbg
    iter_count::UInt64
    actual_pos::MatrixType
    actual_or::MatrixType
end

struct ik_pos_cfg
    dt::Float64
    eps_pos_err_norm::Float64
    eps_or_err_norm::Float64
    max_iter::UInt64
    ls_damping::Float64
end

struct AxisAngle
    axis
    angle::Float64
end

function omega(aa::AxisAngle)
    return aa.axis*aa.angle
end

# The difference between two rotation matrices, as an axis-angle.
# @return the rotation required to go from the second to the first argument,
#         that is, the first "minus" the second.

function orientationDistance(_R_desired, _R_actual)
    thresh = 1e-6   # TODO #magic-number

    # TODO checks on the size 3x3
    R = transpose(_R_actual) *  _R_desired   #  this is 'actual_R_desired'

    x = R[3,2] - R[2,3]
    y = R[1,3] - R[3,1]
    z = R[2,1] - R[1,2]
    norm = sqrt(x*x + y*y + z*z)

    if norm < thresh
        return AxisAngle([0.0; 0.0; 1.0], 0.0) # arbitrary choice of (0,0,1) axis
    end

    theta = atan( norm, LinearAlgebra.tr(R)-1 )
    return AxisAngle([x/norm, y/norm, z/norm], theta)
    # TODO check corner cases, theta close to 0/PI, bad numerical
end

function leastSquaresSolve(A, b)
    fact = LinearAlgebra.svd(A)
    return (fact.V * diagm(0 => (1 ./ fact.S)) ) * transpose(fact.U) * b
end

function ct_twist(newframe_H_current, twist_current, twist_newframe)
    ## We need to perform a spatial motion vector coordinate transform, given
    ## the homogeneous coordinate transform for the same frames.
    ## Check e.g. chapter 2 of Roy's RBDA, for the relation between the two
    ## transforms, to understand the last line of this function

    R = rotationView( newframe_H_current );      # 3x3 rotation matrix
    w = angularCoords( twist_current );          # omega in the current coordinates
    angularCoords( twist_newframe )[:] = R * w;  # omega in the new coordinates
     linearCoords( twist_newframe )[:] = R * linearCoords( twist_current ) +
                           cross( positionView( newframe_H_current),
                                  angularCoords( twist_newframe )    );
end

end # module

max_stars_repo_name: ESROCOS/ilk-compiler

max_issues_repo_name: ESROCOS/ilk-compiler

max_forks_repo_name: ESROCOS/ilk-compiler


text:
using Pkg

Pkg.activate(dirname(@__FILE__))
Pkg.instantiate()

using Kuber
ctx = Kuber.KuberContext()

kubqueens_pod = kuber_obj(ctx, """{
                  "kind": "Pod",
                  "metadata":{
                      "name": "kubqueens-pod",
                      "namespace": "default",
                      "labels": {
                          "name": "kubqueens-pod"
                      }
                  },
                  "spec": {
                      "containers": [{
                          "name": "kubqueens",
                          "image": "hojpreg.azurecr.io/kubqueens",
                          "ports": [{"containerPort": 8888}]
                      }]
                  }
              }""")

kubqueens_service = kuber_obj(ctx, """{
                  "kind": "Service",
                  "metadata": {
                      "name": "kubqueens-service",
                      "namespace": "default",
                      "labels": {"name": "kubqueens-service"}
                  },
                  "spec": {
                      "type": "LoadBalancer",
                      "ports": [{"port": 8888}],
                      "selector": {"name": "kubqueens-pod"}
                  }
              }""")

put!(ctx, kubqueens_pod)
put!(ctx, kubqueens_service)

while true
    println("waiting for loadbalancer to be configured...")
    status = get(ctx, :Service, "kubqueens-service").status
    if nothing !== status.loadBalancer.ingress && !isempty(status.loadBalancer.ingress)
        println(status.loadBalancer.ingress[1])
        break
    end
    sleep(5)
end

# To clean up the services you can use the commented commands below.

println("To clean up the services you can use the commented commands below:")

println("delete!(ctx, kubqueens_service)")
println("delete!(ctx, kubqueens_pod)")

max_stars_repo_path: Chapter 14/kubernetes/kubserver.jl

max_issues_repo_path: Chapter 14/kubernetes/kubserver.jl

max_forks_repo_path: Chapter 14/kubernetes/kubserver.jl


text:
abstract type AbstractCommand end

struct SetObject{O <: AbstractObject} <: AbstractCommand
    object::O
    path::Path
end

struct SetTransform{T <: Transformation} <: AbstractCommand
    tform::T
    path::Path
end

struct Delete <: AbstractCommand
    path::Path
end

struct SetProperty{T} <: AbstractCommand
    path::Path
    property::String
    value::T
end

struct SetAnimation{A <: Animation} <: AbstractCommand
    animation::A
    play::Bool
    repetitions::Int
end

SetAnimation(anim::Animation; play=true, repetitions=1) = SetAnimation(anim, play, repetitions)


text:
type Line
    rect::RectangleShape
end

function Line(p1::Vector2, p2::Vector2, thickness::Real=2)
    rect = RectangleShape()
    set_position(rect, to_vec2f(p1))
    set_size(rect, Vector2f(distance(p1, p2), thickness))
    set_rotation(rect, rad2deg(atan2(p2.y - p1.y, p2.x - p1.x)))
    set_origin(rect, Vector2f(0, thickness / 2))

    Line(rect)
end

function copy(l::Line)
    Line(l.p1, l.p2, l.thickness)
end

function set_points(l::Line, p1::Vector2, p2::Vector2)
    set_position(l.rect, to_vec2f(p1))
    set_rotation(l.rect, rad2deg(atan2(p2.y - p1.y, p2.x - p1.x)))
    set_size(l.rect, Vector2f(distance(p1, p2), get_thickness(l)))
end

function get_points(l::Line)
    pos = get_position(l.rect)
    rot = get_rotation(l.rect)
    len = get_size(l.rect).x

    p2 = Vector2f(pos.x + len*cosd(rot), pos.y + len*sind(rot))
    pos, p2
end

function set_thickness(l::Line, thickness::Real)
    p1, p2 = get_points(l)
    set_size(l.rect, Vector2f(distance(p1, p2), thickness))
    set_origin(l.rect, Vector2f(0, thickness / 2))
end

function get_thickness(l::Line)
    get_size(l.rect).y
end

function set_fillcolor(l::Line, c::Color)
    set_fillcolor(l.rect, c)
end

function get_fillcolor(l::Line)
    get_fillcolor(l.rect)
end

function set_outlinecolor(l::Line, c::Color)
    set_outlinecolor(l.rect, c)
end

function get_outlinecolor(l::Line)
    get_outlinecolor(l.rect)
end

function set_outline_thickness(l::Line, thickness::Real)
    set_outline_thickness(l.rect, thickness)
end

function get_outline_thickness(l::Line)
    get_outline_thickness(l.rect)
end

function draw(window::RenderWindow, object::Line)
    draw(window, object.rect)
end

function draw(texture::RenderTexture, object::Line)
    draw(texture, object.rect)
end


text:
using Distributions
dist = TriangularDist(0,2,1)

println("Parameters: \t\t\t",params(dist))
println("Central descriptors: \t\t",mean(dist),"\t",median(dist))
println("Dispersion descriptors: \t", var(dist),"\t",std(dist))
println("Higher moment shape descriptors: ",skewness(dist),"\t",kurtosis(dist))
println("Range: \t\t\t\t", minimum(dist),"\t",maximum(dist))
println("Mode: \t\t\t\t", mode(dist), "\tModes: ",modes(dist))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook