text:
#########################################################################################
### Some functions to print output and monitor the progress and termination status of CALiPPSO
#########################################################################################
#Memory usage in GB. It's not used by CALiPPSO, but it might be useful at some point to monitor the resources consumed during the execution of the main function
memory_usage() = round(parse(Int, split(read(`ps -p $(getpid()) -o rss`, String))[2]) / 1024^2, digits=3)

"Print the state of CALiPPSO maximum displacement and inflation factor when maximum number of iterations is exceeded"
function print_failed_max_iters(max_iters::Int64, sqrΓ, tol_conv_Γ, max_Si, tol_S_conv; color::Symbol=:red)
    printstyled("Maximum number of iterations (=$max_iters) reached and failed to produce a jammed packing!!\n", bold=true, color=color)
    println("Convergence criterion for √Γ-1: ", tol_conv_Γ, "\t and for max s_i: ", tol_S_conv)
    printstyled("Final values before stopping: \t √Γ-1 = ", sqrΓ-1, ", \t max s_i = ", max_Si, "\n\n", bold=true)
end

"Print information about CALiPPSO progress and current status of the system whenever a non-isostatic configuration is created."
function print_non_isostatic(d::Int64, zs::Vector{Int64}, non_rattlers::Vector{Int64}, t::Int64, max_Si::Float64, bound_s::Float64, f_mismatch::Float64, small_fs::Vector{Float64})

    println("Iteration: ", t)
    printstyled("Non-isostatic packing!!!\n", color=:yellow, bold=true)

    Nc = 0.5*sum(zs[non_rattlers]);
    N = length(zs); Nnr = length(non_rattlers)
    rattlers = setdiff(1:N, non_rattlers); Ncr = 0.5*sum(zs[rattlers])
    zs_rats = zs[rattlers]; inds_max_zs = sortperm(zs_rats, rev=true)[1:min(10, length(zs_rats))]
    iso_gap = d*(Nnr-1)+1 - Nc

    println("N_non-rattlers = ", Nnr, "\t N_c [stable particles] = ", Nc, "\t Isostaticity gap = ", iso_gap, "\t N_c [rattlers] = ", Ncr)
    println("Amount of rattlers = ", length(rattlers))
    println("(rattlers, z_rattlers) with contacts = ", [rattlers[inds_max_zs] zs_rats[inds_max_zs]])

    printstyled("max |sᵢ| (stable particles) = ", max_Si, "; \t bound on |sᵢ| = ", bound_s, "\t Difference of these quantities = ", max_Si - bound_s , "\n", bold=true)
    printstyled("----------------------------------\n", color=:yellow)
    println("Force mismatch = ", f_mismatch)
    println("Smallest forces = ", small_fs)

    printstyled("__________________________________________________________________________\n\n", color=:yellow)
end

function print_non_isostatic(d::Int64, zs::Vector{Int64}, non_rattlers::Vector{Int64})

    printstyled("Non-isostatic packing!!!\n", color=:yellow, bold=true)

    Nc = 0.5*sum(zs[non_rattlers])
    N = length(zs); Nnr = length(non_rattlers)
    rattlers = setdiff(1:N, non_rattlers); Ncr = 0.5*sum(zs[rattlers])
    zs_rats = zs[rattlers]; inds_max_zs = sortperm(zs_rats, rev=true)[1:min(10, length(zs_rats))]
    iso_gap = d*(Nnr-1)+1 - Nc

    println("N_non-rattlers = ", Nnr, "\t N_c [stable particles] = ", Nc, "\t Isostaticity gap = ", iso_gap, "\t N_c [rattlers] = ", Ncr)
    println("Amount of rattlers = ", length(rattlers))
    println("(rattlers, z_rattlers) with contacts = ", [rattlers[inds_max_zs] zs_rats[inds_max_zs]])
    printstyled("__________________________________________________________________________\n\n", color=:yellow)
end

"Print info about CALiPPSO progress, including new values of ``|s_i|``, density and ``R``, statistics of constraints and contacts, forces, etc."
function print_monitor_progress(max_Si::T, R::T, N::I, L::T, d::I, zs::Vector{I}, non_rattlers::Vector{I}, n_constr::Vector{I}, iso_cond::Bool, f_mismatch::T, small_fs::Vector{T};
    dig_R::I=9, dig_z::I = Int(ceil(log10(N*d))), dig_S::I=13, color::Symbol=:cyan) where {I<:Int, T<:AbstractFloat}

     ϕ = packing_fraction(d, R, N, L)

     Nnr = length(non_rattlers)
     rattlers = setdiff(1:N, non_rattlers)
     zs_rat = zs[rattlers]

     # Statistics number of constraints
     max_constr = maximum(n_constr); avg_constr = round(mean(n_constr), digits=2); std_constr = round(std(n_constr), digits=2)

     # Parameters for printing output
     max_Si_print = round(max_Si, digits = dig_S)

    printstyled("\t Max |sᵢ| (stable particles) = ", max_Si_print ,"\n New  values of (φ, R) = ", round.([ϕ, R],  digits=dig_R), "\n", bold=true)

     if Nnr>0
        zs_non_rat = zs[non_rattlers]

        if iso_cond
            println("(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr], "\t Isostatic: ", iso_cond, "\t (Mean, Max, total in rattlers) z = ", [round(mean(zs_non_rat), digits=dig_z), maximum(zs), sum(zs_rat) ] , "\t\tNon_rattlers= ", Nnr)
        else
            iso_gap = (d*(Nnr-1)+1) - 0.5*sum(zs_non_rat)
            println("(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr], "\t Isostatic: ", iso_cond, ";  isostaticity gap = ", iso_gap, "\t (Mean, Max, total in rattlers) z = ", [round(mean(zs_non_rat), digits=dig_z), maximum(zs), sum(zs_rat) ] , "\t\tNon-rattlers= ", Nnr)
        end

     else
        printstyled("All particles are rattlers!!\n", color=:yellow)
        println("Max constraints per particle per d: ", max_constr, "\t Isostatic: NA", "\tNon-rattlers= ", 0, "\tTotal z of rattlers = ", sum(zs_rat))
     end

    printstyled("----------------------------------\n", color=color)
    println("Force mismatch = ", f_mismatch)
    println("Sample of smallest forces = \t", small_fs)

    printstyled("__________________________________________________________________________\n\n", color=color)
end

"Print info about CALiPPSO progress, including new values of ``|s_i|``, density and some statistics of ``R``, statistics of constraints and contacts, forces, etc."
function print_monitor_progress(max_Si::T, Rs::Vector{T}, L::T, d::I, zs::Vector{I}, non_rattlers::Vector{I}, n_constr::Vector{I}, iso_cond::Bool, f_mismatch::T, small_fs::Vector{T};
     dig_R::I=9, dig_z::I = Int(ceil(log10(length(Rs)*d))), dig_S::I=13, color::Symbol=:cyan) where {I<:Int, T<:AbstractFloat}

     ϕ = packing_fraction(d, Rs, L)
     N = length(Rs)

     Nnr = length(non_rattlers)
     rattlers = setdiff(1:N, non_rattlers)
     zs_rat = zs[rattlers]

     # Statistics number of constraints
     max_constr = maximum(n_constr); avg_constr = round(mean(n_constr), digits=2); std_constr = round(std(n_constr), digits=2)

     # Parameters for printing output
     max_Si_print = round(max_Si, digits = dig_S)

     Rmin = minimum(Rs); Ravg = mean(Rs); Rmax = maximum(Rs)
     printstyled("\t Max |sᵢ| (stable particles) = ", max_Si_print ,"\n New  values of φ = ", round(ϕ, digits=dig_R)," and (Rₘᵢₙ, ⟨R⟩, Rₘₐₓ) = ", round.([Rmin, Ravg, Rmax],  digits=dig_R), "\n", bold=true)

     if Nnr>0
        zs_non_rat = zs[non_rattlers]

        if iso_cond
            println("(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr], "\t Isostatic: ", iso_cond, "\t (Mean, Max, total in rattlers) z = ", [round(mean(zs_non_rat), digits=dig_z), maximum(zs), sum(zs_rat) ] , "\t\tNon_rattlers= ", Nnr)
        else
            iso_gap = (d*(Nnr-1)+1) - 0.5*sum(zs_non_rat)
            println("(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr], "\t Isostatic: ", iso_cond, ";  isostaticity gap = ", iso_gap, "\t (Mean, Max, total in rattlers) z = ", [round(mean(zs_non_rat), digits=dig_z), maximum(zs), sum(zs_rat) ] , "\t\tNon-rattlers= ", Nnr)
        end

     else
        printstyled("All particles are rattlers!!\n", color=:yellow)
        println("Max constraints per particle per d: ", max_constr, "\t Isostatic: NA", "\tNon-rattlers= ", 0, "\tTotal z of rattlers = ", sum(zs_rat))
     end

    printstyled("----------------------------------\n", color=color)
    println("Force mismatch = ", f_mismatch)
    println("Sample of smallest forces = \t", small_fs)

    printstyled("__________________________________________________________________________\n\n", color=color)
end

"Print that CALiPPSO has converged, and some basic info (final ϕ, number of non rattlers, etc.)."
function print_converged(t::I, status, sqrΓ::T, max_Si::T, R::T, N::I, L::T, d::I, n_constr::Vector{I}, Nnr::I; dig_R::I=8, dig_S::I=15, color::Symbol=:green) where {I<:Int, T<:AbstractFloat}
    ϕ = packing_fraction(d, R, N, L)

     # Statistics number of constraints
     max_constr = maximum(n_constr); avg_constr = round(mean(n_constr), digits=2); std_constr = round(std(n_constr), digits=2)

    printstyled("\n\n\tCALiPPSO converged! \t \t Status of last LP optimization: ", status, "\n", bold=true, color=color)
    printstyled("Iterations to convergence = ", t, ",\t√Γ-1 = ", sqrΓ-1, ",\t Max displacement = ", round(max_Si, digits=dig_S), ",\t (ϕ, R) = ",  round.([ϕ, R], digits=dig_R), "\n", bold=true, color=color)
    println("\tNon_rattlers= ", Nnr, "\t% of rattlers = ", round(100*(1-Nnr/N), digits=3),  "\t(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr])
end

function print_converged(t::I, status, sqrΓ::T, max_Si::T, Rs::Vector{T}, N::I, L::T, d::I, n_constr::Vector{I}, Nnr::I; dig_R::I=8, dig_S::I=15, color::Symbol=:green) where {I<:Int, T<:AbstractFloat}
    ϕ = packing_fraction(d, Rs, L)
    # Statistics of radii
    Rmin = minimum(Rs); Ravg = mean(Rs); Rmax = maximum(Rs)

     # Statistics number of constraints
     max_constr = maximum(n_constr); avg_constr = round(mean(n_constr), digits=2); std_constr = round(std(n_constr), digits=2)

    printstyled("\n\n\tCALiPPSO converged! \t \t Status of last LP optimization: ", status, "\n", bold=true, color=color)
    printstyled("Iterations to convergence = ", t, ",\t√Γ-1 = ", sqrΓ-1, ",\t Max displacement = ", round(max_Si, digits=dig_S), ",\t Final  values of φ = ", round(ϕ, digits=dig_R), " and (Rₘᵢₙ, ⟨R⟩, Rₘₐₓ) = ", round.([Rmin, Ravg, Rmax],  digits=dig_R), "\n", bold=true, color=color)
    println("\tNon_rattlers= ", Nnr, "\t% of rattlers = ", round(100*(1-Nnr/N), digits=3),  "\t(Max, mean±std) constraints per particle per d: ", [max_constr, avg_constr, std_constr])
end

"""
    print_info_convergence(final_packing::MonoPacking, isostatic::Bool, time, memory; digs::Int64=2)

Print extra details about the configuration obtained once CALiPPSO converges.

The maximum mismatch in the force balance condition is computed (and printed), and also the
isostaticity of the packing is assessed. Besides, some information about performance
(execution time and allocated memory) is also printed out.

See also [`print_converged`](@ref), [`print_monitor_progress`](@ref).
"""
function print_info_convergence(final_packing::AbstractPacking{d, T}, isostatic::Bool, time, memory; digs::Int64=2) where {d, T<:AbstractFloat}

    isostatic ? (col_p = :green) : (col_p = :magenta)

    N = length(final_packing)
    non_rattlers = get_non_rattlers(final_packing); Nnr = length(non_rattlers)
    rattlers = setdiff(1:N, non_rattlers)
    zs = get_coordination_number(final_packing, only_stable=false)

    force_mismatch = maximum(norm.(total_force(final_packing)))
    minutes = round(time/60, digits=digs)
    memory_r = round(memory, digits=digs)

    if isostatic
        printstyled("Isostaticity achieved: ", isostatic, "\t Non-rattlers = ", Nnr, "\t% of rattlers = ", round(100*(1-Nnr/N), digits=3), "\tN_contacts [in stable particles] = ", 0.5*sum(zs[non_rattlers]),"\t z in rattlers = ", sum(zs[rattlers]), "\n", color=col_p)
        printstyled("Maximum force equilibrium mismatch = ", force_mismatch, "\n", color=col_p)
    else
        iso_gap = (d*(Nnr-1)+1) - 0.5*sum(zs)
        printstyled("Isostaticity achieved: ", isostatic, "  ; isostaticity gap = ", iso_gap, "\t Non-rattlers = ", Nnr, "\t% of rattlers = ", round(100*(1-Nnr/N), digits=3), "\tN_contacts [all particles] = ", 0.5*sum(zs),"\t z in rattlers = ", sum(zs[rattlers]), "\n", color=col_p)
        printstyled("Maximum force equilibrium mismatch = ", force_mismatch, "\n", color=col_p)
    end
    println("Time to finish = ", minutes, " minutes;\t Memory allocated (GB): ",  memory_r)
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
function bandec_trans(a::Array,n::Int64,m1::Int64,m2::Int64,al::Array,indx::Vector{Int64})
#      INTEGER m1,m2,mp,mpl,n,np,indx(n)
#      REAL d,a(np,mp),al(np,mpl),TINY
mm=m1+m2+1
@assert(size(a)  == (mm,n))
@assert(size(al) == (m1,n))
@assert(length(indx) == n)
TINY = one(eltype(a))*1.e-20
#      INTEGER i,j,k,l,mm
#      REAL dum
l=m1
czero = zero(eltype(a))
#dum = czero
for i=1:m1
  for j=m1+2-i:mm
    a[j-l,i]=a[j,i]
  end
  l -= 1
  for j=mm-l:mm
    a[j,i]=czero
  end
end
d = one(eltype(a))
l=m1
for k=1:n
  dum=a[1,k]
  i=k
  if (l < n)
    l += 1
  end
  for j=k+1:l
    if (abs(a[1,j]) > abs(dum))
      dum=a[1,j]
      i=j
    end
  end
  indx[k]=i
  if (dum == czero)
    a[1,k]=TINY
  end
  if (i != k)
    d = -d
    for j=1:mm
      dum=a[j,k]
      a[j,k]=a[j,i]
      a[j,i]=dum
    end
  end
  for i=k+1:l
    dum=a[1,i]/a[1,k]
    al[i-k,k]=dum
    for j=2:mm
      a[j-1,i] = a[j,i] - dum*a[j,k]
    end
    a[mm,i]=czero
  end
end
return d
end

ext: jl


text:
"""
    Cartpole w/ internal friction
"""
struct Friction end
struct Frictionless end

struct Cartpole{T,X} <: Model{T}
    # dimensions
    nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
	nc::Int # contact points

    mc::T     # mass of the cart in kg
    mp::T     # mass of the pole (point mass at the end) in kg
    l::T      # length of the pole in m
    g::T      # gravity m/s^2

    friction::Vector{T} # friction coefficients for slider and arm joints
end

function kinematics(model::Cartpole, q)
    [q[1] + model.l * sin(q[2]); -model.l * cos(q[2])]
end

lagrangian(model::Cartpole, q, q̇) = 0.0

function M_func(model::Cartpole, x)
    H = [model.mc + model.mp model.mp * model.l * cos(x[2]);
		 model.mp * model.l * cos(x[2]) model.mp * model.l^2.0]
    return H
end

function B_func(model::Cartpole, x)
    [1.0; 0.0]
end

function P_func(model::Cartpole, x)
    [1.0 0.0;
     0.0 1.0]
end

function C_func(model::Cartpole, q, q̇)
    C = [0.0 -1.0 * model.mp * q̇[2] * model.l * sin(q[2]);
	 	 0.0 0.0]
    G = [0.0,
		 model.mp * model.g * model.l * sin(q[2])]
    return -C * q̇ + G
end

function dynamics(model::Cartpole{T,Friction}, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2) where T
	# evalutate at midpoint
    qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

    D1L1, D2L1 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm1, vm1)
    D1L2, D2L2 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm2, vm2)

    d = 0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2 # variational integrator (midpoint)

	return d + B_func(model, qm2) * u1[1] + transpose(P_func(model, q2)) * λ1
end

function dynamics(model::Cartpole{T,Frictionless}, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2) where T
	# evalutate at midpoint
    qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

    D1L1, D2L1 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm1, vm1)
    D1L2, D2L2 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm2, vm2)

    d = 0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2 # variational integrator (midpoint)

	return d + B_func(model, qm2) * u1[1]
end

function residual(model::Cartpole{T,Friction}, z, θ, κ) where T
    nq = model.nq
    nu = model.nu
    nc = model.nc

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    μ_slider = θ[2nq + nu .+ (1:1)]
    μ_angle = θ[2nq + nu + 1 .+ (1:1)]
    h = θ[2nq + nu + 1 + 1 .+ (1:1)]

    q2 = z[1:nq]
    ψ = z[nq .+ (1:nc)]
    b = z[nq + nc .+ (1:nc)]
    sψ = z[nq + nc + nc .+ (1:nc)]
    sb = z[nq + nc + nc + nc .+ (1:nc)]

    vT1 = (q2[1] - q1[1]) / h
    vT2 = (q2[2] - q1[2]) / h

    λ1 = b

    [
     dynamics(model, a -> M_func(model, a), (a, b) -> C_func(model, a, b),
        h, q0, q1, u1, zeros(model.nw), λ1, q2);
     sb[1] - vT1[1];
     ψ[1] .- μ_slider[1] * (model.mp + model.mc) * model.g * h[1];
     sb[2] - vT2[1];
     ψ[2] .- μ_angle[1] * (model.mp * model.g * model.l) * h[1];
     cone_product([ψ[1]; b[1]], [sψ[1]; sb[1]]) - [κ[1]; 0.0];
     cone_product([ψ[2]; b[2]], [sψ[2]; sb[2]]) - [κ[1]; 0.0];
    ]
end

function residual(model::Cartpole{T,Frictionless}, z, θ, κ) where T
    nq = model.nq
    nu = model.nu

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]

    q2 = z[1:nq]

    return dynamics(model, a -> M_func(model, a), (a, b) -> C_func(model, a, b),
        h, q0, q1, u1, zeros(model.nw), zeros(0), q2);
end

# models
cartpole_friction = Cartpole{Float64,Friction}(2, 1, 0, 2, 1.0, 0.2, 0.5, 9.81, [0.1; 0.1])
cartpole_frictionless = Cartpole{Float64,Frictionless}(2, 1, 0, 2, 1.0, 0.2, 0.5, 9.81, [0.0; 0.0])


text:
function PhaseAmbiguityModifier(arg0::jint, arg1::jdouble)
    return PhaseAmbiguityModifier((jint, jdouble), arg0, arg1)
end

function get_parameters_drivers(obj::PhaseAmbiguityModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::PhaseAmbiguityModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
mutable struct SurvProblem{Tpar <: AbstractParam, Tdata <: AbstractData, Tprior <: AbstractPrior} <: NUTSGibbs.AbstractProblem
    theta::Tpar
    data::Tdata
    prior::Tprior
end

mutable struct SurvProblem2{Tpar <: AbstractParam, Tdata <: AbstractData, Tprior <: AbstractPrior} <: NUTSGibbs.AbstractProblem
    theta::Tpar
    data::Tdata
    prior::Tprior
end

## class-specific delta
function (problem::SurvProblem{Param, Tdata, Prior})(x) where Tdata <: AbstractData
    @unpack theta, data, prior = problem
    @unpack delta, bl_shape, bl_scale = x

    lp = 0.0
    for g in 1:size(delta,1)
        sel = theta.cind .== g
        lp += loglik_p_surv(data.T[sel], data.E[sel], data.Xe[sel,:], bl_shape[g], bl_scale[g], delta[g,:])
        lp += logpdf(prior.bl_shape, bl_shape[g]) + logpdf(prior.bl_scale, bl_scale[g]) + logpdf(prior.delta_mv, delta[g,:])
    end

    isfinite(lp) || return -Inf
    lp
end

function (problem::SurvProblem{ParamSpline, Tdata, PriorSpline})(x) where Tdata <: AbstractData
    @unpack theta, data, prior = problem
    @unpack delta, w = x
    lp = 0.0
    for g in 1:size(delta,1)
        sel = theta.cind .== g
        lp += loglik_p_surv(data.T[sel], data.E[sel], data.Xe[sel,:], data.haz_basis[sel,:], data.cumhaz_basis[sel,:], w[g,:], delta[g,:])
        lp += sum(logpdf.(prior.w, w[g,:])) + logpdf(prior.delta_mv, delta[g,:])
    end
    isfinite(lp) || return -Inf
    lp
end

## common delta across classes
function (problem::SurvProblem2{Param, Tdata, Prior})(x) where Tdata <: AbstractData
    @unpack theta, data, prior = problem
    @unpack delta, bl_shape, bl_scale = x

    lp = 0.0
    for g in 1:length(bl_scale)
        sel = theta.cind .== g
        lp += loglik_p_surv(data.T[sel], data.E[sel], data.Xe[sel,:], bl_shape[g], bl_scale[g], delta)
        lp += logpdf(prior.bl_shape, bl_shape[g]) + logpdf(prior.bl_scale, bl_scale[g])
    end
    lp += logpdf(prior.delta_mv, delta)

    isfinite(lp) || return -Inf
    lp
end

function (problem::SurvProblem2{ParamSpline, Tdata, PriorSpline})(x) where Tdata <: AbstractData
    @unpack theta, data, prior = problem
    @unpack delta, w = x

    v = var(prior.w)

    G = size(w,1)
    lp = 0.0
    for g in 1:G
        sel = theta.cind .== g
        lp += loglik_p_surv(data.T[sel], data.E[sel], data.Xe[sel,:], data.haz_basis[sel,:], data.cumhaz_basis[sel,:], w[g,:], delta)
        lp -= sum(w[g,:].^2)/(2*v)
        ##lp += sum(logpdf.(prior.w, w[g,:]))
    end
    lp += logpdf(prior.delta_mv, delta)

    isfinite(lp) || return -Inf
    lp
end

function transformation(p::SurvProblem{Param, Tdata, Prior}) where Tdata <: AbstractData
    G, r = size(p.theta.delta)
    as((delta = as(Matrix, G, r), bl_shape = as(Vector, asℝ₊, G), bl_scale = as(Vector, asℝ₊, G)))
end

function transformation(p::SurvProblem{ParamSpline, Tdata, PriorSpline}) where Tdata <: AbstractData
    G, r = size(p.theta.delta)
    k = size(p.theta.w,2)
    as((delta = as(Matrix, G, r), w = as(Matrix, G, k)))
end

function transformation(p::SurvProblem2{Param, Tdata, Prior}) where Tdata <: AbstractData
    G, r = size(p.theta.delta)
    as((delta = as(Vector, r), bl_shape = as(Vector, asℝ₊, G), bl_scale = as(Vector, asℝ₊, G)))
end

function transformation(p::SurvProblem2{ParamSpline, Tdata, PriorSpline}) where Tdata <: AbstractData
    G, r = size(p.theta.delta)
    k = size(p.theta.w,2)
    as((delta = as(Vector, r), w = as(Matrix, G, k)))
end

function surv_model_init(theta, data, prior, warmup)
    problem = SurvProblem(theta, data, prior)
    NUTSGibbs.init(problem, transformation, warmup)
end

function surv_model_init2(theta, data, prior, warmup)
    problem = SurvProblem2(theta, data, prior)
    NUTSGibbs.init(problem, transformation, warmup)
end

function surv_model_step!(x::NUTSGibbsSampler, theta::Param)
    x.problem.theta = theta
    res, ac = NUTSGibbs.step(x, (delta=theta.delta, bl_shape=theta.bl_shape, bl_scale=theta.bl_scale))
    theta.delta = res.delta
    theta.bl_shape = res.bl_shape
    theta.bl_scale = res.bl_scale
    theta, ac
end

function surv_model_step!(x::NUTSGibbsSampler, theta::ParamSpline)
    x.problem.theta = theta
    res, ac = NUTSGibbs.step(x, (delta=theta.delta, w=theta.w))
    theta.delta = res.delta
    theta.w = res.w
    theta, ac
end

function surv_model_step2!(x::NUTSGibbsSampler, theta::Param)
    x.problem.theta = theta
    res, ac = NUTSGibbs.step(x, (delta=theta.delta[1,:], bl_shape=theta.bl_shape, bl_scale=theta.bl_scale))

    for g in 1:size(theta.delta,1)
        theta.delta[g,:] = res.delta
    end

    theta.bl_shape = res.bl_shape
    theta.bl_scale = res.bl_scale
    theta, ac
end

function surv_model_step2!(x::NUTSGibbsSampler, theta::ParamSpline)
    x.problem.theta = theta
    res, ac = NUTSGibbs.step(x, (delta=theta.delta[1,:], w=theta.w))

    for g in 1:size(theta.delta,1)
        theta.delta[g,:] = res.delta
    end

    theta.w = res.w
    theta, ac
end


text:
# BiotSav(PointPath,r;Current=1,MinThreshold = 0.01)
# DCore_PointPath(r₁, r₂; dr = (r₂-r₁)/100)
# MakeEllip(r₁,r₂;Center = [0,0,0],NPts=100)
# VecDist(X::Array) = [ X[2,1]-X[1,1], X[2,2]-X[1,2] , X[2,3]-X[1,3]]

"""
This function takes in the following (required) variables:
    IRad: Inner radius of the D-shaped toroid in meters
    ORad: Outer radius of the D-shaped toroid in meters (α*IRad)
    N: Number of (evenly distributed) turns

    kwargs:
    DownSampleFac: this is the amount the D-core point path will be down-sampled  to make the ``Test points''
        It is more accurate to have a highly sampled point path and less finely sampled interior.
    PlotsOn is a boolean that, if true, plots the wires as they are simulated
    NPtsPath is the number of sections on the point path
    NLayers is the number of interior test point layers
    WireRadius is passed to Mutual_L_TwoLoops and (slightly) changes the way test points are found

"""
function eval_Induct_DToroid(IRad, ORad, N::Int;DownSampleFac=1,PlotOn=false,NPtsPath=100,NLayers=20,WireRadius=0.001)
    ΘperTurn = 360 / N
    if PlotOn
        pygui(true)
        gcf()
    end
    PointPath_SingleLoop =  DCore_PointPath(IRad, ORad;NPts=NPtsPath)
    SliceArea = SymmetricPP2SliceArea(PointPath_SingleLoop)
    LMat = zeros(N, N)
    if PlotOn
        scatter3D(PointPath_SingleLoop[:,1], PointPath_SingleLoop[:,2], PointPath_SingleLoop[:,3])
    end
    NewPP = Rot_PointPath_Y(PointPath_SingleLoop, ΘperTurn)
    if PlotOn
        scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3])
    end
    Mut, Self, SavedBSelfArr = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10,WireRadius=WireRadius,IncludeWireInduct=true,SaveΦ₁=true)
    LMat .+= Self .* OffDiagOnes(N, 1)
    LMat .+= Mut .* OffDiagOnes(N, 2)
    LMat .+= Mut .* OffDiagOnes(N, N)
    for i in 2:Int(ceil((N - 1) / 2))
        NewPP = Rot_PointPath_Y(PointPath_SingleLoop, i * ΘperTurn)
        if PlotOn
            scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3])
        end
        Mut = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP, SavedBSelfArr;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10,WireRadius=WireRadius)
        println(Mut)

        LMat .+= Mut .* OffDiagOnes(N, i + 1)
        LMat .+= Mut .* OffDiagOnes(N, N - (i - 1))

        println(i)
    end
    LMat[N,N] = LMat[1,1]
    return LMat, sum(LMat)
end
"""
This function takes in the following (required) variables:
    IRad: Inner radius of the circular-core toroid in meters
    ORad: Outer radius of the circular-core toroid in meters (α*IRad)
    N: Number of (evenly distributed) turns

    kwargs:
    DownSampleFac: this is the amount the D-core point path will be down-sampled  to make the ``Test points''
        It is more accurate to have a highly sampled point path and less finely sampled interior.
    PlotsOn is a boolean that, if true, plots the wires as they are simulated
    NPtsPath is the number of sections on the point path
    NLayers is the number of interior test point layers
    WireRadius is passed to Mutual_L_TwoLoops and (slightly) changes the way test points are found

"""
function eval_Induct_CircToroid(IRad,ORad,N::Int;DownSampleFac=1,
    PlotOn=false,
    NPtsPath=100,
    NLayers=20,WireRadius=nothing)
    ΘperTurn = 360 / N
    if PlotOn
        pygui(true)
        gcf()
    end
    CoreRad = (ORad - IRad) / 2
    CoreOffset = (ORad + IRad) / 2
    PointPath_SingleLoop =  MakeEllip(CoreRad, CoreRad;Center=[CoreOffset,0,0], NPts=NPtsPath)
    SliceArea = SymmetricPP2SliceArea(PointPath_SingleLoop)
    LMat = zeros(N, N)
    if PlotOn
        scatter3D(PointPath_SingleLoop[:,1], PointPath_SingleLoop[:,2], PointPath_SingleLoop[:,3])
    end

    for i in 1:(N - 1)
        NewPP = Rot_PointPath_Y(PointPath_SingleLoop, i * ΘperTurn)
        if PlotOn
            scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3])
        end
        Mut, Self = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10,WireRadius=WireRadius,IncludeWireInduct=true)
        println(Mut)
        LMat[i,i] = Self
        LMat .+= Mut .* OffDiagOnes(N, i + 1)

        println(i)
    end
    LMat[N,N] = LMat[1,1]
    return LMat, sum(LMat)
end

"""
This function takes in the following (required) variables:
    IRad: Inner radius of the D-shaped toroid in meters
    ORad: Outer radius of the D-shaped toroid in meters (α*IRad)
    N₁: Number of (evenly distributed) turns on the primary
    N₂: Number of (evenly distributed) turns on the secondary

    kwargs:
    DownSampleFac: this is the amount the D-core point path will be down-sampled  to make the ``Test points''
        It is more accurate to have a highly sampled point path and less finely sampled interior.
    PlotsOn is a boolean that, if true, plots the wires as they are simulated
    NPtsPath is the number of sections on the point path
    NLayers is the number of interior test point layers
    WireRadius is passed to Mutual_L_TwoLoops and (slightly) changes the way test points are found
    SelfInductScaleFactor is, as the name implies, a small ``fudge factor'' if you want to artificially scale up the self-inductance. This may be useful if there is a known offset to account for.

"""
function eval_Induct_DToroidTransformer(IRad, ORad, N₁::Int, N₂::Int;DownSampleFac=1,PlotOn=false,NPtsPath=100,NLayers=20,SelfInductScaleFactor=1,UseAnalyticalSelfL=false)
    if PlotOn
        pygui(true)
        gcf()
    end

    N = N₁ + N₂

    N₁List = 0:(N / N₁):(N - N / N₁)
    N₂List = 0.0001:(N / N₂):(N - N / N₂ + 0.0001)

    NIndexingArray = hcat(vcat(N₁List, N₂List), vcat(ones(length(N₁List), 1), zeros(length(N₂List), 1)))
    NIndexingArray = sortslices(NIndexingArray;dims=1) # This creates an array that will allow for determining which turn is associated with primary or secondary
    println(NIndexingArray)
    ΘperTurn = 360 / N

    PointPath_SingleLoop =  DCore_PointPath(IRad, ORad;NPts=NPtsPath)
    SliceArea = SymmetricPP2SliceArea(PointPath_SingleLoop)
    LMat = zeros(N, N)
    if PlotOn
        scatter3D(PointPath_SingleLoop[:,1], PointPath_SingleLoop[:,2], PointPath_SingleLoop[:,3];color="b")
    end
    NewPP = Rot_PointPath_Y(PointPath_SingleLoop, ΘperTurn)
    Mut, Self, SavedBSelfArr = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10,IncludeWireInduct=true,SaveΦ₁=true)
    LMat .+= Self .* OffDiagOnes(N, 1)
    LMat .+= Mut .* OffDiagOnes(N, 2)
    LMat .+= Mut .* OffDiagOnes(N, N)
    for i in 1:Int(ceil((N - 1) / 2))
        isPrimary = (NIndexingArray[i + 1,2] == 1)
        NewPP = Rot_PointPath_Y(PointPath_SingleLoop, i * ΘperTurn)
        if PlotOn
            if isPrimary
                scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3];color="b",s=2)
            else
                scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3];color="r",s=2)
            end
        end
        Mut = Mutual_L_TwoLoops(PointPath_SingleLoop, NewPP, SavedBSelfArr;DownSampleFac=DownSampleFac,MeasLayers=NLayers,MinThreshold=1e-10)
        println(Mut)

        LMat .+= Mut .* OffDiagOnes(N, i + 1)
        LMat .+= Mut .* OffDiagOnes(N, N - (i - 1))

        println(i)
    end
    LMat[N,N] = LMat[1,1]
    Primary_Self, Secondary_Self, MutualInduct = LMatMutualIndexing(N, NIndexingArray[:,2])
    L₁₁ = LMat[Primary_Self] .* SelfInductScaleFactor
    L₂₂ = LMat[Secondary_Self] .* SelfInductScaleFactor
    LMut = LMat[MutualInduct]
    if UseAnalyticalSelfL
        println("Using analyical formula for self inductance. Make sure ORad/IRad is an integer")
        L₁₁ = DCore_DetermineIdealInduct(IRad * 2, N₁, Int(round(ORad / IRad)))
        L₂₂ = DCore_DetermineIdealInduct(IRad * 2, N₂, Int(round(ORad / IRad)))
    end
    Lμ = sum(LMut) * N₁ / N₂ / 2 # THe division by 2 is due to the double-counting of M in a symmetric matrix
    LLeak1 = sum(L₁₁) - Lμ
    LLeak2 = sum(L₂₂) - Lμ * (N₂ / N₁)^2

    K = (sum(LMut) / 2) / √(sum(L₁₁) * sum(L₂₂))
    println("L₁₁ = $(round(sum(L₁₁);sigdigits=3))")
    println("L₂₂ = $(round(sum(L₂₂);sigdigits=3))")
    println("Lₘ = $(round(sum(LMut) / 2;sigdigits=3))")
    println("K = $(round(sum(K);sigdigits=3))")
    println("Lμ = $(round(Lμ;sigdigits=3))")
    println("LLeak1 = $(round(LLeak1;sigdigits=3))")
    println("LLeak2 = $(round(LLeak2;sigdigits=3))")
    return LMat, NIndexingArray, L₁₁, L₂₂, LMut, Lμ, LLeak1, LLeak2, K
end

function eval_Induct_DToroid_SplitWindings(IRad, ORad, N₁Half::Int)

    figure()
    LMat, NIndexingArray, L₁₁, L₂₂, LMut, Lμ, LLeak1, LLeak2, K = CircModDesign.eval_Induct_DToroidTransformer(IRad,ORad,N₁Half,N₁Half;
    DownSampleFac=10,
    PlotOn=true,NLayers=45)
    Primary_Self, Secondary_Self, MutualInduct = CircModDesign.LMatMutualIndexing(N₁Half * 2, NIndexingArray[:,2])
    xlabel("X, meters")
    ylabel("Y, meters")
    zlabel("Z, meters")
    xlim([-1.1 * ORad, 1.1 * ORad])
    zlim([-1.1 * ORad, 1.1 * ORad])
    ylim([-1.1 * ORad, 1.1 * ORad])

    BundleScaling = N₁Half
    figure()
    CircModDesign.ScatterPlotToroid(IRad,ORad,Int(N₁Half),Int(N₁Half),BundleScaling)
    Primary_Self, Secondary_Self, MutualInduct = CircModDesign.LMatMutualIndexing(N₁Half * 2, CircModDesign.makeIndexinVec(NIndexingArray[:,2], BundleScaling))
    xlabel("X, meters")
    ylabel("Y, meters")
    zlabel("Z, meters")
    xlim([-1.1 * ORad, 1.1 * ORad])
    zlim([-1.1 * ORad, 1.1 * ORad])
    ylim([-1.1 * ORad, 1.1 * ORad])
    yticks([])
    ylabel("")
    L₁₁ = LMat[Primary_Self]
    L₂₂ = LMat[Secondary_Self]
    LMut = LMat[MutualInduct]
    Lμ = sum(LMut) / 2 # THe division by 2 is due to the double-counting of M in a symmetric matrix
    LLeak1 = sum(L₁₁) - Lμ
    LLeak2 = sum(L₂₂) - Lμ
    K = (sum(LMut) / 2) / √(sum(L₁₁) * sum(L₂₂))
    figure()
    subplot(131)
    imshow(Primary_Self)
    title("Self-Inductance Indices - Primary")
    subplot(132)
    imshow(Secondary_Self)
    title("Self-Inductance Indices - Secondary")
    subplot(133)
    imshow(MutualInduct)
    title("Mutual Inductance Indices")
    return LMat, NIndexingArray, L₁₁, L₂₂, LMut, Lμ, LLeak1, LLeak2, K
end

"""
This builds a rotation matrix, to rotate around the z axis
Radians input.
"""
function BuildRotMatrix_deg_Z(Θ)
    Θ = Θ * π / 180
   Mat =  [cos(Θ) -1 * sin(Θ) 0
           sin(Θ)  cos(Θ)   0
           0       0        1]
end
"""
This builds a rotation matrix, to rotate around the Y axis
Radians input.
"""
function BuildRotMatrix_deg_Y(Θ)
    Θ = Θ * π / 180
   Mat =  [cos(Θ)  0   -1 * sin(Θ)
           0       1   0
           sin(Θ)       0   cos(Θ)]
end

function Rot_PointPath_Y(Point_Path, Θ)
    Point_PathTrans = transpose(Point_Path)
    return transpose(BuildRotMatrix_deg_Y(Θ) * Point_PathTrans)
end

function Rot_PointPath_Z(Point_Path, Θ)
    Point_PathTrans = transpose(Point_Path)
    return transpose(BuildRotMatrix_deg_Z(Θ) * Point_PathTrans)
end

"""
This function takes in two point paths and simulates the self inductance of the first, and mutual inductance between them.

    kwargs:
    DownSampleFac: this is the amount the point path will be down-sampled  to make the ``Test points''
        It is more accurate to have a highly sampled point path and less finely sampled interior.

    MeasLayers is the number of interior test point layers
    WireRadius (slightly) changes the way test points are found. It mostly just adds a set of test points near the wire if the distance is large between the outer layer of test points and wire edge
    IncludeWireInduct is a boolean. If true, adds the wire's internal inductance to the self-inductance of the first
    MinThreshold should not be modified to be larger than the wire radius. It essentially makes the bior savart function ignore test points under that treshold, which can be useful to prevent divisions by zero in edge cases.

"""
function Mutual_L_TwoLoops(PP₁,PP₂;DownSampleFac=1,
    MeasLayers=55,# Number of concentric layers to calc field
    MinThreshold=0.000001,
    WireRadius=0.001,
    IncludeWireInduct=false,
    SaveΦ₁=false
    )

    TestPoints, Weights = PP_2_TestPoints_drdΘ(PP₁, MeasLayers, WireRadius;DownSampleFac=DownSampleFac)
    PtsPerLayer = length(PP₁[:,1])
    NPtsPerSlice = MeasLayers * PtsPerLayer
    SliceArea = sum(Weights)
        Weights = Weights / sum(Weights) * NPtsPerSlice  # normalizing
        if SaveΦ₁
            B_Self = zeros(length(TestPoints[:,1]), 3)
            Φ₁₁ = 0.
        else
            Φ₁₁ = 0.
        end

    Φ₂₁ = 0.
    for i in 1:length(TestPoints[:,1])

        if SaveΦ₁
            B_Self[i,:] =  BiotSav(PP₁, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice
            Φ₁₁ += √(sum(B_Self[i,:].^2))
            Φ₂₁ += sum(B_Self[i,:]  .* BiotSav(PP₂, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice) /    √(sum(B_Self[i,:].^2))
        else
            B_Self =  BiotSav(PP₁, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice
            Φ₁₁ += √(sum(B_Self.^2))

            Φ₂₁ += sum(B_Self .* BiotSav(PP₂, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice) /    √(sum(B_Self.^2))
        end

    end

    CumulativeDist = 0
    if IncludeWireInduct
        for i in 1:(length(PP₁[:,1]) - 1)
            CumulativeDist += sqrt(sum((PP₁[i,:] .- PP₁[i + 1,:]).^2))
        end
        println(CumulativeDist)
        WireInductance = internalInductance(CumulativeDist)
        Φ₁₁ += WireInductance
    end

    return Φ₂₁, Φ₁₁, B_Self

end

function Mutual_L_TwoLoops(PP₁,PP₂,PriorBSelf;DownSampleFac=1,
    MeasLayers=55,# Number of concentric layers to calc field
    MinThreshold=0.000001,
    WireRadius=0.001,

    )

    TestPoints, Weights = PP_2_TestPoints_drdΘ(PP₁, MeasLayers, WireRadius;DownSampleFac=DownSampleFac)
    PtsPerLayer = length(PP₁[:,1])
    NPtsPerSlice = MeasLayers * PtsPerLayer
    SliceArea = sum(Weights)
        Weights = Weights / sum(Weights) * NPtsPerSlice  # normalizing

    Φ₂₁ = 0.

    PP₂ = vcat(PP₂,PP₂[1,:]')
    dL = ([VecDist(PP₂[I-1:I,:]) for I in 2:length(PP₂[:,1])])
    PP₂_RS = ([PP₂[I,:] for I in 1:length(PP₂[:,1])])
    L = length(PP₂_RS[:,1])

    #BiotSav(PointPath,dL,r,L) is the fast version
    for i in 1:length(TestPoints[:,1])
        Φ₂₁ += sum(PriorBSelf[i,:] .* BiotSav(PP₂_RS,dL, TestPoints[i,:],L) .* Weights[i] .* SliceArea ./ NPtsPerSlice) /    √(sum(PriorBSelf[i,:].^2))

    end

    return Φ₂₁

end

function SymmetricPP2SliceArea(PP)
Total = 0
    for i in 1:(Int(round(length(PP[:,1]) / 2)))
        dX = PP[i + 1,1] - PP[i,1]
        MeanY = (PP[i + 1,2] + PP[i,2]) / 2
        Total += dX * MeanY
    end

    return Total * 2 # integrating over half, then multiplying by 2 to account for symmetry
end

function OffDiagOnes(FullSize, Dist)
    M = zeros(FullSize, FullSize)
    j = 0
    for i in Dist:(FullSize)
        j += 1
        M[i,j] = 1
        M[j,i] = 1
    end
    return M
end

function LMatMutualIndexing(NTotal, IndexingVector_Primary)

    Primary_Self = zeros(NTotal, NTotal)
    Secondary_Self = zeros(NTotal, NTotal)
    MutualInduct = zeros(NTotal, NTotal)

    for i in 1:NTotal, j in 1:NTotal
        if ((IndexingVector_Primary[i] == 1) & (IndexingVector_Primary[j] == 1))
            Primary_Self[i,j] = 1
        elseif ((IndexingVector_Primary[i] == 0) & (IndexingVector_Primary[j] == 0))
            Secondary_Self[i,j] = 1
        end
    end
    MutualInduct = (Primary_Self .+ Secondary_Self) .== 0
    Primary_Self = Primary_Self .== 1
    Secondary_Self = Secondary_Self .== 1
    return Primary_Self, Secondary_Self, MutualInduct

end

function Mutual_L_CircTwoLoops_Separation(D,SeparationVec  ;DownSampleFac=1,
    MeasLayers=55,# Number of concentric layers to calc field
    MinThreshold=1e-15, NPts=100,PlotOn=false
    )

    PP₁ = MakeEllip(D / 2, D / 2;NPts=NPts)
    if PlotOn
        scatter3D(PP₁[:,1], PP₁[:,2], PP₁[:,3];color="r")
    end
    SliceArea = π / 4 * D^2

    MutualLMat = zeros(length(SeparationVec))
    Counter = 0
    for j in SeparationVec
        Counter += 1
        PP₂ = MakeEllip(D / 2, D / 2;Center=[0,0,j])
        if PlotOn
            scatter3D(PP₂[:,1], PP₂[:,2], PP₂[:,3];color="b")
        end
        TestPoints, Weights = PP_2_TestPoints_drdΘ(PP₁, MeasLayers;DownSampleFac=DownSampleFac)
        PtsPerLayer = length(PP₁[:,1])
        NPtsPerSlice = MeasLayers * PtsPerLayer
        Weights = Weights / sum(Weights) * NPtsPerSlice  # normalizing

        Φ₁₁ = 0.
        Φ₂₁ = 0.
        for i in 1:length(TestPoints[:,1])
            B_Self =  BiotSav(PP₁, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice
            Φ₁₁ += √(sum(B_Self.^2))
            Φ₂₁ += sum(B_Self .* BiotSav(PP₂, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice) /    √(sum(B_Self.^2))

        end
        MutualLMat[Counter] = Φ₂₁
    end

    return MutualLMat

end

function Self_L_CircLoop(D  ; DownSampleFac=1,
    MeasLayers=55,# Number of concentric layers to calc field
    MinThreshold=1e-15, NPts=100,PlotOn=false,WireRadius=0.0005,IncludeWireInduct=true
    )

    PP₁ = MakeEllip(D / 2, D / 2;NPts=NPts)
    if PlotOn
        scatter3D(PP₁[:,1], PP₁[:,2], PP₁[:,3];color="r")
    end
    SliceArea = π / 4 * D^2

        PP₂ = MakeEllip(D / 2, D / 2)
        if PlotOn
            scatter3D(PP₂[:,1], PP₂[:,2], PP₂[:,3];color="b")
        end
        DistBetweenLayers = D / 2 / MeasLayers
        if DistBetweenLayers > WireRadius
            TestPoints, Weights = PP_2_TestPoints_drdΘ(PP₁, MeasLayers, WireRadius;DownSampleFac=DownSampleFac)
        else
            TestPoints, Weights = PP_2_TestPoints_drdΘ(PP₁, MeasLayers, WireRadius;DownSampleFac=DownSampleFac)
        end
        PtsPerLayer = length(PP₁[:,1])
        NPtsPerSlice = MeasLayers * PtsPerLayer
        Weights = Weights / sum(Weights) * NPtsPerSlice  # normalizing

        Φ₁₁ = 0.
        Φ₂₁ = 0.
        for i in 1:length(TestPoints[:,1])
            B_Self =  BiotSav(PP₁, TestPoints[i,:];MinThreshold) .* Weights[i] .* SliceArea ./ NPtsPerSlice
            Φ₁₁ += √(sum(B_Self.^2))
            # Φ₂₁ += sum(B_Self .* BiotSav(PP₂,TestPoints[i,:];MinThreshold) .* Weights[i].*SliceArea ./NPtsPerSlice) /    √(sum(B_Self.^2))

        end

        if IncludeWireInduct
            if DistBetweenLayers < WireRadius
                WireInductance = 0
            else
                WireInductance = internalInductance(π * D)
            end
            println("Adding wire internal inductance")
            println("Distance between layers = : $(round(DistBetweenLayers;sigdigits=3))")
            Φ₁₁ += WireInductance
        end

    return Φ₁₁

end

function internalInductance(L)
    μ₀ = 4 * pi * 1e-7
    return μ₀ * L / (8 * π)
end

function makeIndexinVec(BasisVec, Scale)
    NewVec = []
    N = length(BasisVec)
    for i in 1:length(BasisVec), j in 1:Scale

        push!(NewVec, BasisVec[i])

    end
    NewVec = NewVec[1:N]
    return NewVec
end

function ScatterPlotToroid(IRad, ORad, N₁, N₂, BundleScale;NPtsPath=100,NLayers=20)

        pygui(true)
        gcf()

    N = N₁ + N₂

    N₁List = 0:(N / N₁):(N - N / N₁)
    N₂List = 0.0001:(N / N₂):(N - N / N₂ + 0.0001)

    NIndexingArray = hcat(vcat(N₁List, N₂List), vcat(ones(length(N₁List), 1), zeros(length(N₂List), 1)))
    NIndexingArray = sortslices(NIndexingArray;dims=1) # This creates an array that will allow for determining which turn is associated with primary or secondary
    PrimaryIndex = makeIndexinVec(NIndexingArray[:,2], BundleScale)
    ΘperTurn = 360 / N

    PointPath_SingleLoop =  DCore_PointPath(IRad, ORad;NPts=NPtsPath)

        scatter3D(PointPath_SingleLoop[:,1], PointPath_SingleLoop[:,2], PointPath_SingleLoop[:,3];color="b",s=3)

    for i in 1:(N - 1)
        isPrimary = (PrimaryIndex[i + 1] == 1)
        NewPP = Rot_PointPath_Y(PointPath_SingleLoop, i * ΘperTurn)

            if isPrimary
                scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3];color="b",s=3)
            else
                scatter3D(NewPP[:,1], NewPP[:,2], NewPP[:,3];color="r",s=3)
            end

    end

    return
end

"""
This function takes in a point path and produces test points and weights for each.

    Inputs are a point path and the number of test point layers. A ``layer'' is just a concentric set of points to the point path with a smaller radius. Using this sampling approach gives an uneven distribution of points, but ensures they are all on the interior.

"""
function PP_2_TestPoints_drdΘ(PP, Layers=10, WireRadius=0.001; DownSampleFac::Int=1)

    PP = PP[1:DownSampleFac:end,:]
    PP_Mean = [sum(PP[:,i]) / length(PP[:,i]) for i in 1:3]
    PP_Cent = hcat([PP[:,1] .- PP_Mean[1] , PP[:,2] .- PP_Mean[2], PP[:,3] .- PP_Mean[3]]...)

    TestPoints_Cent = vcat([PP_Cent .* √(NLayer / (Layers + 1))  for NLayer in 1:Layers]...)
    TestPoints_Cent = reverse(TestPoints_Cent;dims=1)
    MaxDist = maximum(sqrt.(sum(PP_Cent.^2;dims=2)))
    TestPoints_Cent = vcat(reverse(PP_Cent .* ((MaxDist - WireRadius) / MaxDist);dims=1), TestPoints_Cent)
    Layers = Layers + 1 # Because adding an additional layer near wire

    TestPoints_Polar = Cart2Polar(TestPoints_Cent)
    PP_Cent_Polar = Cart2Polar(PP_Cent)
    dΘ = zeros(length(TestPoints_Polar[:,1]), 1)
    dr = similar(dΘ)

    # dΘ = abs.(TestPoints_Polar[:,2] .- circshift(TestPoints_Polar[:,2],1))
    # dΘ = minimum(hcat(dΘ,  abs.(TestPoints_Polar[:,2] .- circshift(TestPoints_Polar[:,2],-1)));dims=2)
    L = length(PP_Cent[:,1])

    Counter = 0
    for j in 1:Layers

        for k in 1:L
            Counter += 1

            if k == 1
                dΘ[Counter] = AbsAngBetw(TestPoints_Polar[Counter,2], TestPoints_Polar[Counter + 1,2])

            elseif k == L
                dΘ[Counter] = AbsAngBetw(TestPoints_Polar[Counter,2], TestPoints_Polar[Counter - 1,2])
            else
                dΘ[Counter] = AbsAngBetw(TestPoints_Polar[Counter - 1,2], TestPoints_Polar[Counter + 1,2]) / 2
            end

        end
        StartInd = ((j - 1) * L + 1)
        StopInd = j * L

        StartInd2 = (j * L + 1)
        StopInd2 = (j + 1) * L

        if j !== Layers
            dr[StartInd:StopInd] = TestPoints_Polar[StartInd:StopInd,1] .- TestPoints_Polar[StartInd2:StopInd2,1]
        else
            dr[StartInd:StopInd] = TestPoints_Polar[StartInd:StopInd,1]
        end
    end

    TestPoints = hcat([TestPoints_Cent[:,1] .+ PP_Mean[1] , TestPoints_Cent[:,2] .+ PP_Mean[2], TestPoints_Cent[:,3] .+ PP_Mean[3]]...)

    Weights = abs.(TestPoints_Polar[:,1] .* dr[:] .* dΘ[:])
    r = TestPoints_Polar[:,1]
    return TestPoints, Weights, dr, dΘ, r
end

function AbsAngBetw(θ₁, θ₂)
    if θ₁ == θ₂
        return 0.
    elseif (abs(θ₁ - θ₂) > π)
        return (2 * π - mod2pi(abs(θ₁ - θ₂)))
    else
        return abs(θ₁ - θ₂)
    end
end

function PP_2_TestPoints(PP, Layers=10;  DownSampleFac::Int=1,WireRadius=nothing)
    PP = PP[1:DownSampleFac:end,:]
    PP_Mean = [sum(PP[:,i]) / length(PP[:,i]) for i in 1:3]
    PP_Cent = hcat([PP[:,1] .- PP_Mean[1] , PP[:,2] .- PP_Mean[2], PP[:,3] .- PP_Mean[3]]...)

    TestPoints_Cent = vcat([PP_Cent .* (NLayer / (Layers + 1))  for NLayer in 1:Layers]...)

    if ~(WireRadius === nothing)
        println("Adding inner test layer")
        MaxDist = maximum(sqrt.(sum(PP_Cent.^2;dims=2)))
        TestPoints_Cent = vcat(TestPoints_Cent, PP_Cent .* ((MaxDist - WireRadius) / MaxDist))

    end
    TestPoints = hcat([TestPoints_Cent[:,1] .+ PP_Mean[1] , TestPoints_Cent[:,2] .+ PP_Mean[2], TestPoints_Cent[:,3] .+ PP_Mean[3]]...)

    Weights = [√(TestPoints_Cent[i,1]^2 + TestPoints_Cent[i,2]^2 + TestPoints_Cent[i,3]^2) for i in 1:length(TestPoints_Cent[:,1]) ]

    return TestPoints, Weights
end

function Cart2Polar(Arr)
    PolarArr = similar(Arr)
    for i in 1:length(Arr[:,1])
        PolarArr[i,1] = √(Arr[i,1]^2 + Arr[i,2]^2)
        PolarArr[i,2] = atan(Arr[i,2], Arr[i,1])
        PolarArr[i,3] = Arr[i,3]
    end

    return PolarArr
end


text:
# RCA -----------------------------------------------------
# [ ] Test result against results from groupby approach

function get_rca(M::Array{Float64, 2}; binary::Bool=true)

	# Matrix dimensions
	C = size(M, 1)
	P = size(M, 2)

	# Total country export = row sum
	country_exports = M * ones(P, 1)

	# World product export = col sum
	world_product_exports = M' * ones(C, 1)

	# World total export = mat sum
	total_world_export = sum(M)

	# Calculate:
	RCA = (M ./ country_exports) ./ (world_product_exports ./ total_world_export)'

	# Return:
	if binary == true
		RCA = [i >= 1 ? 1 : 0 for i in RCA]
	end

	return RCA
end

# RpcA ----------------------------------------------------
# TODO: Test

function get_rpca(M::Array{Float64, 2}, POP::Array{Float64, 1}; binary::Bool=true)

	# Check inputs:
	size(M, 1) == size(POP, 1) || return "M and POP must have the same number of rows."

	C = size(M, 1)

	# World export in p
	world_product_exports = M' * ones(C, 1)

	# World pop
	world_pop = sum(POP)

	RPCA = (M ./ POP) ./ (world_product_exports ./ world_pop)'

	if binary == true
		RPCA = [i >= 1 ? 1 : 0 for i in RPCA]
	end

	return RPCA
end

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


text:
#####
##### Counter gradient flux model proposed by Large et al (1994)
#####

Base.@kwdef struct LMDCounterGradientFlux{T} <: AbstractParameters
    CNL :: T = 6.33 # Mass flux proportionality constant
end

const CGModel = Model{K, <:LMDCounterGradientFlux} where K

update_nonlocal_flux!(model::CGModel) = nothing

instantiate_plume(::LMDCounterGradientFlux, grid) = (T=nothing, S=nothing, W²=nothing)

mass_flux(m::CGModel, i) = 0

∂z_explicit_nonlocal_flux_T(m::Model{K, <:LMDCounterGradientFlux}, i) where K =
    KPP.∂NL∂z(m.nonlocalflux.CNL, m.state.Qθ, d(m, i+1), d(m, i), Δf(m.grid, i), m)

∂z_explicit_nonlocal_flux_S(m::Model{K, <:LMDCounterGradientFlux}, i) where K =
    KPP.∂NL∂z(m.nonlocalflux.CNL, m.state.Qs, d(m, i+1), d(m, i), Δf(m.grid, i), m)

#####
##### Diagnostic plume model
#####

Base.@kwdef struct DiagnosticPlumeModel{T} <: AbstractParameters
     Ca :: T = 0.1
    Cbw :: T = 2.86
     Ce :: T = 0.4
    Cew :: T = 0.572
     Cα :: T = 1.0
    Cστ :: T = 2.2
    Cσb :: T = 1.32
end

instantiate_plume(::DiagnosticPlumeModel, grid) =
    (T=CellField(grid), S=CellField(grid), W²=FaceField(grid))

#####
##### Empirical standard deviation
#####

@inline function w_standard_dev(ωb, ωτ, Cστ, Cσb, d::T) where T
    if 0 < d < 1
        return (Cστ * ωτ^3 + Cσb * ωb^3 * d)^(1/3) * (1 - d)^(1/2)
    else
        return zero(T)
    end
end

@inline w_standard_dev(m, i) = @inbounds w_standard_dev(ωb(m), ωτ(m), m.nonlocalflux.Cστ,
                                                        m.nonlocalflux.Cσb, -m.grid.zc[i] / m.state.h)

#####
##### Entrainment
#####

@inline entrainment(z, Δz, Ce, h) = Ce * (1 / (Δz - z) + 1 / (Δz + z + h))

@inline entrainment(z, model) =
    @inbounds entrainment(z, Δc(model.grid, model.grid.N), model.nonlocalflux.Ce, model.state.h)

#####
##### Plume boundary conditions
#####

function set_tracer_plume_bc!(ϕ̆, ϕ, Qϕ, Cα, model)
    n = ϕ̆.grid.N

    # Surface layer model: √w² Δϕ̆ = - C Qϕ, where Δϕ̆ is plume excess.
    @inbounds ϕ̆[n] = ϕ[n] - Cα * Qϕ / w_standard_dev(model, n)

    return nothing
end

function set_vertical_momentum_plume_bc!(W²)
    @inbounds W²[W².grid.N+1] = 0
    return nothing
end

#####
##### Plume equations
#####

#####
##### Buoyancy
#####

@inline plume_buoyancy_excess(i, grid, T̆, S̆, T, S, α, β, g) =
    @inbounds g * (α*(T̆[i] - T[i]) - β*(S̆[i] - S[i]))

@inline plume_buoyancy_excess(i, grid, model) =
    plume_buoyancy_excess(i, grid, model.state.plume.T, model.state.plume.S,
                                   model.solution.T, model.solution.S,
                                   model.constants.α, model.constants.β, model.constants.g)

function update_nonlocal_flux!(model::Model{K, <:DiagnosticPlumeModel}) where K

    set_tracer_plume_bc!(model.state.plume.T, model.solution.T,
                         model.state.Qθ, model.nonlocalflux.Cα, model)

    set_tracer_plume_bc!(model.state.plume.S, model.solution.S,
                         model.state.Qs, model.nonlocalflux.Cα, model)

    set_vertical_momentum_plume_bc!(model.state.plume.W²)

    integrate_plume_equations!(model.state.plume.T, model.state.plume.S, model.state.plume.W²,
                               model.solution.T, model.solution.S, model.grid, model)

    return nothing
end

function integrate_plume_equations!(T̆, S̆, W̆², T, S, grid, model)

    n = grid.N

    # Vertical momentum at the nᵗʰ cell interface, approximating excess buoyancy at
    # interface with excess at top cell center:
    ΔB̆ᵢ₊₁ = plume_buoyancy_excess(n, grid, model)
    @inbounds W̆²[n] = -model.nonlocalflux.Cbw * Δf(grid, n) * ΔB̆ᵢ₊₁

    # Integrate from surface cell `N-1` downwards
    for i in n-1 : -1 : 1
        if W̆²[i+1] <= 0 # plume is stopped
            @inbounds W̆²[i] = 0
            @inbounds T̆[i] = T[i]
            @inbounds S̆[i] = S[i]
        else # plume still lives

            # Integrate temperature and salinity
            @inbounds T̆[i] = T̆[i+1] + Δc(grid, i+1) * entrainment(grid.zc[i+1], model) * (T̆[i+1] - T[i+1])
            @inbounds S̆[i] = S̆[i+1] + Δc(grid, i+1) * entrainment(grid.zc[i+1], model) * (S̆[i+1] - S[i+1])

            # Integrate vertical momentum
            ΔB̆ᵢ₊₁ = onface(i+1, grid, plume_buoyancy_excess, model)

            @inbounds W̆²[i] = W̆²[i+1] - Δf(grid, i+1) * (
                                  model.nonlocalflux.Cbw * ΔB̆ᵢ₊₁
                                - model.nonlocalflux.Cew * entrainment(grid.zf[i+1], model) * W̆²[i+1])
        end
    end

    return nothing
end

#####
##### Mass flux contribution to tracer budgets
#####

maxzero(ϕ::T) where T = max(zero(T), ϕ)

mass_flux(m::Model{K, <:DiagnosticPlumeModel}, i) where K =
    @inbounds -m.nonlocalflux.Ca * sqrt(maxzero(m.state.plume.W²[i]))

@inline M_Φ(i, grid, Φ, model) = @inbounds mass_flux(model, i) * Φ[i]

# Use upwards-biased difference to effect upwind differencing for a downward-travelling plume:
∂z_explicit_nonlocal_flux_T(m::Model{K, <:DiagnosticPlumeModel}, i) where K =
    @inbounds ∂z⁺(i, m.grid, M_Φ, m.state.plume.T, m)

∂z_explicit_nonlocal_flux_S(m::Model{K, <:DiagnosticPlumeModel}, i) where K =
    @inbounds ∂z⁺(i, m.grid, M_Φ, m.state.plume.S, m)


text:
"""
    Fermi.HartreeFock

Module for running Hartree--Fock computations in Fermi.

# Methods

    > Fermi.HartreeFock.RHF
    > Fermi.HartreeFock.UHF
"""
module HartreeFock
# Import Fermi basics
using LinearAlgebra: hermitian, vcat
using Fermi
using Fermi.Options
using Fermi.Integrals
using Fermi.Orbitals

function hf_header()
    output(repeat("=",80))
    output("|{:33}{:^12}{:33}|", "", "Hartree-Fock", "")
    output("|{:34}{:^9}{:34}|", "", "Module  by","")
    output("|{:25}{:^28}{:25}|", "", "G.J.R Aroeira and M.M. Davis", "")
    output(repeat("=",80))
end

function uhf_header()
    output(repeat("=",80))
    output("|{:33}{:^12}{:33}|", "", "Hartree-Fock", "")
    output("|{:34}{:^9}{:34}|", "", "Module  by","")
    output("|{:17}{:^44}{:17}|", "", "G.J.R Aroeira, M.M. Davis, and S.M. Goodlett", "")
    output(repeat("=",80))
end

"""
    Fermi.HartreeFock.AbstractHFWavefunction

Abstract type common to all Hartree-Fock wave functions.

Struct tree

**AbstractHFWavefunction** <: AbstractWavefunction
"""
abstract type AbstractHFWavefunction <: Fermi.AbstractWavefunction end

"""
    Fermi.HartreeFock.get_scf_alg()

Returns a singleton type corresponding to a RHF implementation based on the options.
"""

function get_scf_alg()
    ref = Options.get("reference")
    if ref == "rhf"
        implemented = [RHFa()]
    elseif ref == "uhf"
        implemented = [UHFa()]
    else
        throw(FermiException("$ref is not a supported reference energy"))
    end

    N = Options.get("scf_alg")
    try
        return implemented[N]
    catch BoundsError
        throw(FermiException("implementation number $N not available for $ref."))
    end
end

# Different Hartree-Fock methods are included here:

# Restricted Hartree--Fock
include("RHF/RHF.jl")

# UHF
include("UHF/UHF.jl")

end #module


text:
"Perform Matrix Operations Such as a convolutional layer followed by a pooling layer applied on the filters"
function MatrixWise(layer::LayerGenerator)
    return LayerGenerator(
        function(pb, sett, input_shape)
                filter_count = input_shape[end]
                design = layer(pb, sett, (input_shape[1:(end - 1)]))

                LayerDesign(Sig(:MatrixWise), input_shape, input_shape, (design.constant_shape..., filter_count), nothing,
                    LayerInitializer(
                        function (ins, outs, cons)
                            collection = [design.layer_init(view(ins, design.input_shape..., i), view(outs, design.output_shape..., i), view(cons, design.constant_shape..., i)) for i in 1:filter_count]
                            return function ()
                                for item in collection
                                    item()
                                end
                            end
                        end))
        end)
end

export MatrixWiseLayer
MatrixWiseLayer(layer::LayerGenerator) = MatrixWise(fun)


text:
abstract type AbstractBuffer end
abstract type AbstractNetworkBuffer <: AbstractBuffer end

"Store indexing on target device"
struct IndexingCache{IVT} <: AbstractBuffer
    index_pv::IVT
    index_pq::IVT
    index_ref::IVT
    index_generators::IVT
    index_pv_to_gen::IVT
    index_ref_to_gen::IVT
end

struct PolarNetworkState{VT} <: AbstractNetworkBuffer
    vmag::VT
    vang::VT
    pinj::VT
    qinj::VT
    pg::VT
    qg::VT
    balance::VT
    dx::VT
end


text:
using Completion: unify, Unifier

@testset "unify" begin
    @test unify(x, x) == Unifier()
    @test unify(x, y) == Unifier(x => y)
    @test unify(y, x) == Unifier(y => x)

    @test unify(@term(2), @term(2)) == Unifier()
    @test unify(@term(1), @term(2)) === nothing
    @test unify(@term(1), @term(1.0)) === nothing

    @test unify(x, f(y)) == Unifier(x => f(y))
    @test unify(f(y), x) == Unifier(x => f(y))
    @test unify(x, f(x)) == nothing
    @test unify(f(x), f(y)) == Unifier(x => y)
    @test unify(f(x), g(y)) == nothing
    @test unify(g(x, x), g(y, z)) == Unifier(x => z, y => z)
    @test unify(g(f(x), x), g(f(y), z)) == Unifier(x => z, y => z)
    @test unify(g(f(x), x), g(y, y)) == nothing
    @test unify(g(f(x), x), g(y, z)) == Unifier(y => f(z), x => z)
    @test unify(f(x), f(x, y)) == nothing
end


text:
@testset "Check For Test Param Function Tests" begin

  @test isdefined(Julz, :check_for_test_param) == true

end


text:
global radiance_datacube
global clouds_datacube
global mumbai_map
"""
The datacubes and the district level shapefile for Mumbai city are provided with the package as examples. These can be loaded using the function. The radiance datacube, cloud-free observations datacube and the shapefile of Mumbai districts are loaded as global variables.

    """
    function load_example()
        println("Loading sample data")
        package_path = pathof(SatelliteImages)
        path_len = length(package_path)
        assets_path = package_path[1:path_len-22] * "assets"
        map_path = assets_path * "/mumbai_map/mumbai_districts.shp"
        radiance_jld_path = assets_path * "/mumbai_ntl/datacube/mumbai_radiance.jld"
        clouds_jld_path = assets_path * "/mumbai_ntl/datacube/mumbai_clouds.jld"

        global radiance_datacube = load_datacube(radiance_jld_path)
        global clouds_datacube = load_datacube(clouds_jld_path)
        global mumbai_map = load_shapefile(map_path)
        """
        i) Distict level shapefile of Mumbai is loaded as mumbai_map.
        ii) Radiance datacube of Mumbai is loaded as radiance_datacube.
        iii) Cloud-free observations data is loaded as clouds_datacube.
        iv) MUMBAI_COORDINATE_SYSTEM should be used as the coordinate system.
        """
        println("
        Follow the tutorial on:

        https://xkdr.github.io/SatelliteImages.jl/tutorial.html")

end


text:
using Oceananigans.Advection

#####
##### Legend:
##### ρ  -> density fields
##### ρũ -> momentum fields
#####  ũ -> velocity fields
##### ρc̃ -> conservative tracer fields
#####  c̃ -> tracer fields
#####  K̃ -> diffusivity fields
#####

#####
##### Convert conservative variables to primitive variables
#####

@inline U_over_ρ(i, j, k, grid, ρ, ρu) = @inbounds ρu[i, j, k] / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
@inline V_over_ρ(i, j, k, grid, ρ, ρv) = @inbounds ρv[i, j, k] / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
@inline W_over_ρ(i, j, k, grid, ρ, ρw) = @inbounds ρw[i, j, k] / ℑzᵃᵃᶠ(i, j, k, grid, ρ)
@inline C_over_ρ(i, j, k, grid, ρ, ρc) = @inbounds ρc[i, j, k] / ρ[i, j, k]

#####
##### Kinetic energy
#####

@inline kinetic_energy(i, j, k, grid::AbstractGrid{FT}, ρ, ρũ) where FT =
    @inbounds FT(0.5) * (  (ℑxᶜᵃᵃ(i, j, k, grid, ρũ.ρu) / ρ[i, j, k])^2
                         + (ℑyᵃᶜᵃ(i, j, k, grid, ρũ.ρv) / ρ[i, j, k])^2
                         + (ℑzᵃᵃᶜ(i, j, k, grid, ρũ.ρw) / ρ[i, j, k])^2)

#####
##### Pressure work
#####

@inline p∂u∂x(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
       (  ℑxᶠᵃᵃ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
        * Axᵃᵃᶠ(i, j, k, grid) * U_over_ρ(i, j, k, grid, ρ, ρũ.ρu))

@inline p∂v∂y(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
       (  ℑyᵃᶠᵃ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
        * Ayᵃᵃᶠ(i, j, k, grid) * V_over_ρ(i, j, k, grid, ρ, ρũ.ρv))

@inline p∂w∂z(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
       (  ℑzᵃᵃᶠ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
        * Azᵃᵃᵃ(i, j, k, grid) * W_over_ρ(i, j, k, grid, ρ, ρũ.ρw))

@inline ∂ⱼpuⱼ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, p∂u∂x, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
           + δyᵃᶜᵃ(i, j, k, grid, p∂v∂y, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
           + δzᵃᵃᶜ(i, j, k, grid, p∂w∂z, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)))

#####
##### Tracer advection
#####

@inline tracer_flux_ρuc(i, j, k, grid, scheme, ρ, ρu, ρc) = advective_tracer_flux_x(i, j, k, grid, scheme, ρu, ρc) / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
@inline tracer_flux_ρvc(i, j, k, grid, scheme, ρ, ρv, ρc) = advective_tracer_flux_y(i, j, k, grid, scheme, ρv, ρc) / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
@inline tracer_flux_ρwc(i, j, k, grid, scheme, ρ, ρw, ρc) = advective_tracer_flux_z(i, j, k, grid, scheme, ρw, ρc) / ℑzᵃᵃᶠ(i, j, k, grid, ρ)

@inline div_ρUc(i, j, k, grid, scheme, density, momenta, ρc) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, tracer_flux_ρuc, scheme, density, momenta.ρu, ρc)
           + δyᵃᶜᵃ(i, j, k, grid, tracer_flux_ρvc, scheme, density, momenta.ρv, ρc)
           + δzᵃᵃᶜ(i, j, k, grid, tracer_flux_ρwc, scheme, density, momenta.ρw, ρc)))

#####
##### Diffusion
##### FIXME: need a better way to enforce boundary conditions on diffusive fluxes at rigid boundaries
#####

@inline diffusive_flux_x(i, j, k, grid, κᶠᶜᶜ, ρ, ρc) =
    ℑxᶠᵃᵃ(i, j, k, grid, ρ) * κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, C_over_ρ, ρ, ρc)

@inline diffusive_flux_y(i, j, k, grid, κᶜᶠᶜ, ρ, ρc) =
    ℑyᵃᶠᵃ(i, j, k, grid, ρ) * κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, C_over_ρ, ρ, ρc)

@inline diffusive_flux_z(i, j, k, grid, κᶜᶜᶠ, ρ, ρc) =
    ℑzᵃᵃᶠ(i, j, k, grid, ρ) * κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, C_over_ρ, ρ, ρc)

@inline diffusive_tvar_flux_x(i, j, k, grid, κᶠᶜᶜ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
    (ℑxᶠᵃᵃ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
     * κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, C_over_ρ, ρ, ρc))

@inline diffusive_tvar_flux_y(i, j, k, grid, κᶜᶠᶜ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
    (ℑyᵃᶠᵃ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
     * κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, C_over_ρ, ρ, ρc))

@inline diffusive_tvar_flux_z(i, j, k, grid, κᶜᶜᶠ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
    (ℑzᵃᵃᶠ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
     * κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, C_over_ρ, ρ, ρc))

@inline diffusive_pressure_flux_x(i, j, k, grid, κᶠᶜᶜ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) =
    (ℑxᶠᵃᵃ(i, j, k, grid, ρ) * κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid)
     * ∂xᶠᵃᵃ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))

@inline diffusive_pressure_flux_y(i, j, k, grid, κᶜᶠᶜ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) =
    (ℑyᵃᶠᵃ(i, j, k, grid, ρ) * κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid)
     * ∂yᵃᶠᵃ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))

@inline function diffusive_pressure_flux_z(i, j, k, grid::AbstractGrid{FT}, κᶜᶜᶠ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) where FT
    (k <= 1 || k > grid.Nz) && return zero(FT)
    return (ℑzᵃᵃᶠ(i, j, k, grid, ρ) * κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid)
            * ∂zᵃᵃᶠ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))
end

@inline function ∂ⱼDᶜⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tracer_index, ρ, ρc, args...)
    @inbounds κ = closure.κ[tracer_index]
    return (1/Vᵃᵃᶜ(i, j, k, grid)
            * (  δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, κ, ρ, ρc)
               + δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, κ, ρ, ρc)
               + δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, κ, ρ, ρc)))
end

@inline function ∂ⱼtᶜDᶜⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tvar_diag,
                         tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc, args...)

    @inbounds κ = closure.κ[tracer_index]
    return (1/Vᵃᵃᶜ(i, j, k, grid)
               * (  δxᶜᵃᵃ(i, j, k, grid, diffusive_tvar_flux_x, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)
                  + δyᵃᶜᵃ(i, j, k, grid, diffusive_tvar_flux_y, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)
                  + δzᵃᵃᶜ(i, j, k, grid, diffusive_tvar_flux_z, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)))
end

@inline function ∂ⱼDᵖⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tracer_index,
                       p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃, args ...)

    @inbounds κ = closure.κ[tracer_index]
    return (1/Vᵃᵃᶜ(i, j, k, grid)
               * (  δxᶜᵃᵃ(i, j, k, grid, diffusive_pressure_flux_x, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)
                  + δyᵃᶜᵃ(i, j, k, grid, diffusive_pressure_flux_y, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)
                  + δzᵃᵃᶜ(i, j, k, grid, diffusive_pressure_flux_z, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)))
end

#####
##### Momentum advection
##### Note the convention "momentum_flux_Ua" corresponds to the advection _of_ a _by_ U.
#####

@inline momentum_flux_ρuu(i, j, k, grid, scheme, ρ, ρu) = @inbounds momentum_flux_uu(i, j, k, grid, scheme, ρu, ρu) / ρ[i, j, k]
@inline momentum_flux_ρvv(i, j, k, grid, scheme, ρ, ρv) = @inbounds momentum_flux_vv(i, j, k, grid, scheme, ρv, ρv) / ρ[i, j, k]
@inline momentum_flux_ρww(i, j, k, grid, scheme, ρ, ρw) = @inbounds momentum_flux_ww(i, j, k, grid, scheme, ρw, ρw) / ρ[i, j, k]

@inline momentum_flux_ρuv(i, j, k, grid, scheme, ρ, ρu, ρv) = momentum_flux_uv(i, j, k, grid, scheme, ρv, ρu) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
@inline momentum_flux_ρuw(i, j, k, grid, scheme, ρ, ρu, ρw) = momentum_flux_uw(i, j, k, grid, scheme, ρw, ρu) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρvu(i, j, k, grid, scheme, ρ, ρu, ρv) = momentum_flux_vu(i, j, k, grid, scheme, ρu, ρv) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
@inline momentum_flux_ρvw(i, j, k, grid, scheme, ρ, ρv, ρw) = momentum_flux_vw(i, j, k, grid, scheme, ρw, ρv) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρwu(i, j, k, grid, scheme, ρ, ρu, ρw) = momentum_flux_wu(i, j, k, grid, scheme, ρu, ρw) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρwv(i, j, k, grid, scheme, ρ, ρv, ρw) = momentum_flux_wv(i, j, k, grid, scheme, ρv, ρw) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)

@inline div_ρuũ(i, j, k, grid, scheme, ρ, ρũ) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶠᵃᵃ(i, j, k, grid, momentum_flux_ρuu, scheme, ρ, ρũ.ρu)
           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρuv, scheme, ρ, ρũ.ρu, ρũ.ρv)
           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρuw, scheme, ρ, ρũ.ρu, ρũ.ρw)))

@inline div_ρvũ(i, j, k, grid, scheme, ρ, ρũ) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρvu, scheme, ρ, ρũ.ρu, ρũ.ρv)
           + δyᵃᶠᵃ(i, j, k, grid, momentum_flux_ρvv, scheme, ρ, ρũ.ρv)
           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρvw, scheme, ρ, ρũ.ρv, ρũ.ρw)))

@inline div_ρwũ(i, j, k, grid, scheme, ρ, ρũ) =
    (1/Vᵃᵃᶠ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρwu, scheme, ρ, ρũ.ρu, ρũ.ρw)
           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρwv, scheme, ρ, ρũ.ρv, ρũ.ρw)
           + δzᵃᵃᶠ(i, j, k, grid, momentum_flux_ρww, scheme, ρ, ρũ.ρw)))

#####
##### Viscous dissipation
#####

@inline strain_rate_tensor_ux(i, j, k, grid::AbstractGrid{FT}, ρ, ρu) where FT = FT(1/3) * ∂xᶜᵃᵃ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_vy(i, j, k, grid::AbstractGrid{FT}, ρ, ρv) where FT = FT(1/3) * ∂yᵃᶜᵃ(i, j, k, grid, V_over_ρ, ρ, ρv)
@inline strain_rate_tensor_wz(i, j, k, grid::AbstractGrid{FT}, ρ, ρw) where FT = FT(1/3) * ∂zᵃᵃᶜ(i, j, k, grid, W_over_ρ, ρ, ρw)

@inline strain_rate_tensor_uy(i, j, k, grid, ρ, ρu, ρv) = ∂yᵃᶠᵃ(i, j, k, grid, U_over_ρ, ρ, ρu) + ∂xᶠᵃᵃ(i, j, k, grid, V_over_ρ, ρ, ρv)
@inline strain_rate_tensor_uz(i, j, k, grid, ρ, ρu, ρw) = ∂zᵃᵃᶠ(i, j, k, grid, U_over_ρ, ρ, ρu) + ∂xᶠᵃᵃ(i, j, k, grid, W_over_ρ, ρ, ρw)
@inline strain_rate_tensor_vx(i, j, k, grid, ρ, ρv, ρu) = ∂xᶠᵃᵃ(i, j, k, grid, V_over_ρ, ρ, ρv) + ∂yᵃᶠᵃ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_vz(i, j, k, grid, ρ, ρv, ρw) = ∂zᵃᵃᶠ(i, j, k, grid, V_over_ρ, ρ, ρv) + ∂yᵃᶠᵃ(i, j, k, grid, W_over_ρ, ρ, ρw)
@inline strain_rate_tensor_wx(i, j, k, grid, ρ, ρw, ρu) = ∂xᶠᵃᵃ(i, j, k, grid, W_over_ρ, ρ, ρw) + ∂zᵃᵃᶠ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_wy(i, j, k, grid, ρ, ρw, ρv) = ∂yᵃᶠᵃ(i, j, k, grid, W_over_ρ, ρ, ρw) + ∂zᵃᵃᶠ(i, j, k, grid, V_over_ρ, ρ, ρv)

@inline viscous_flux_ux(i, j, k, grid, νᶜᶜᶜ, ρ, ρu) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Axᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_ux(i, j, k, grid, ρ, ρu)
@inline viscous_flux_vy(i, j, k, grid, νᶜᶜᶜ, ρ, ρv) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Ayᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_vy(i, j, k, grid, ρ, ρv)
@inline viscous_flux_wz(i, j, k, grid, νᶜᶜᶜ, ρ, ρw) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_wz(i, j, k, grid, ρ, ρw)

@inline viscous_flux_uy(i, j, k, grid, νᶠᶠᶜ, ρ, ρu, ρv) = ℑxyᶠᶠᵃ(i, j, k, grid, ρ) * νᶠᶠᶜ * Ayᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_uy(i, j, k, grid, ρ, ρu, ρv)
@inline viscous_flux_uz(i, j, k, grid, νᶠᶜᶠ, ρ, ρu, ρw) = ℑxzᶠᵃᶠ(i, j, k, grid, ρ) * νᶠᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_uz(i, j, k, grid, ρ, ρu, ρw)
@inline viscous_flux_vx(i, j, k, grid, νᶠᶠᶜ, ρ, ρv, ρu) = ℑxyᶠᶠᵃ(i, j, k, grid, ρ) * νᶠᶠᶜ * Axᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_vx(i, j, k, grid, ρ, ρv, ρu)
@inline viscous_flux_vz(i, j, k, grid, νᶜᶠᶠ, ρ, ρv, ρw) = ℑyzᵃᶠᶠ(i, j, k, grid, ρ) * νᶜᶠᶠ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_vz(i, j, k, grid, ρ, ρv, ρw)
@inline viscous_flux_wx(i, j, k, grid, νᶠᶜᶠ, ρ, ρw, ρu) = ℑxzᶠᵃᶠ(i, j, k, grid, ρ) * νᶠᶜᶠ * Axᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wx(i, j, k, grid, ρ, ρw, ρu)
@inline viscous_flux_wy(i, j, k, grid, νᶜᶠᶠ, ρ, ρw, ρv) = ℑyzᵃᶠᶠ(i, j, k, grid, ρ) * νᶜᶠᶠ * Ayᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wy(i, j, k, grid, ρ, ρw, ρv)

@inline ∂ⱼτ₁ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, closure.ν, ρ, ρũ.ρu)
           + δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, closure.ν, ρ, ρũ.ρu, ρũ.ρv)
           + δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, closure.ν, ρ, ρũ.ρu, ρũ.ρw)))

@inline ∂ⱼτ₂ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    (1/Vᵃᵃᶜ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, closure.ν, ρ, ρũ.ρv, ρũ.ρu)
           + δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, closure.ν, ρ, ρũ.ρv)
           + δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, closure.ν, ρ, ρũ.ρv, ρũ.ρw)))

@inline ∂ⱼτ₃ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    (1/Vᵃᵃᶠ(i, j, k, grid)
        * (  δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, closure.ν, ρ, ρũ.ρw, ρũ.ρu)
           + δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, closure.ν, ρ, ρũ.ρw, ρũ.ρv)
           + δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, closure.ν, ρ, ρũ.ρw)))

@inline u₁∂ⱼτ₁ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    U_over_ρ(i, j, k, grid, ρ, ρũ.ρu) * ∂ⱼτ₁ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline u₂∂ⱼτ₂ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    V_over_ρ(i, j, k, grid, ρ, ρũ.ρv) * ∂ⱼτ₂ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline u₃∂ⱼτ₃ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    W_over_ρ(i, j, k, grid, ρ, ρũ.ρw) * ∂ⱼτ₃ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline Q_dissipation(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
    (  ℑxᶜᵃᵃ(i, j, k, grid, u₁∂ⱼτ₁ⱼ, closure, ρ, ρũ, args...)
     + ℑyᵃᶜᵃ(i, j, k, grid, u₂∂ⱼτ₂ⱼ, closure, ρ, ρũ, args...)
     + ℑzᵃᵃᶜ(i, j, k, grid, u₃∂ⱼτ₃ⱼ, closure, ρ, ρũ, args...))


text:
module BlinkDisplay

using Blink, Media, Lazy, Requires

include("display.jl")
include("objects.jl")
include("mimes.jl")

export docs

docpane(url = "http://docs.julialang.org/en/latest/") =
  Window(d(:url => url,
           :alwaysOnTop => true,
           :width => 500,
           :height => 700))

docs() = docpane()

docs(search) = docpane("http://docs.julialang.org/en/latest/search/?q=$search")

end


text:
function make_logo()
circle = """
    oooo
 oooooooooo
oooooooooooo
oooooooooooo
 oooooooooo
    oooo    """

    # create circles
    green = Panel(
        RenderableText(circle, "#389826 bold");
        fit=true,
        style = "dim #389826",
        justify = :center,
        padding=(0, 0, 0, 0),
        title = "[italic]Made",
        title_style = "bold red",
    )
    red = Panel(
        RenderableText(circle, "#CB3C33 bold");
        fit=true,
        style = "dim #CB3C33",
        justify = :center,
        padding=(0, 0, 0, 0),
        subtitle = "[italic]with",
        subtitle_style = "bold #b656e3",
        subtitle_justify = :right,
    )
    purple = Panel(
        RenderableText(circle, "#9558B2 bold");
        fit=true,
        style = "dim #9558B2",
        justify = :center,
        padding=(0, 0, 0, 0),
        subtitle = "[italic]Term",
        subtitle_style = "bold #389826",
    )

    indigo = "#42A5F5"

    as_code(x) = "[orange1 italic]`$x`[/orange1 italic]"

    main = TextBox(
        """Term.jl is a [#9558B2]Julia[/#9558B2] package for creating styled terminal outputs.

        Term provides a simple [italic green4 bold]markup language[/italic green4 bold] toadd [bold bright_blue]color[/bold bright_blue] and [bold underline]styles[/bold underline] to your text.
        More complicated text layout can be created using $(as_code("Renderable")) objects such
        as $(as_code("Panel")) and $(as_code("TextBox")).
        These can also be nested and stacked to create[italic pink3]fancy[/italic pink3]and [underline]informative[/underline] terminal ouputs for your Julia code""";
        width = 35,
    )

    # create "spacers" and stack renderables
    hspacer = Spacer(green.measure.w / 2 + 1, green.measure.h)
    line = Spacer(green.measure.w * 2 + 6, 1)
    circles =  line / (hspacer * green * hspacer) / (red * Spacer(2, purple.measure.h) * purple)

    content =
        circles *
        vLine(main.measure.h; style = indigo * " dim") *
        main

    # add second message
    second_message = RenderableText(
        """
        Term.jl can also be used to create [underline]fancy[/underline] $(as_code("logging")) and $(as_code("error")) messages.

        Check the examples and documentation for more information!
        Term.jl is under [bold]active[/bold] development, get in touch for questions or ideas on how to  improve it!""";
        width = content.measure.w,
    )
    hline = hLine(content.measure.w; style = indigo * " dim")
    content = content / hline / second_message

    # add a final message
    msg = """[#75b6e0]Term.jl is based on the [underline]Python[/underline] library [orange1 italic]Rich[/orange1 italic] by [/#75b6e0]Will McGugan.
                        [dim]https://github.com/Textualize/rich[/dim]"""

    logo =
        Panel(
            content;
            title = "Term.jl",
            title_style = "bold $indigo",
            style = "dim $indigo",
            subtitle = "https://github.com/FedeClaudi/Term.jl",
            subtitle_justify = :right,
            subtitle_style = "dim",
            fit = true,
        ) / "" / (Spacer(12, 2) * msg)

    return logo
end

# print(make_logo())


text:
# Huber loss function

using LinearAlgebra

export HuberLoss

"""
**Huber loss**

    HuberLoss(ρ=1, μ=1)

Returns the function
```math
f(x) = \\begin{cases}
    \\tfrac{μ}{2}\\|x\\|^2 & \\text{if}\\ \\|x\\| ⩽ ρ \\\\
    ρμ(\\|x\\| - \\tfrac{ρ}{2}) & \\text{otherwise},
\\end{cases}
```
where `ρ` and `μ` are positive parameters.
"""
struct HuberLoss{R <: Real, S <: Real} <: ProximableFunction
    rho::R
    mu::S
    function HuberLoss{R, S}(rho::R, mu::S) where {R <: Real, S <: Real}
        if rho <= 0 || mu <= 0
            error("parameters rho and mu must be positive")
        else
            new(rho, mu)
        end
    end
end

is_convex(f::HuberLoss) = true
is_smooth(f::HuberLoss) = true

HuberLoss(rho::R=1, mu::S=1) where {R <: Real, S <: Real} = HuberLoss{R, S}(rho, mu)

function (f::HuberLoss)(x::AbstractArray{T}) where {R <: Real, T <: Union{R, Complex{R}}}
    normx = norm(x)
    if normx <= f.rho
        return f.mu / R(2) * normx^2
    else
        return f.rho * f.mu * (normx - f.rho / R(2))
    end
end

function gradient!(y::AbstractArray{T}, f::HuberLoss, x::AbstractArray{T}) where {R <: Real, T <: Union{R, Complex{R}}}
    normx = norm(x)
    if normx <= f.rho
        y .= f.mu .* x
        v = f.mu / R(2) * normx^2
    else
        y .= (f.mu * f.rho) / normx .* x
        v = f.rho * f.mu * (normx - f.rho / R(2))
    end
    return v
end

function prox!(y::AbstractArray{T}, f::HuberLoss, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: Union{R, Complex{R}}}
    normx = norm(x)
    mugam = f.mu*gamma
    scal = (1 - min(mugam / (1 + mugam), mugam * f.rho / normx))
    for k in eachindex(y)
        y[k] = scal*x[k]
    end
    normy = scal*normx
    if normy <= f.rho
        return f.mu / R(2) * normy^2
    else
        return f.rho * f.mu * (normy - f.rho / R(2))
    end
end

fun_name(f::HuberLoss) = "Huber loss"
fun_dom(f::HuberLoss) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::HuberLoss) = "x ↦ (μ/2)||x||² if ||x||⩽ρ, μρ(||x||-ρ/2) otherwise"
fun_params(f::HuberLoss) = string("ρ = $(f.rho), μ = $(f.mu)")

function prox_naive(f::HuberLoss, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: Union{R, Complex{R}}}
    y = (1 - min(f.mu * gamma / (1 + f.mu * gamma), f.mu * gamma * f.rho / norm(x))) * x
    if norm(y) <= f.rho
        return y, f.mu / R(2) * norm(y)^2
    else
        return y, f.rho * f.mu * (norm(y) - f.rho / R(2))
    end
end


text:
include("./util.jl")

RegOrValue = Union{AbstractString,UInt16}

struct LShift
    value::RegOrValue
    amount::UInt16
end
struct RShift
    value::RegOrValue
    amount::UInt16
end
struct And
    lhs::RegOrValue
    rhs::RegOrValue
end
struct Or
    lhs::RegOrValue
    rhs::RegOrValue
end
struct Not
    value::RegOrValue
end
struct Store
    value::RegOrValue
end

Op = Union{LShift,RShift,And,Or,Not,Store}

struct Instruction
    op::Op
    target::String
end

function reg_or_value(s::AbstractString)::RegOrValue
    tried = tryparse(UInt16, s)
    if isnothing(tried)
        s
    else
        tried
    end
end

function Instruction(s::AbstractString)::Instruction
    re =
        r"(?:(?:(?<lhs>[a-z0-9]+) (?<op>AND|OR|RSHIFT|LSHIFT) (?<rhs>[a-z0-9]+))|(?:(?<nop>NOT) (?<nvalue>[a-z0-9]+))|(?<value>[a-z0-9]+)) -> (?<target>[a-z]+)"

    m = match(re, s)
    if isnothing(m)
        throw(DomainError("Unknown instruction $s"))
    end
    op = if isnothing(m[:op]) && isnothing(m[:nop])
        Store(reg_or_value(m[:value]))
    elseif m[:nop] == "NOT"
        Not(reg_or_value(m[:nvalue]))
    elseif m[:op] == "AND"
        And(reg_or_value(m[:lhs]), reg_or_value(m[:rhs]))
    elseif m[:op] == "OR"
        Or(reg_or_value(m[:lhs]), reg_or_value(m[:rhs]))
    elseif m[:op] == "RSHIFT"
        RShift(reg_or_value(m[:lhs]), reg_or_value(m[:rhs]))
    elseif m[:op] == "LSHIFT"
        LShift(reg_or_value(m[:lhs]), reg_or_value(m[:rhs]))
    else
        throw(DomainError("Unknown op $(m[:op])"))
    end
    Instruction(op, m[:target])
end

instructions = Dict()
cache = Dict()

inputs = map(Instruction, split(util.load_input(7), '\n'))

for i in inputs
    instructions[i.target] = i
end

function resolve(register::AbstractString)::UInt16
    if register in keys(cache)
        cache[register]
    else
        result = compute(instructions[register].op)
        cache[register] = result
        result
    end
end

resolve(value::UInt16)::UInt16 = value

compute(op::Store)::UInt16 = resolve(op.value)

compute(op::Not)::UInt16 = ~resolve(op.value)

compute(op::Or)::UInt16 = resolve(op.lhs) | resolve(op.rhs)

compute(op::And)::UInt16 = resolve(op.lhs) & resolve(op.rhs)

compute(op::RShift)::UInt16 = resolve(op.value) >> op.amount

compute(op::LShift)::UInt16 = resolve(op.value) << op.amount

reg_a = resolve("a")
println("Part 1: $reg_a")

for k in keys(cache)
    delete!(cache, k)
end
instructions["b"] = Instruction(Store(UInt16(16076)), "b")

reg_a = resolve("a")
println("Part 2: $reg_a")


text:
using IRTools, MacroTools, InteractiveUtils, Test
using IRTools: @dynamo, IR, meta, isexpr, xcall, self, insertafter!, recurse!,
  argument!, return!, func, var, functional

@dynamo roundtrip(a...) = IR(a...)

@dynamo function passthrough(a...)
  ir = IR(a...)
  ir == nothing && return
  recurse!(ir)
  return ir
end

add(a, b) = a+b
@test roundtrip(add, 2, 3) == 5
@test passthrough(add, 2, 3) == 5
@test passthrough(sin, 1.0) == sin(1.0)

f_returnnode(::T) where {T} = T
@test passthrough(f_returnnode, 1) == Int
f_gotoifnot(::Val{x}) where {x} = x ? 1 : 0
@test passthrough(f_gotoifnot, Val(true)) == 1

@test @code_ir(passthrough, add(2, 3)) isa IR

IRTools.refresh(roundtrip)

add(a, b) = a*b
@test roundtrip(add, 2, 3) == 6

relu(x) = x > 0 ? x : 0
@test roundtrip(relu, 1) == 1
@test passthrough(relu, 1) == 1

foo(x) = y = x > 0 ? x + 1 : x - 1
@test roundtrip(foo, 1) == 2
@test roundtrip(foo, -1) == -2
@test passthrough(foo, 1) == 2

@test passthrough(() -> [1, 2, 3]) == [1, 2, 3]

function err(f)
  try
    f()
  catch e
    e
  end
end

@test passthrough(err, () -> 2+2) == 4
@test passthrough(err, () -> 0//0) isa ArgumentError

function err2(f)
    x = 1
    y = 1
    try
        y = f()
    catch e
      x, y
    end
    x, y
end

@test passthrough(err2, () -> 2+2) == (1, 4)
@test passthrough(err2, () -> 0//0) == (1, 1)

function err3(f)
  y = 1
  try
    y = f()
  catch
  end
  return y
end

@test passthrough(err3, () -> 2+2) == 4
@test_broken passthrough(err3, () -> 0//0) == 1

@dynamo function mullify(a...)
  ir = IR(a...)
  ir == nothing && return
  ir = MacroTools.prewalk(ir) do x
    x isa GlobalRef && x.name == :(*) && return GlobalRef(Base, :+)
    return x
  end
  for (x, st) in ir
    isexpr(st.expr, :call) || continue
    ir[x] = xcall(Main, :mullify, st.expr.args...)
  end
  return ir
end

@test mullify(prod, [5, 10]) == 15

@dynamo err(a...) = 0//0

@test_throws IRTools.Inner.CompileError err(+, 2, 3)

mutable struct Context
  calls::Int
end

@dynamo function (cx::Context)(a...)
  ir = IR(a...)
  ir == nothing && return
  recurse!(ir)
  calls = pushfirst!(ir, xcall(:getfield, self, QuoteNode(:calls)))
  calls = insertafter!(ir, calls, xcall(:+, calls, 1))
  insertafter!(ir, calls, xcall(:setfield!, self, QuoteNode(:calls), calls))
  return ir
end

@code_lowered Context(0)(add, 2, 3)

cx = Context(0)
@test cx(add, 2, 3.0) == 6
@test cx.calls > 5

ir = IR()
x = argument!(ir)
y = push!(ir, xcall(:*, x, x))
return!(ir, y)

@test IRTools.evalir(ir, 5) == 25

ir = IR()
x = argument!(ir)
y = argument!(ir)
z = push!(ir, xcall(:*, x, y))
return!(ir, z)

@test IRTools.evalir(ir, 5, 3) == 15

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

ir = @code_ir pow(2, 3)

ir[var(8)] = xcall(:+, var(5), var(2))

mul = func(ir)

@test mul(nothing, 10, 3) == 31

@dynamo function ir_add(_, _)
  ir = IR()
  args = argument!(ir)
  x = push!(ir, xcall(:getindex, args, 1))
  y = push!(ir, xcall(:getindex, args, 2))
  return!(ir, xcall(:+, x, y))
end

@test ir_add(5, 2) == 7

@dynamo function test_lambda(x)
  λ = IR()
  self = argument!(λ)
  y = argument!(λ)
  x = push!(λ, xcall(:getindex, self, 1))
  return!(λ, xcall(:+, x, y))
  ir = IR()
  args = argument!(ir)
  x = push!(ir, xcall(:getindex, args, 1))
  return!(ir, Expr(:lambda, λ, x))
end

let
  f = test_lambda(3)
  @test f(6) == 9
end

anf(f::Core.IntrinsicFunction, args...) = f(args...)

@dynamo function anf(args...)
  ir = IR(args...)
  ir == nothing && return
  functional(recurse!(ir, anf))
end

@test anf(pow, 2, 3) == 8

function f94(x)
    i = x
    while i > 0
        i -= 1
    end
    return x
end
@test f94(1) == 1
@test passthrough(f94, 1) == 1


text:
_width = 500
_height = 500

"Return the current canvas size (width and height)."
canvassize() = _width, _height

"Set the canvas size (width and height)."
canvassize(w::Integer, h::Integer) = global _width, _height = w, h

width() = _width
height() = _height
width(w) = global _with = w
height(h) = global _height = h

_xmin = -2
_xmax = 2
_ymin = -2
_ymax = 2

for fun = (:xmin, :xmax, :ymin, :ymax)
    var = Symbol("_",fun)
    @eval begin
        $fun() = $var
        $fun(t::Real) = global $var = t
    end # eval
end # for

"Return the current X limits of the rectangular region of drawing."
xlims() = _xmin, _xmax

"Set the X limits of the rectangular region of drawing."
xlims(a::Real, b::Real) = global _xmin, _xmax = a, b

"Return the current Y limits of the rectangular region of drawing."
ylims() = _ymin, _ymax

"Set the Y limits of the rectangular region of drawing."
ylims(a::Real, b::Real) = global _ymin, _ymax = a, b

"Set the rectangular region limits."
rectregion(x0::Real, x1::Real, y0::Real,y1 ::Real) =
global _xmin, _xmax, _ymin, _ymax = x0, x1, y0, y1

"Create a square (**Tuple** of four real numbers) to set the rectangular region."
square(l::Real, x0::Real=0, y0::Real=0) = x0-l, x0+l, y0-l, y0+l

"Check if (x,y) is inside of the rectangular region."
insiderectregion(x::Real, y::Real) = _xmin < x < _xmax && _ymin < y < _ymax

"Check if z is inside of the rectangular region."
insiderectregion(z::Complex) = insiderectregion(z.re, z.im)

"Check if x is inside of the rectangular region."
insiderectregion(x::Real) = insiderectregion(x,0)

"Adjust canvas size to rectangular region to have the same aspect ratio."
function adjustcanvas()
    aspectratio = (_ymax-_ymin)/(_xmax-_xmin)
    if aspectration <= 1
        global _width = Int(ceil(_height/aspectratio))
    else
        global _height = Int(ceil(_width*aspectratio))
    end # if
end

"Adjust rectangular region to canvas size to have the same aspect ratio."
function adjustrectregion() # INCORRECT!
    aspectratio = _height/_width
    if aspectration <= 1
        x0 = (_xmin+_xmax)/2
        Δx = (_xmax-_xmin)/aspectratio
        global _xmin = x0 - Δx
        global _xmax = x0 + Δx
    else
        y0 = (_ymin+_ymax)/2
        Δy = (_ymax-_ymin)*aspectratio
        global _ymin = y0 - Δy
        global _ymax = y0 + Δy
    end # if
end

"Convert plane coordinates to image coordinates."
pointtopixel(x::Real, y::Real) =
    Integer(ceil(_height*(_ymax-y)/(_ymax-_ymin))), # j
    Integer(ceil(_width*(x-_xmin)/(_xmax-_xmin)))   # i

"Convert complex plane coordinates to image coordinates."
pointtopixel(z::Number) = pointtopixel(real(z), imag(z))

"Convert plane coordinates to canvas coordinates."
pointtocanvas(x::Real, y::Real) =
    _width*(x-_xmin)/(_xmax-_xmin), _height*(_ymax-y)/(_ymax-_ymin)

"Convert complex plane coordinates to canvas coordinates."
pointtocanvas(z::Number) = pointtocanvas(real(z), imag(z))

"Convert image coordinates to plane coordinates."
pixeltopoint(i::Integer, j::Integer) =
    _xmin+(i-1)*(_xmax-_xmin)/(_width-1), _ymax-(j-1)*(_ymax-_ymin)/(_height-1)

"Convert canvas coordinates to plane coordinates."
canvastopoint(i::Real, j::Real) =
    _xmin+(i-1)*(_xmax-_xmin)/(_width-1), _ymax-(j-1)*(_ymax-_ymin)/(_height-1)


text:
# Is this usable data?
function isconcrete{T<:Number}(x::T)
    !isna(x) && isfinite(x)
end

function isconcrete(x::MathConst)
    return true
end

function isconcrete(x)
    !isna(x)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    n = 0
    for x in xs
        if isconcrete(x)
            n += 1
        end
    end
    n
end

function concrete_length{T}(xs::DataArray{T})
    n = 0
    for i = 1:length(xs)
        if !xs.na[i] && isfinite(xs.data[i]::T)
            n += 1
        end
    end
    n
end

function concrete_length(xs::Iterators.Chain)
    n = 0
    for obj in xs.xss
        n += concrete_length(obj)
    end
    n
end

function concrete_minimum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_min = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_max = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function concrete_minmax{T<:Real}(xs, xmin::T, xmax::T)
    if eltype(xs) <: Base.Callable
        return xmin, xmax
    end

    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T}(xs, xmin::T, xmax::T)
    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real}(xs::Iterators.Chain, xmin::T, xmax::T)
    for obj in xs.xss
        xmin, xmax = concrete_minmax(obj, xmin, xmax)
    end
    xmin, xmax
end

function nonzero_length(xs)
    n = 0
    for x in xs
        if x != 0
            n += 1
        end
    end
    n
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit{T}(a::T, b::T)
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!{T}(a::T, b::T)
    for field in T.names
        aval = getfield(a, field)
        bval = getfield(b, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if aval === nothing || aval === string || aval == showoff
            setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = is(u, nothing)
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has{T,N}(xs::AbstractArray{T,N}, y::T)
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

Maybe(T) = Union(T, Nothing)

function lerp(x::Float64, a::Float64, b::Float64)
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::String)
    s = replace(s, r"<[^>]*>", "")
    s = replace(s, r"\s", "_")
    s
end

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable{T <: FloatingPoint}(x::T) = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    @assert n > 1

    step = (b - a) / (n - 1)
    xs = Array(typeof(a + step), n)
    for i in 1:n
        xs[i] = a + (i-1) * step
    end

    df = DataFrame(xs, map(f, xs))
    # NOTE: 'colnames!' is the older deprecated name. 'names!' was also defined
    # but threw an error.
    try
        names!(df, [:x, :f_x])
    catch
        colnames!(df, ["x", "f_x"])
    end
    df
end

# Construct a jscall to store arbitrary data in the element
function jsdata(key::String, value::String, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        data("$(escape_string(key))", $(value))
        """, arg)
end

# Construct jscall properties to store arbitrary data in plotroot elements.
function jsplotdata(key::String, value::String, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        plotroot().data("$(escape_string(key))", $(value))
        """, arg)
end

function svg_color_class_from_label(label::String)
    return @sprintf("color_%s", escape_id(label))
end

# TODO: Delete when 0.3 compatibility is dropped
if VERSION < v"0.4-dev"
    using Dates

    function Showoff.showoff{T <: Date}(ds::AbstractArray{T}, style=:none)
        const month_names = [
            "January", "February", "March", "April", "May", "June", "July",
            "August", "September", "October", "November", "December"
        ]

        const month_abbrevs = [
            "Jan", "Feb", "Mar", "Apr", "May", "Jun",
            "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"
        ]

        day_all_1   = all(map(d -> Dates.day(d) == 1, ds))
        month_all_1 = all(map(d -> Dates.month(d) == 1, ds))
        years = Set()
        for d in ds
            push!(years, Dates.year(d))
        end

        buf = IOBuffer()
        labels = Array(String, length(ds))
        if day_all_1 && month_all_1
            # only label years
            for (i, d) in enumerate(ds)
                print(buf, Dates.year(d))
                labels[i] = takebuf_string(buf)
            end
        elseif day_all_1
            # label months and years
            for (i, d) in enumerate(ds)
                if d == ds[1] || Dates.month(d) == 1
                    print(buf, month_abbrevs[Dates.month(d)], " ", Dates.year(d))
                else
                    print(buf, month_abbrevs[Dates.month(d)])
                end
                labels[i] = takebuf_string(buf)
            end
        else
            for (i, d) in enumerate(ds)
                if d == ds[1] || (Dates.month(d) == 1 && Dates.day(d) == 1)
                    print(buf, month_abbrevs[Dates.month(d)], " ", Dates.day(d), " ", Dates.year(d))
                elseif Dates.day(d) == 1
                    print(buf, month_abbrevs[Dates.month(d)], " ", Dates.day(d))
                else
                    print(buf, Dates.day(d))
                end
                labels[i] = takebuf_string(buf)
            end
        end

        return labels
    end
else
    using Base.Dates
end

if !method_exists(/, (Dates.Day, Dates.Day))
    /(a::Dates.Day, b::Dates.Day) = a.value / b.value
end

if !method_exists(/, (Dates.Day, Real))
    /(a::Dates.Day, b::Real) = Dates.Day(iround(a.value / b))
end

#if !method_exists(*, (FloatingPoint, Dates.Day))
    *(a::FloatingPoint, b::Dates.Day) = Dates.Day(iround(a * b.value))
    *(a::Dates.Day, b::FloatingPoint) = b * a
#end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
struct _2UKP
    nSize::Int
    P1::Vector{Int}
    P2::Vector{Int}
    W::Vector{Int}
    C::Int
end

Base.show(io::IO, id::_2UKP) = print("Bi-Objective Knapsack Problem with $(id.nSize) variables.")

set2UKP(n::Int, p1::Vector{Int}, p2::Vector{Int}, w::Vector{Int}, c::Int) = begin
    @assert n == length(p1) == length(p2) == length(w)
    _2UKP(n, p1, p2, w, c)
end
set2UKP(p1::Vector{Int}, p2::Vector{Int}, w::Vector{Int}, c::Int) = set2UKP(length(p1), p1, p2, w, c)

struct UKPsolver
    solve::Function #::Function(id::_2LAP) -> ...
end

function vSolve(id::_2UKP, solver::UKPsolver = UKP_Jorge2010())
    solver.solve(id)
end

function UKP_Jorge2010(output::Bool = true)::UKPsolver
    f = (id::_2UKP) -> begin
        pb = Bi01KP.problem(id.P1, id.P2, id.W, id.C)
        return Bi01KP.solve(pb, output)
    end

    return UKPsolver(f)
end

function load2UKP(fname)
    f = open(fname)
    n = parse(Int, readline(f))
    P = parse(Int, readline(f))
    nbC = parse(Int, readline(f))
    p1 = parse.(Int, split(readline(f)))
    p2 = parse.(Int, split(readline(f)))
    w = parse.(Int, split(readline(f)))
    c = parse(Int, readline(f))
    _2UKP(n, p1, p2, w, c)
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
# This script provides an example of simulating changepoint data,
# using PELT to find these changepoints, and plotting the results
#
# Requires Gadfly plotting package

using Distributions, Gadfly
using Changepoints

n = 1000        # Sample size
λ = 100         # Frequencey of changes

μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)

# Find changepoints via PELT
pelt_cps, pelt_cost = @PELT sample Normal(?, σ)

# Find changepoints via BS
bs_cps = @BS sample Normal(?, σ)

# Plot data with true changepoints
p = plot(sample, cps)