text:
### FlatMap and FlatFourier types
struct FlatMap{P<:Flat,T<:Real,M<:AbstractRank2or3Array{T}} <: Field{Map,S0,P,T}
    Ix :: M
end
struct FlatFourier{P<:Flat,T<:Real,M<:AbstractRank2or3Array{Complex{T}}} <: Field{Fourier,S0,P,Complex{T}}
    Il :: M
end
const FlatS0{P,T,M} = Union{FlatMap{P,T,M},FlatFourier{P,T,M}}

### convenience constructors
for (F, X, T) in [
    (:FlatMap,     :Ix, :T),
    (:FlatFourier, :Il, :(Complex{T})),
]
    doc = """
        # main constructor:
        $F(
            $X::AbstractArray; $(F==:FlatFourier ? "\n        Nside, # required, size of the map in pixels" : "")
            θpix,  # optional, resolution in arcmin (default: 1)
            ∂mode, # optional, fourier∂ or map∂ (default: fourier∂)
        )

        # more low-level:
        $F{P}($X::AbstractArray) # specify pixelization P explicilty
        $F{P,T}($X::AbstractArray) # additionally, convert elements to type $T
        $F{P,T,M<:AbstractArray{$T}}($X::M) # specify everything explicilty

    Construct a `$F` object. The top form of the constructor is most convenient
    for interactive work, while the others may be more useful for low-level code.
    """
    @eval begin
        @doc $doc $F
        $F{P}($X::M) where {P,T,M<:AbstractRank2or3Array{$T}} = $F{P,T,M}($X)
        $F{P,T}($X::AbstractRank2or3Array) where {P,T} = $F{P}($T.($X))
    end
    T!=:T && @eval $F{P}($X::M) where {P,T,M<:AbstractRank2or3Array{T}} = $F{P,T}($X)
end
FlatMap(Ix::AbstractRank2or3Array; kwargs...) =
    FlatMap{Flat(Nside=(size(Ix,1)==size(Ix,2) ? size(Ix,1) : size(Ix)[1:2]), D=size(Ix,3); kwargs...)}(Ix)
FlatFourier(Il::AbstractRank2or3Array; Nside, kwargs...) =
    FlatFourier{Flat(Nside=Nside, D=size(Il,3); kwargs...)}(Il)

### array interface
# Base.size and Base.lastindex refer to the splayed-out vector representation of the field:
size(f::FlatS0) = (length(firstfield(f)),)
lastindex(f::FlatS0, i::Int) = lastindex(f.Ix, i)
# content_size and content_ndims refer to the actual array holding the field content:
content_size(::Type{<:FlatMap{    <:Flat{N,<:Any,<:Any,D}}}) where {N,D} = (N .* (1,1)         ..., (D==1 ? () : (D,))...)
content_size(::Type{<:FlatFourier{<:Flat{N,<:Any,<:Any,D}}}) where {N,D} = (N .÷ (2,1) .+ (1,0)..., (D==1 ? () : (D,))...)
content_ndims(::Type{<:FlatS0{<:Flat{<:Any,<:Any,<:Any,D}}}) where {D} = D==1 ? 3 : 2
@propagate_inbounds @inline getindex(f::FlatS0, I...) = getindex(firstfield(f), I...)
@propagate_inbounds @inline setindex!(f::FlatS0, X, I...) = (setindex!(firstfield(f), X, I...); f)
adapt_structure(to, f::F) where {P,F<:FlatS0{P}} = basetype(F){P}(adapt(to,firstfield(f)))
adapt_structure(::Type{T}, f::F) where {T<:Union{Float32,Float64},P,F<:FlatS0{P}} = T(f)
function similar(f::F,::Type{T},dims::Dims) where {P,F<:FlatS0{P},T<:Number}
    @assert size(f)==dims "Tried to make a field similar to $F but dims should have been $(size(f)), not $dims."
    basetype(F){P}(similar(firstfield(f),T))
end

### broadcasting
struct FlatS0Style{F,M} <: AbstractArrayStyle{1} end
(::Type{FS})(::Val{1}) where {FS<:FlatS0Style} = FS()
(::Type{FS})(::Val{2}) where {FS<:FlatS0Style} = error("Broadcast expression would create a dense Field operator.")
@generated BroadcastStyle(::Type{F}) where {P,T,M,F<:FlatS0{P,T,M}} = FlatS0Style{basetype(F){P},basetype(M)}()
# both orders needed bc of https://github.com/JuliaLang/julia/pull/35948:
BroadcastStyle(S::FlatS0Style{<:FlatS0{Flat{N,θ,∂m,D}}},  ::FlatS0Style{<:FlatS0{Flat{N,θ,∂m,1}}}) where {N,θ,∂m,D} = S
BroadcastStyle( ::FlatS0Style{<:FlatS0{Flat{N,θ,∂m,1}}}, S::FlatS0Style{<:FlatS0{Flat{N,θ,∂m,D}}}) where {N,θ,∂m,D} = S
BroadcastStyle(S::FieldTupleStyle, ::FlatS0Style) = S
BroadcastStyle(S::FieldOrOpArrayStyle, ::FlatS0Style) = S
BroadcastStyle(S1::FlatS0Style{<:Field{B1}}, S2::FlatS0Style{<:Field{B2}}) where {B1,B2} =
    invalid_broadcast_error(B1,S1,B2,S2)
instantiate(bc::Broadcasted{<:FlatS0Style}) = bc
similar(::Broadcasted{FS}, ::Type{T}) where {T<:Number,FS<:FlatS0Style} = similar(FS,T)
similar(::Type{FlatS0Style{F,M}}, ::Type{T}) where {F<:FlatS0,M,T<:Number} =
    F(basetype(M){eltype(F{real(T)})}(undef,content_size(F)))
@inline preprocess(dest::F, bc::Broadcasted) where {F<:FlatS0} =
    Broadcasted{DefaultArrayStyle{content_ndims(F)}}(bc.f, preprocess_args(dest, bc.args), map(OneTo,content_size(F)))
preprocess(dest::F, arg) where {F<:FlatS0} = broadcastable(F, arg)
broadcastable(::Type{<:FlatS0}, f::FlatS0) = firstfield(f)
broadcastable(::Type{<:FlatS0{<:Flat,T}}, r::Real) where {T} = convert(T,r)
broadcastable(::Any, x) = x
@inline function Broadcast.copyto!(dest::FlatS0, bc::Broadcasted{Nothing})
    bc′ = preprocess(dest, bc)
    @simd for I in eachindex(bc′)
        @inbounds dest[I] = bc′[I]
    end
    return dest
end

### basis conversion
Fourier(f::FlatMap{P}) where {P} = FlatFourier{P}(fieldinfo(f).FFT * f.Ix)
Map(f::FlatFourier{P}) where {P} = FlatMap{P}(fieldinfo(f).FFT \ f.Il)

### inplace conversion
Fourier(f′::FlatFourier, f::FlatMap) =  (mul!(f′.Il, fieldinfo(f).FFT, f.Ix); f′)
Map(f′::FlatMap, f::FlatFourier)     = (ldiv!(f′.Ix, fieldinfo(f).FFT, maybecopy(f.Il)); f′)
# need this for FFTW (but not MKL) see https://github.com/JuliaLang/julia/issues/37000
@static FFTW.fftw_vendor==:fftw ? maybecopy(x) = copy(x) : maybecopy(x) = x

### properties
getproperty(f::FlatS0, s::Symbol) = getproperty(f,Val(s))
getproperty(f::FlatS0, ::Val{s}) where {s} = getfield(f,s)
getproperty(f::FlatS0, ::Val{:I}) = f
function getindex(f::FlatS0, k::Symbol; full_plane=false)
    maybe_unfold = full_plane ? x->unfold(x,fieldinfo(f).Ny) : identity
    @match k begin
        :I  => f
        :Ix => Map(f).Ix
        :Il => maybe_unfold(Fourier(f).Il)
        _   => throw(ArgumentError("Invalid FlatS0 index: $k"))
    end
end

### dot products
# do in Map space for simplicity, and use sum_kbn to reduce roundoff error
dot(a::FlatS0{<:Flat{N,θ}}, b::FlatS0{<:Flat{N,θ}}) where {N,θ} = batch(sum_kbn(Map(a).Ix .* Map(b).Ix, dims=(1,2)))

### isapprox
≈(a::F, b::F) where {P,T,F<:FlatS0{P,T}} = all(.≈(a[:], b[:], atol=sqrt(eps(T)), rtol=sqrt(eps(T))))

### simulation and power spectra
function white_noise(rng::AbstractRNG, ::Type{F}) where {N,P<:Flat{N},T,M,F<:FlatS0{P,T,M}}
    FlatMap{P}(randn!(rng, basetype(M){T}(undef, content_size(FlatMap{P}))))
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S0}, Cℓ::InterpolatedCℓs; units=fieldinfo(P).Ωpix) where {P,T}
    Diagonal(FlatFourier{P}(Cℓ_to_2D(P,T,Cℓ)) / units)
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S0}, (Cℓ, ℓedges, θname)::Tuple; units=fieldinfo(P).Ωpix) where {P,T}
    C₀ = Cℓ_to_Cov(P, T, S0, Cℓ, units=units)
    Cbins = Diagonal.(MidPasses(ℓedges) .* [diag(C₀)])
    BinRescaledOp(C₀,Cbins,θname)
end

function cov_to_Cℓ(L::DiagOp{<:FlatS0{P}}; units=fieldinfo(P).Ωpix) where {P}
    ii = sortperm(fieldinfo(L.diag).kmag[:])
    InterpolatedCℓs(fieldinfo(L.diag).kmag[ii], real.(unfold(L.diag.Il, fieldinfo(L.diag).Ny))[ii] * units, concrete=false)
end

function get_Cℓ(f::FlatS0{P}, f2::FlatS0{P}=f; Δℓ=50, ℓedges=0:Δℓ:16000, Cℓfid=ℓ->1, err_estimate=false) where {P}
    @unpack Nx, Ny ,Δx,kmag = fieldinfo(f)
    α = Nx*Ny/Δx^2

    # faster to excise unused parts:
    kmask = (kmag .> minimum(ℓedges)) .&  (kmag .< maximum(ℓedges))
    L = Float64.(kmag[kmask])
    CLobs = real.(dot.(unfold(Float64(f)[:Il],Ny)[kmask], unfold(Float64(f2)[:Il],Ny)[kmask])) ./ α
    w = @. nan2zero((2*Cℓfid(L)^2/(2L+1))^-1)

    sum_in_ℓbins(x) = fit(Histogram, L, Weights(x), ℓedges).weights

    local A, Cℓ, ℓ, N, Cℓ²
    Threads.@sync begin
        Threads.@spawn A  = sum_in_ℓbins(w)
        Threads.@spawn Cℓ = sum_in_ℓbins(w .* CLobs)
        Threads.@spawn ℓ  = sum_in_ℓbins(w .* L)
        if err_estimate
            Threads.@spawn N   = sum_in_ℓbins(one.(w)) / 2
            Threads.@spawn Cℓ² = sum_in_ℓbins(w .* CLobs.^2)
        end
    end

    if err_estimate
        σℓ  = sqrt.((Cℓ² ./ A .- Cℓ.^2) ./ N)
        InterpolatedCℓs(ℓ./A,  Cℓ./A .± σℓ)
    else
        InterpolatedCℓs(ℓ./A,  Cℓ./A)
    end
end

"""
    ud_grade(f::Field, θnew, mode=:map, deconv_pixwin=true, anti_aliasing=true)

Up- or down-grades field `f` to new resolution `θnew` (only in integer steps).
Two modes are available specified by the `mode` argument:

* `:map`     — Up/downgrade by replicating/averaging pixels in map-space
* `:fourier` — Up/downgrade by extending/truncating the Fourier grid

For `:map` mode, two additional options are possible. If `deconv_pixwin` is
true, deconvolves the pixel window function from the downgraded map so the
spectrum of the new and old maps are the same. If `anti_aliasing` is true,
filters out frequencies above Nyquist prior to down-sampling.

"""
function ud_grade(f::FlatS0{P,T,M}, θnew; mode=:map, deconv_pixwin=(mode==:map), anti_aliasing=(mode==:map)) where {T,M,θ,N,∂mode,P<:Flat{N,θ,∂mode}}
    θnew==θ && return f
    (mode in [:map,:fourier]) || throw(ArgumentError("Available modes: [:map,:fourier]"))

    fac = θnew > θ ? θnew÷θ : θ÷θnew
    (round(Int, fac) ≈ fac) || throw(ArgumentError("Can only ud_grade in integer steps"))
    fac = round(Int, fac)
    Nnew = @. round(Int, N * θ ÷ θnew)
    Pnew = Flat(Nside=Nnew, θpix=θnew, ∂mode=∂mode)
    @unpack Δx,kx,ky,Nx,Ny,nyq = fieldinfo(Pnew,T,M)

    if deconv_pixwin
        PWF = @. T((pixwin(θnew,ky[1:end÷2+1])*pixwin(θnew,kx)')/(pixwin(θ,ky[1:end÷2+1])*pixwin(θ,kx)'))
    else
        PWF = 1
    end

    if θnew>θ
        # downgrade
        if anti_aliasing
            AA = Diagonal(FlatFourier{P}(
                ifelse.((abs.(fieldinfo(P,T,M).ky[1:end÷2+1]) .> nyq) .| (abs.(fieldinfo(P,T,M).kx') .> nyq), 0, 1)
            ))
        else
            AA = 1
        end
        if mode==:map
            fnew = FlatMap{Pnew}(mapslices(mean,reshape((AA*f)[:Ix],(fac,Ny,fac,Nx)),dims=(1,3))[1,:,1,:])
            deconv_pixwin ? FlatFourier{Pnew}(fnew[:Il] ./ PWF) : fnew
        else
            @assert fieldinfo(f).Nside isa Int "Downgrading resolution with `mode=:fourier` only implemented for maps where `Nside isa Int`"
            FlatFourier{Pnew}((AA*f)[:Il][1:(Nnew÷2+1), [1:(isodd(Nnew) ? Nnew÷2+1 : Nnew÷2); (end-Nnew÷2+1):end]])
        end
    else
        # upgrade
        @assert fieldinfo(f).Nside isa Int "Upgrading resolution only implemented for maps where `Nside isa Int`"
        if mode==:map
            fnew = FlatMap{Pnew}(permutedims(hvcat(N,(x->fill(x,(fac,fac))).(f[:Ix])...)))
            deconv_pixwin ? FlatFourier{Pnew}(fnew[:Il] .* PWF' .* PWF[1:Nnew÷2+1]) : fnew
        else
            fnew = FlatFourier{Pnew}(zeros(Nnew÷2+1,Nnew))
            setindex!.(Ref(fnew.Il), f[:Il], 1:(N÷2+1), [findfirst(fieldinfo(fnew).k .≈ fieldinfo(f).k[i]) for i=1:N]')
            deconv_pixwin ? fnew * fac^2 : fnew
        end
    end
end


text:
using ForwardDiff

"""
Create a isotropic Hooke material matrix C

More information: http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
                  https://en.wikipedia.org/wiki/Hooke's_law
                  http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf
Parameters
----------
    E: Float
        Elastic modulus
    ν: Float
        Poisson constant

Returns
-------
    Array{Float64, (6,6)}
"""
function stiffnessTensor(E, ν)
    a = 1 - ν
    b = 1 - 2*ν
    c = 1 + ν
    multiplier = E / (b * c)
    return Float64[a ν ν 0 0 0;
                   ν a ν 0 0 0;
                   ν ν a 0 0 0;
                   0 0 0 b 0 0;
                   0 0 0 0 b 0;
                   0 0 0 0 0 b].*multiplier
end

# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
function find_root!(f, df, x; max_iter=50, norm_acc=1e-10)
    converged = false
    for i=1:max_iter
        dx = df(x) \ -f(x)
        x += dx
        norm(dx) < norm_acc && (converged = true; break)
    end
    converged || error("no convergence!")
    x
end

type State
    C :: Array{Float64, 2}
    stress_y :: Float64
    stress :: Array{Float64, 1}
    strain :: Array{Float64, 1}
end

"""
Equivalent tensile stress.

More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion
    Section: Reduced von Mises equation for different stress conditions

Parameters
----------
    σ: Array{Float64, 6}
        Stress in Voigt notation

Returns
-------
    Float
"""
function stress_eq(stress)
    stress_ten = [stress[1] stress[6] stress[5];
                  stress[6] stress[2] stress[4];
                  stress[5] stress[4] stress[3]]
    stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3)
    s = vec(stress_dev)
    return sqrt(3/2 * dot(s, s))
end

"""
Von Mises Yield criterion

More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf

Parameters
----------
    σ: Array{Float64, 6}
        Stress in Voigt notation
    k: Float64
        Material constant, Yield limit

Returns
-------
    Float
"""
function vonMisesYield(stress, stress_y)
    stress_eq(stress) - stress_y
end

"""
Function for NLsolve. Inside this function are the equations which we want to find root.
Ψ is the yield function below. Functions defined here:

    dσ - C (dϵ - dλ*dΨ/dσ) = 0
                 σₑ(σ) - k = 0

Parameters
----------
    params: Array{Float64, 7}
        Array containing values from solver
    dϵ: Array{Float64, 6}
        Strain rate vector in Voigt notation
    C: Array{Float64, (6, 6)}
        Material tensor
    k: Float
        Material constant, yield limit
    Δt: Float
        time increment
    σ_begin:Array{Float64, 6}
        Stress vector in Voigt notation

Returns
-------
    Array{Float64, 7}, return values for solver
"""
function vonMisesRoot(params, dstrain, C, stress_y, stress_base)

    # Creating wrapper for gradient
    vm_wrap(stress_) = vonMisesYield(stress_, stress_y)
    dfds = ForwardDiff.gradient(vm_wrap)

    # Stress rate and total strain
    dstress = params[1:6]
    stress_tot = vec(stress_base) + params[1:6]

    # Calculating plastic strain rate
    dstrain_p = params[end] * dfds(stress_tot)

    # Calculating equations
    function_1 = dstress - C * (dstrain - dstrain_p)
    function_2 = vm_wrap(stress_tot)
    [vec(function_1); function_2]
end

"""
Stress for ideal plastic von Mises material model

Parameters
----------
    dϵ: Array{Float64, 6}
        Strain rate vector in Voigt notation
    Δt: Float
        time increment
    σ: Array{Float64, 6}
        Last stress vector in Voigt notation
    C: Array{Float64, (6, 6)}
        Material tensor
    k: Float
        Material constant, yield limit

Returns
-------
    Tuple
    Plastic strain rate dϵᵖ and new stress vector σ
"""
function calculate_stress!(dstrain, mat::State, ::Type{Val{:vonMises}})
    stress = mat.stress
    C = mat.C
    stress_y = mat.stress_y
    # Test stress
    stress_tria = stress + C * dstrain

    # Calculating and checking for yield
    yield = vonMisesYield(stress_tria, stress_y)
    if isless(yield, 0.0)
        mat.stress = vec(stress_tria)
    else
        # Yielding happened
        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
        initial_guess = Float64[vec(stress_tria - stress); 0.1]
        f(stress_)  = vonMisesRoot(stress_, dstrain, C, stress_y, stress)
        df     = ForwardDiff.jacobian(f)

        # Calculating root
        result = nlsolve(not_in_place(f, df), initial_guess).zero
        mat.stress  += result[1:6]
    end
end

function calculate_stress(dstrain, stress, C, stress_y,
     ::Type{Val{:vonMises}},
     ::Type{Val{:ElasticPlasticProblem}})
    # Test stress
    stress_tria = stress + C * dstrain

    # Calculating and checking for yield
    yield = vonMisesYield(stress_tria, stress_y)
    if isless(yield, 0.0)
        # stress[i] = stress_tria[i]
        return 0.0
    else
        # Yielding happened
        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
        x = [vec(stress_tria - stress); 0.0]
        f(stress_)  = vonMisesRoot(stress_, dstrain, C, stress_y, stress)
        df     = ForwardDiff.jacobian(f)

        # Calculating root
        # result = nlsolve(not_in_place(f, df), initial_guess).zero
        max_iter = 10
        converged = false
        for i=1:5
            dx = df(x) \ -f(x)
            x += dx
         #   println(x)
            norm(dx) < 1e-10 && (converged = true; break)
        end
        converged || error("no convergence!")
        # stress[:] += x[1:6]
        return x[end]
    end
end

##################################################################################
# ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION          ----- #
##################################################################################

"""
http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm
"""
function stiffnessTensorPlaneStress(E, ν)
    a = 1 - ν^2
    b = 1 - ν
    multiplier = E / a
    return Float64[1 ν 0;
                   ν 1 0;
                   0 0 b].*multiplier
end

# von mises: plane stress
# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
function stress_eq_plane_stress(stress)
    s1, s2, t12  = stress
    # Calculating principal stresses
    # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm
    se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
    se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2)
    return sqrt(se1^2 -se1*se2 + se2^2)
end

# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
function vonMisesYieldPlaneStress(stress, stress_y)
    stress_eq_plane_stress(stress) - stress_y
end

function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base)

    # Creating wrapper for gradient
    vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
    dfds = ForwardDiff.gradient(vm_wrap)

    # Stress rate and total strain
    dstress = params[1:3]
    stress_tot = vec(stress_base) + params[1:3]

    # Calculating plastic strain rate
    dstrain_p = params[end] * dfds(stress_tot)

    # Calculating equations
    function_1 = dstress - C * (dstrain - dstrain_p)
    function_2 = vm_wrap(stress_tot)
    [vec(function_1); function_2]
end

function calculate_stress(dstrain, stress, C, stress_y,
                          ::Type{Val{:vonMises}},
                          ::Type{Val{:PlaneStressElasticPlasticProblem}})
    # Test stress
    dstress = C * dstrain
    stress_tria = stress + dstress

    # Calculating and checking for yield
    yield = vonMisesYieldPlaneStress(stress_tria, stress_y)
    if isless(yield, 0.0)
        return dstress, zeros(3)
    else
        info("yielded")
        # Yielding happened
        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
        x = [vec(stress_tria - stress); 0.0]
        f(stress_)  = vonMisesRootPlaneStress(stress_, dstrain, C, stress_y, stress)
        df = ForwardDiff.jacobian(f)
        # Calculating root
        results = find_root!(f, df, x)
        dstress = results[1:3]
        stress_tot = stress + dstress
        plastic_multiplier = results[end]
        vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
        dfds = ForwardDiff.gradient(vm_wrap)
        dep = plastic_multiplier * dfds(vec(stress_tot))
        info("II ", stress_tot)
        info(vm_wrap(stress_tot))
        return dstress, dep
    end
end

ext: jl


text:
@testset "AbstractIndex" begin
    float_offset = Index(2.0:11.0)
    static_index = Index(UnitSRange(1, 10))
    mutable_index = Index(UnitMRange(1, 10))

    @test @inferred(valtype(float_offset)) == Int
    @test @inferred(keytype(float_offset)) == Float64
    @test @inferred(isempty(float_offset)) == false
    @test @inferred(length(float_offset)) == 10
    @test @inferred(step(float_offset)) == 1
    @test @inferred(first(float_offset)) == 1
    @test @inferred(firstindex(float_offset)) == 2.0
    @test @inferred(last(float_offset)) == 10
    @test @inferred(lastindex(float_offset)) == 11.0

    @test @inferred(AbstractIndices.is_fixed(float_offset)) == true
    @test @inferred(AbstractIndices.is_static(static_index)) == true
    @test @inferred(AbstractIndices.is_dynamic(mutable_index)) == true

    AbstractIndices.set_length!(mutable_index, 12)
    @test length(mutable_index) == 12

    #=
    @testset "iterate vs pairs" begin
        i, state = iterate(symbol_index)
        p = IndexPosition(symbol_index)
        @test (i, state) == (p, p)

        i, state = iterate(symbol_index, state)
        p._state += CartesianIndex(1)
        @test (i, state) == (p, p)

        i, state = iterate(symbol_index, state)
        p._state += CartesianIndex(1)
        @test (i, state) == (p, p)

        @test collect(pairs(symbol_index)) == [:one => 1, :two => 2, :three => 3]
    end

    @testset "NamedIndex" begin
        r1 = range(.1, stop = .2, length=5)
        r2 = ["a", "b", "c", "d"]
        r3 = 2:4

        ind1 = asindex(r1, :a)
        ind1_uint = asindex(asindex(r1, OneTo(UInt(5))), :a)
        ind2 = asindex(r2, :b)
        ind3 = asindex(r3, :c)

        @test asindex(ind1, OneTo(5)) == ind1
        @test asindex(ind1, OneTo(UInt(5))) == ind1_uint
        @test asindex(ind1, ind1) == ind1
        @test asindex(r1, ind1) == ind1
        @test asindex(ind1, ind2) == ind1

        ind4 = asindex(r1, ind1)
        @test values(ind4) == OneTo(5)
    end

    =#
    #=
    @testset "asindex" begin
        @test asindex((:one, :two, :three)) == asindex(symbol_index)
        ks = NamedIndex{:a}(symbol_index)
        @test asindex(ks, nothing) = asindex(symbol_index, :a)
    end
    =#
    x = Index{:a}(1:10)
    y = Index(1:10)
    z = Index{:b}(1:10)
    @testset "combine" begin
        @test 0 == @allocated combine_names(x, y)
        @test combine_names(x, y) == :a
        @test 0 == @allocated combine_names(y, x)
        @test combine_names(y, x) == :a

        @test 0 == @allocated combine_names(x, y)
        @test combine_names(x, z) == :a
        @test 0 == @allocated combine_names(y, x)
        @test combine_names(z, x) == :b

        @test combine_values(x, y) == OneTo(10)
        @test combine_keys(x, y) == 1:10
    end
end


text:
__precompile__()
function exec()
    resize!(ARGS, 2)
    ARGS[1] = "../benchmarks/3QP/toy1"
    ARGS[2] = "toy1.out"
    include("qp.jl")
end

@time exec()
@time exec()

max_stars_repo_name: alsam/julia-samples

max_issues_repo_name: alsam/julia-samples

max_forks_repo_name: alsam/julia-samples


text:
"""
Boundary defining a line, rectangle or box.

# Constructors

    function BoundaryFlat(boundX::BoundaryFlatSubtypes=Free(),
        boundY::BoundaryFlatSubtypes=Free(),
        boundZ::BoundaryFlatSubtypes=Free())

### Arguments
 - **boundX::BoundaryFlatSubtypes** Bound type of the first dimension. (Default Free())
 - **boundY::BoundaryFlatSubtypes** Bound type of the second dimension. (Default Free())
 - **boundZ::BoundaryFlatSubtypes** Bound type of the third dimension. (Default Free())
"""
struct BoundaryFlat<:Boundary

    dims::Int
    boundaries::Array{BoundaryFlatSubtypes,1}
    addSymbols::Array{Symbol}

end

function BoundaryFlat(n::Int=0,boundX::BoundaryFlatSubtypes=Free(),
                        boundY::BoundaryFlatSubtypes=Free(),
                        boundZ::BoundaryFlatSubtypes=Free())

    bound = Array{BoundaryFlatSubtypes,1}([])
    symbols = Array{Symbol,1}([])
    for (pos,i) in enumerate([boundX,boundY,boundZ][1:n])
        push!(bound,i)
        for j in keys(i.addSymbols)
            append!(symbols,i.addSymbols[j])
        end

        if :fMin in fieldnames(typeof(i))
            if n == 1 && !(hasmethod(i.fMin,(FLOAT,FLOAT)))
                error("Function ",i.fMin," must take two arguments (x,t). No function with such signature has been declared.")
            elseif n == 2 && !(hasmethod(i.fMin,(FLOAT,FLOAT,FLOAT)))
                error("Function ",i.fMin," must take three arguments (x,y,t). No function with such signature has been declared.")
            elseif n == 3 && !(hasmethod(i.fMin,(FLOAT,FLOAT,FLOAT,FLOAT)))
                error("Function ",i.fMin," must take four arguments (x,y,z,t). No function with such signature has been declared.")
            end
        end

        if typeof(boundX) != Free
            push!(symbols,POSITIONSYMBOLS[pos])
            unique!(symbols)
        end
    end

    return BoundaryFlat(n,bound,symbols)
end

"""
Boundary defining periodic boundary conditions.

# Constructors

    function Free(;medium::BoundaryMedium=DirichletBoundaryCondition())

### Keyword arguments
- **additional::Array{Symbol,1}** (Default Array{Symbol,1}()) Additional symbols to be updated when the symbol `s` crosses the boundary. (Empty)
- **medium::BoundaryMedium** Boundary method of the medium. (Default DirichletBoundaryCondition())

# Example: A free particle moving in the hole space.

```
bound = Free()
```
"""
struct Free<:BoundaryFlatSubtypes
    addSymbols::Dict{String,Array{Symbol,1}}
    medium::BoundaryMedium
end

function Free(;medium::BoundaryMedium=DirichletBoundaryCondition())
    return Free(Dict{String,Array{Symbol,1}}(),medium)
end

"""
Boundary defining periodic boundary conditions.

# Constructors

    function Periodic(;additional::Array{Symbol,1}=Array{Symbol,1}(),medium=MediumPeriodic())

### Keyword arguments
 - **additional::Array{Symbol,1}** (Default Array{Symbol,1}()) Additional symbols to be updated when the symbol `s` crosses the boundary.
 - **medium::BoundaryMedium** Boundary method of the medium. (Default PeriodicBoundaryCondition())

# Example: An elongated particle with two degrees of freedom `x` and `x1` moving in periodic conditions.

```
bound = Periodic(:x,0,10,additional=[:x1])
```
"""
struct Periodic<:BoundaryFlatSubtypes
    addSymbols::Dict{String,Array{Symbol,1}}
    medium::BoundaryMedium
end

function Periodic(;additional::Array{Symbol,1}=Array{Symbol,1}(),medium::BoundaryMedium=PeriodicBoundaryCondition())
    return Periodic(Dict("add"=>additional),medium)
end

"""
Default boundary defining a non-periodic, bounded space. The parameters specified in keyword arguments can have three behaviours:

 - **stop**: If the symbol `s` crosses the boundary, the specified parameter will be set to zero.
 - **bounce**: If the symbol `s` crosses the boundary, the specified parameter will be set to the position with respect to the boundary as it have bounced back.
 - **reflect**: If the symbol `s` crosses the boundary, the specified parameter will change sign.

All the behaviours can be specified for the minimum, the maximum or both boundaries.

# Constructors

    function Bounded(;
        stop::Array{Symbol,1}=Array{Symbol,1}(),
        stopMin::Array{Symbol,1}=Array{Symbol,1}(),
        stopMax::Array{Symbol,1}=Array{Symbol,1}(),
        bounce::Array{Symbol,1}=Array{Symbol,1}(),
        bounceMin::Array{Symbol,1}=Array{Symbol,1}(),
        bounceMax::Array{Symbol,1}=Array{Symbol,1}(),
        reflect::Array{Symbol,1}=Array{Symbol,1}(),
        reflectMin::Array{Symbol,1}=Array{Symbol,1}(),
        reflectMax::Array{Symbol,1}=Array{Symbol,1}(),
        medium::BoundaryMedium=DirichletBoundaryCondition()
        )

### Keyword arguments
 - **stop::Array{Symbol,1}** (Default Array{Symbol,1}()))
 - **stopMin::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **stopMax::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **bounce::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **bounceMin::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **bounceMax::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **reflect::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **reflectMin::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **reflectMax::Array{Symbol,1}** (Default Array{Symbol,1}())
 - **medium::BoundaryMedium** Boundary method of the medium. (Default DirichletBoundaryCondition())

Example: A particle with some velocity that bounces when touches a border and changes direction.

```
boundX = Bounded(:x,0,10,
            bounce=[:x],
            reflect=[:vx])
boundY = Bounded(:y,0,10,
            bounce=[:y],
            reflect=[:vy])
```
"""
struct Bounded<:BoundaryFlatSubtypes
    addSymbols::Dict{String,Array{Symbol,1}}
    medium::BoundaryMedium
end

function Bounded(;
                stop::Array{Symbol,1}=Array{Symbol,1}(),
                stopMin::Array{Symbol,1}=Array{Symbol,1}(),
                stopMax::Array{Symbol,1}=Array{Symbol,1}(),
                bounce::Array{Symbol,1}=Array{Symbol,1}(),
                bounceMin::Array{Symbol,1}=Array{Symbol,1}(),
                bounceMax::Array{Symbol,1}=Array{Symbol,1}(),
                reflect::Array{Symbol,1}=Array{Symbol,1}(),
                reflectMin::Array{Symbol,1}=Array{Symbol,1}(),
                reflectMax::Array{Symbol,1}=Array{Symbol,1}(),
                medium::BoundaryMedium=DirichletBoundaryCondition())

        d = Dict{String,Array{Symbol,1}}()

        d["stop"]=stop
        d["stopMin"]=stopMin
        d["stopMax"]=stopMax
        d["bounce"]=bounce
        d["bounceMin"]=bounceMin
        d["bounceMax"]=bounceMax
        d["reflect"]=reflect
        d["reflectMin"]=reflectMin
        d["reflectMax"]=reflectMax

    return Bounded(d,medium)
end


text:
#ConstGate
I = ConstGate(I_kernel)
X = ConstGate(X_kernel)
Y = ConstGate(Y_kernel)
Z = ConstGate(Z_kernel)
H = ConstGate(H_kernel)
S = ConstGate(S_kernel)
SDAG = ConstGate(SDAG_kernel)
T = ConstGate(T_kernel)
TDAG = ConstGate(TDAG_kernel)
SWAP = ConstGate(SWAP_kernel)
CNOT = ConstGate(CNOT_kernel)
CX = ConstGate(CNOT_kernel)
TOFFOLI = ConstGate(TOFFOLI_kernel)
FREDKIN = ConstGate(FREDKIN_kernel)
CZ = ConstGate(CZ_kernel)


text:
__precompile__()

module Ripemd
using Compat

export ripemd160, update!, digest!, RIPEMD160_CTX

include("consts.jl")
include("types.jl")
include("interface.jl")
include("transform.jl")

if VERSION < v"0.7.0-DEV.3213"
    codeunits(x) = convert(Array{UInt8}, x)
end

function ripemd160(data::T) where T <: Union{DenseArray{UInt8, 1},
                                             NTuple{N, UInt8} where N}
    ctx = RIPEMD160_CTX()
    update!(ctx, data)
    return digest!(ctx)
end

ripemd160(str::AbstractString) = ripemd160(Vector{UInt8}(codeunits(str)))

end # module Ripemd


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 43e8fa28-fda8-11ea-0051-dd3db29697b8
using Pkg, DrWatson

# ╔═╡ 43e9386e-fda8-11ea-245e-57c3be7b7977
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ e0197532-fda6-11ea-2023-27973dec779c
md"## Clip-05-37s.jl"

# ╔═╡ 43e9b9ea-fda8-11ea-06af-87a022dcae90
md"### snippet 5.29"

# ╔═╡ 43f71298-fda8-11ea-31ae-1136eb790910
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
	df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
	df.neocortex_perc = parse.(Float64, df.neocortex_perc)
	df.lmass = log.(df.mass)
	scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])
end;

# ╔═╡ 43f87ef8-fda8-11ea-1764-9b59eb6c3e8b
md"### snippet 5.1"

# ╔═╡ 44037e3e-fda8-11ea-0d38-bdae456569f2
stan5_6 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] K; // Outcome
 vector[N] M; // Predictor
}

parameters {
 real a; // Intercept
 real bM; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         //Priors
  bM ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bM * M;
  K ~ normal(mu , sigma);     // Likelihood
}
";

# ╔═╡ 44041268-fda8-11ea-3918-c1a03b0e6936
md"## Define the SampleModel, etc."

# ╔═╡ 440ddc4e-fda8-11ea-218b-e9d90d16c8d6
begin
	m5_6s = SampleModel("m5.6", stan5_6);
	m5_6_data = Dict("N" => size(df, 1), "M" => df.lmass_s, "K" => df.kcal_per_g_s);
	rc5_6s = stan_sample(m5_6s, data=m5_6_data);
end;

# ╔═╡ 4414afce-fda8-11ea-04bc-c743afb894bc
if success(rc5_6s)

  # Describe the draws

  post5_6s_df = read_samples(m5_6s, :dataframe)

  title = "Kcal_per_g vs. log mass" * "\nshowing 89% predicted and hpd range"
  plotbounds(
    df, :lmass, :kcal_per_g,
    post5_6s_df, [:a, :bM, :sigma];
    title=title
  )
end

# ╔═╡ 4415fff0-fda8-11ea-3a82-332931a21614
md"## End of clip-05-37s.jl"

# ╔═╡ Cell order:
# ╟─e0197532-fda6-11ea-2023-27973dec779c
# ╠═43e8fa28-fda8-11ea-0051-dd3db29697b8
# ╠═43e9386e-fda8-11ea-245e-57c3be7b7977
# ╟─43e9b9ea-fda8-11ea-06af-87a022dcae90
# ╠═43f71298-fda8-11ea-31ae-1136eb790910
# ╟─43f87ef8-fda8-11ea-1764-9b59eb6c3e8b
# ╠═44037e3e-fda8-11ea-0d38-bdae456569f2
# ╟─44041268-fda8-11ea-3918-c1a03b0e6936
# ╠═440ddc4e-fda8-11ea-218b-e9d90d16c8d6
# ╠═4414afce-fda8-11ea-04bc-c743afb894bc
# ╟─4415fff0-fda8-11ea-3a82-332931a21614


text:
####
#### Defines list of Application Programming Interface (APIs)
####

apis = Any[
    "Home" => "APIs/index.md",
    "Atmos" => Any[
        "AtmosModel" => "APIs/Atmos/AtmosModel.md"
        "Microphysics" => "APIs/Atmos/Microphysics.md"
    ],
    "Ocean" => Any[],
    "Land" => Any[],
    "Common" => Any["Surface Fluxes" => "APIs/Common/SurfaceFluxes.md"],
    "Arrays" => "APIs/Arrays/Arrays.md",
    "Diagnostics" => Any[
        "List of variables" => "APIs/Diagnostics/Diagnostics.md"
        "Input/Output" => "APIs/Diagnostics/InputOutput.md"
    ],
    "Numerics" => Any[
        "Meshes" => "APIs/Numerics/Meshes/Mesh.md",
        "LinearSolvers" => "APIs/Numerics/LinearSolvers/LinearSolvers.md",
        "ODESolvers" => "APIs/Numerics/ODESolvers/ODESolvers.md",
        "Balance Law" => "APIs/Numerics/DGmethods/BalanceLawOverview.md",
    ],
]


text:
# HyperbolicModels.jl

"""
Applies the Hyperbolic model proposed by Rodriguez et al.
"""
function friction(mdl::HyperbolicModel, nrm::Number, vel::Number)
    mdl.friction_coefficient * nrm * mdl.normalization_coefficient *
        sign(mdl.dissipation_coefficient * vel) *
        tanh(abs(mdl.dissipation_coefficient * vel))^(2 * mdl.hysteresis_coefficient - 1) /
        (1.0 + atan(abs(mdl.dissipation_coefficient * vel))^(2 * mdl.stribeck_velocity)) +
        mdl.viscous_damping * vel
end

# ------------------------------------------------------------------------------
# Required routines to support friction model fitting.

function model_from_array(mdl::HyperbolicModel, x::Array{T}) where T <: Number
    HyperbolicModel(x[1], x[2], x[3], x[4], x[5], x[6])
end

function model_to_array(x::HyperbolicModel)
    [
        x.friction_coefficient,
        x.normalization_coefficient,
        x.dissipation_coefficient,
        x.hysteresis_coefficient,
        x.stribeck_velocity,
        x.viscous_damping
    ]
end


text:
@generated function evaluate(itp::AbstractInterpolation, point::SVector{N}) where N
    Expr(:ref, :itp, [:(point[$i]) for i in 1:N]...)
end

function evaluate(itp::AbstractInterpolation, point::AbstractArray)
    itp[point...]
end

function evaluate(cell::Cell{D}, point::AbstractArray) where {D <: AbstractInterpolation}
    leaf = findleaf(cell, point)
    evaluate(leaf.data, leaf.boundary, point)
end

function evaluate(interp::AbstractInterpolation, boundary::HyperRectangle, point::AbstractArray)
    coords = (point - boundary.origin) ./ boundary.widths + 1
    evaluate(interp, coords)
end


text:
using StatsModels: hasintercept, omitsintercept
import StatsModels: drop_intercept, implicit_intercept

struct DroppyMod <: StatisticalModel end
drop_intercept(::Type{DroppyMod}) = true

# define structs for testing implicit intercept trait:

# default for StatisticalModel is true
struct DefaultImplicit <: StatisticalModel end

# override default = true for StatisticalModels
struct NoImplicit <: StatisticalModel end
implicit_intercept(::Type{NoImplicit}) = false

# manual override of default = false
struct YesImplicit end
implicit_intercept(::Type{YesImplicit}) = true

@testset "Model traits" begin
    d = (y = rand(10), x = rand(10), z = [:a, :b, :c])
    sch = schema(d)
    f = @formula(y ~ x)
    f1 = @formula(y ~ 1 + x)
    f0 = @formula(y ~ 0 + x)

    @testset "drop_intercept" begin
        @test_throws ArgumentError apply_schema(f1, sch, DroppyMod)
        ff = apply_schema(f, sch, DroppyMod)
        @test !hasintercept(ff)
        @test !omitsintercept(ff)
        ff0 = apply_schema(f0, sch, DroppyMod)
        @test !hasintercept(ff0)
        @test omitsintercept(ff0)

        @test drop_intercept(DroppyMod()) == drop_intercept(DroppyMod)
        # drop_intercept blocks implicit_intercept == true
        @test implicit_intercept(DroppyMod)

        @testset "categorical promotion" begin
            # drop_intercept == true means that model should always ACT like
            # intercept is present even if it's not specified or even ommitted.
            # (pushes intercept term to the FullRank already seen terms list).

            # full dummy coding
            @test width(apply_schema(@formula(y ~ 0 + z), sch, StatisticalModel).rhs) == 3
            # droppy regular coding
            @test width(apply_schema(@formula(y ~ 0 + z), sch, DroppyMod).rhs) == 2
        end
    end

    @testset "implicit_intercept" begin
        @testset "default" begin
            ff, ff0, ff1 = apply_schema.((f, f0, f1), Ref(sch), Any)
            @test !hasintercept(ff)
            @test !hasintercept(ff0)
            @test hasintercept(ff1)
            @test !omitsintercept(ff)
            @test omitsintercept(ff0)
            @test !omitsintercept(ff1)
        end

        @testset "StatisticalModel default" begin
            ff, ff0, ff1 = apply_schema.((f, f0, f1), Ref(sch), DefaultImplicit)
            @test hasintercept(ff)
            @test !hasintercept(ff0)
            @test hasintercept(ff1)
            @test !omitsintercept(ff)
            @test omitsintercept(ff0)
            @test !omitsintercept(ff1)

            @test implicit_intercept(DefaultImplicit()) == implicit_intercept(DefaultImplicit)
        end

        @testset "Override StatisticalModel default" begin
            ff, ff0, ff1 = apply_schema.((f, f0, f1), Ref(sch), NoImplicit)
            @test !hasintercept(ff)
            @test !hasintercept(ff0)
            @test hasintercept(ff1)
            @test !omitsintercept(ff)
            @test omitsintercept(ff0)
            @test !omitsintercept(ff1)

            @test implicit_intercept(NoImplicit()) == implicit_intercept(NoImplicit)
        end

        @testset "Override Any default" begin
            ff, ff0, ff1 = apply_schema.((f, f0, f1), Ref(sch), YesImplicit)
            # broken because traits are not checked during apply_schema for
            # context that is not <:StatisticalModel
            @test_broken hasintercept(ff)
            @test !hasintercept(ff0)
            @test hasintercept(ff1)
            @test !omitsintercept(ff)
            @test omitsintercept(ff0)
            @test !omitsintercept(ff1)

            @test implicit_intercept(YesImplicit()) == implicit_intercept(YesImplicit)
        end

    end
end


text:
function decode(c::Union{Char,String}, encoder::AbstractEncoder)
    return string(encoder.decode[Symbol(c)])
end

"""
    arabic(s::String[, encoder::AbstractEncoder])

Encode the `String` object into Arabic characters. Custom `encoder`
generated from `@transliterator` can be provided, but default is `Buckwalter`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> arabic(verses(crpsdata[114])[1])
"قُلْ أَعُوذُ بِرَبِّ ٱلنَّاسِ"
```
"""
function arabic(s::String)
    trans = Transliterator()
    return arabic(s, trans)
end

function arabic(s::String, encoder::AbstractEncoder)
    words = ""
    for c in s
        if c === ' '
            words *= " "
        else
            words *= decode(c, encoder)
        end
    end
    return words
end

"""
    verses(quran::AbstractQuran; number=false, start_end=true)

Extract the verses of a `AbstractQuran` object.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> verses(crpsdata[1])[7]
"Sira`Ta {l~a*iyna >anoEamota Ealayohimo gayori {lomagoDuwbi Ealayohimo walaA {lD~aA^l~iyna"
julia> verses(crpsdata[113:114], number=true)[1]
"113:(1,5)"
julia> verses(crpsdata[113:114], number=true, start_end=false)[1]
([113], [1, 2, 3, 4, 5])
```
"""
function verses(quran::AbstractQuran; number::Bool=false, start_end::Bool=true)
    try
        endidx = length(rows(quran.data))
        return verses(quran.data, 1, endidx; number=number, start_end=start_end)
    catch
        return select(quran.data, :form)
    end
end

"""
    verses(quran::TanzilData)

Extract the verses of a `TanzilData` object.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> verses(tnzldata)[1]
"بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
```
"""
function verses(quran::TanzilData)
    return select(quran.data, :form)
end

"""
    verses(data::IndexedTable[, a::Int64[, b::Int64]]; number=false, start_end=true)

Extract the verses of a `IndexedTable` object from row `a` to row `b`.
"""
function verses(data::IndexedTable, a::Int64, b::Int64; number::Bool=false, start_end::Bool=true)
    if number
        if start_end
            verse_out = String[]
            verse_nos = ""
        else
            verse_out = Tuple{Array{Int64,1},Array{Int64,1}}[]
            chapter_nos = Int64[]
            verse_nos = Int64[]
        end
    else
        verse_out = String[]; words = ""
    end

    j = 1
    for i in select(data, :form)[a:b]
        if j > 1
            prev_word  = select(data, :word)[j-1]
            next_word  = select(data, :word)[j]
            prev_verse = select(data, :verse)[j-1]
            next_verse = select(data, :verse)[j]
            prev_chapter = select(data, :chapter)[j-1]
            next_chapter = select(data, :chapter)[j]

            if prev_word !== next_word && prev_verse !== next_verse
                if number
                    if prev_chapter !== next_chapter
                        if start_end
                            push!(verse_out, verse_nos * ")")
                            verse_nos = string(next_chapter) * ":(" * string(next_verse)
                        else
                            push!(verse_out, (chapter_nos, verse_nos))
                            chapter_nos = Int64[]
                            verse_nos = Int64[]

                            push!(chapter_nos, next_chapter)
                            push!(verse_nos, next_verse)
                        end
                    else
                        if start_end
                            if occursin(",", verse_nos)
                                verse_nos = replace(verse_nos, "," * string(prev_verse) => "," * string(next_verse))
                            else
                                verse_nos *= "," * string(next_verse)
                            end
                        else
                            push!(verse_nos, next_verse)
                        end
                    end
                else
                    push!(verse_out, words)
                    words = ""
                end
            elseif prev_word !== next_word && prev_verse === next_verse
                if !number
                    words *= " "
                end
            elseif prev_word === next_word && prev_verse !== next_verse
                if number
                    if prev_chapter !== next_chapter
                        if start_end
                            push!(verse_out, verse_nos * ")")
                            verse_nos = string(next_chapter) * ":(" * string(next_verse)
                        else
                            push!(verse_out, (chapter_nos, verse_nos))
                            chapter_nos = Int64[]
                            verse_nos = Int64[]

                            push!(chapter_nos, next_chapter)
                            push!(verse_nos, next_verse)
                        end
                    else
                        if start_end
                            if occursin(",", verse_nos)
                                verse_nos = replace(verse_nos, "," * string(prev_verse) => "," * string(next_verse))
                            else
                                verse_nos *= "," * string(next_verse)
                            end
                        else
                            push!(verse_nos, next_verse)
                        end
                    end
                else
                    push!(verse_out, words)
                    words = ""
                end
            end

            if !number
                words *= i
            end

            if j === b
                if number
                    if start_end
                        push!(verse_out, verse_nos * ")")
                    else
                        push!(verse_out, (chapter_nos, verse_nos))
                    end
                else
                    push!(verse_out, words)
                end
            end
        else
            if number
                if start_end
                    verse_nos = string(select(data, :chapter)[j]) * ":(" * string(select(data, :verse)[j])
                else
                    push!(chapter_nos, select(data, :chapter)[j])
                    push!(verse_nos, select(data, :verse)[j])
                end
            else
                words *= i
            end

            if j === b
                if number
                    if start_end
                        push!(verse_out, verse_nos * ")")
                    else
                        push!(verse_out, (chapter_nos, verse_nos))
                    end
                else
                    push!(verse_out, words)
                end
            end
        end
        j += 1
    end
    return verse_out
end

"""
    words(quran::AbstractQuran)

Extract words of the input quran.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> words(tnzldata[1][7])[1]
"صِرَٰطَ"
```
"""
function words(quran::AbstractQuran)
    return string.(vcat(split.(verses(quran), " ")...))
end

"""
    chapter_name(quran::AbstractQuran, transliterate::Bool=false; lang::Symbol=:arabic)

Extract the chapter name of the input quran, in either `:arabic` (default) or `:english`

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> chapter_name(crpsdata[13][2][1])
"ٱلرَّعْد"
```
"""
function chapter_name(quran::AbstractQuran, transliterate::Bool=false; lang::Symbol=:arabic)
    chapterlabel = ChapterLabel()
    chapterlabel = getproperty(chapterlabel, lang)
    chapterindex = quran isa Chapter ? quran.numbers : quran.chapters
    if transliterate && lang === :arabic
        trans = Transliterator()
        name = ""
        for c in chapterlabel[chapterindex]
            if c == ' '
                name *= " "
            else
                name *= encode(c, trans)
            end
        end
        return name
    elseif transliterate && lang === :english
        @warn "transliterate only applies for :arabic lang."
        return chapterlabel[chapterindex]
    elseif !transliterate && lang === :arabic
        return chapterlabel[chapterindex]
    else
        return chapterlabel[chapterindex]
    end
end


text:
using PlotNested,
      InteractBase,
      Widgets,
      Plots,
      Test

struct TestData{T}
    n::Symbol
    x::Int
    a::Vector
    t::T
end

gr()
testdata = TestData(:parent, 1, [1,2,3], TestData(:child, 2, [4,5,6], [7 8; 9 10]))
@test plottables(testdata) == (:a => [1,2,3], :a => [4,5,6], :t => [7 8; 9 10])
@test plotnames(testdata) == (:a, :a, :t)
@test plotdata(testdata) == ([1,2,3], [4,5,6], [7 8; 9 10])
@test length(plot_selected(testdata, (true, false, true))) == 2
@test typeof(plot_selected(testdata, (true, false, true))) == Array{Plots.Plot{Plots.GRBackend},1}
@test length(plot_all(testdata)) == 3
@test typeof(plot_all(testdata)) == Array{Plots.Plot{Plots.GRBackend},1}
@test typeof(plotchecks(testdata)) == Array{Widgets.Widget{:checkbox,Bool},1}
@test length(plotchecks(testdata)) == 3

p = autoplot(testdata);
@test typeof(p) == Plots.Plot{Plots.GRBackend}


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events_data

using AWS.Compat
using AWS.UUIDs
"""
    BatchPutMessage()

Sends a set of messages to the AWS IoT Events system. Each message payload is transformed into the input you specify (\"inputName\") and ingested into any detectors that monitor that input. If multiple messages are sent, the order in which the messages are processed isn't guaranteed. To guarantee ordering, you must send messages one at a time and wait for a successful response.

# Required Parameters
- `messages`: The list of messages to send. Each message has the following format: '{ \"messageId\": \"string\", \"inputName\": \"string\", \"payload\": \"string\"}'

"""

batch_put_message(messages; aws_config::AWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}("messages"=>messages); aws_config=aws_config)
batch_put_message(messages, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("messages"=>messages), args)); aws_config=aws_config)

"""
    BatchUpdateDetector()

Updates the state, variable values, and timer settings of one or more detectors (instances) of a specified detector model.

# Required Parameters
- `detectors`: The list of detectors (instances) to update, along with the values to update.

"""

batch_update_detector(detectors; aws_config::AWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}("detectors"=>detectors); aws_config=aws_config)
batch_update_detector(detectors, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectors"=>detectors), args)); aws_config=aws_config)

"""
    DescribeDetector()

Returns information about the specified detector (instance).

# Required Parameters
- `detectorModelName`: The name of the detector model whose detectors (instances) you want information about.

# Optional Parameters
- `keyValue`: A filter used to limit results to detectors (instances) created because of the given key ID.
"""

describe_detector(detectorModelName; aws_config::AWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/"; aws_config=aws_config)
describe_detector(detectorModelName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/", args; aws_config=aws_config)

"""
    ListDetectors()

Lists detectors (the instances of a detector model).

# Required Parameters
- `detectorModelName`: The name of the detector model whose detectors (instances) are listed.

# Optional Parameters
- `maxResults`: The maximum number of results to return at one time.
- `nextToken`: The token for the next set of results.
- `stateName`: A filter that limits results to those detectors (instances) in the given state.
"""

list_detectors(detectorModelName; aws_config::AWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)"; aws_config=aws_config)
list_detectors(detectorModelName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)", args; aws_config=aws_config)


text:
# ---
# title: 861. Score After Flipping Matrix
# id: problem861
# author: AquaIndigo
# date: 2020-12-07
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/score-after-flipping-matrix/description/>
# hidden: true
# ---
#
# We have a two dimensional matrix `A` where each value is `0` or `1`.
#
# A move consists of choosing any row or column, and toggling each value in that
# row or column: changing all `0`s to `1`s, and all `1`s to `0`s.
#
# After making any number of moves, every row of this matrix is interpreted as a
# binary number, and the score of the matrix is the sum of these numbers.
#
# Return the highest possible score.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[0,0,1,1],[1,0,1,0],[1,1,0,0]]
#     Output: 39
#     Explanation: Toggled to [[1,1,1,1],[1,0,0,1],[1,1,1,1]].
#     0b1111 + 0b1001 + 0b1111 = 15 + 9 + 15 = 39
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 20`
#   2. `1 <= A[0].length <= 20`
#   3. `A[i][j]` is `0` or `1`.
#
#
## @lc code=start
using LeetCode

function matrix_score(A::Vector{Vector{Int}})::Int
    for vec in A
        vec .⊻= 0 ⊻ vec[1]  ## the first column are all converted to 1
    end
    w, h, res = length(A[1]), length(A), 0
    ## res = h << (w - 1)
    for i in 1:w
        temp_cnt = 0
        for j in 1:h
            temp_cnt += A[j][i]
        end
        res += max(temp_cnt, h - temp_cnt) << (w - i)
    end
    return res
end
## @lc code=end


text:
using Colors
using CoordinateTransformations
using GeometryBasics
using MeshCat
using Rotations
using LinearAlgebra

# second-order cone
function κ_so(z)
    z1 = z[1:end-1]
    z2 = z[end]

    z_proj = zero(z)

    if norm(z1) <= z2
        z_proj = copy(z)
    elseif norm(z1) <= -z2
        z_proj = zero(z)
    else
        a = 0.5 * (1.0 + z2 / norm(z1))
        z_proj[1:end-1] = a * z1
        z_proj[end] = a * norm(z1)
    end

    return z_proj
end

# visualizer
function second_order_cone_viz(; h = 1.0, trans = 0.1)
    # visualizer
    vis = Visualizer()

    # set empty background
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)

    # generate cone
    l = h * 1.0
    w = h * 1.45
    pyramid = Pyramid(Point3(0.0, 0.0, 0.0), l, w)

    n = 500
    for i = 1:n
        setobject!(vis["pyramid$i"], pyramid,
            MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, trans)))
        settransform!(vis["pyramid$i"],
            compose(Translation(0.0, 0.0, l),
                LinearMap(RotX(π) * RotZ(π * i / n))))
    end

    # generate point
    setobject!(vis["point"], Sphere(Point3f0(0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    settransform!(vis["point"], Translation(0.5, 0.5, 0.5))

    return vis
end

function set_point!(vis, p)
    settransform!(vis["point"], Translation(p...))
end

function project_point!(vis, p)
    p_proj = κ_so(p)
    settransform!(vis["point"], Translation(p_proj...))
    return p_proj
end

vis = second_order_cone_viz(h = 2.0)
render(vis)

p = [1.0, 0.0, 1.0]

set_point!(vis, p)
p_proj = project_point!(vis, p)


text:
# This file contains code that helps to query from the registry to determine the format

"""
`unknown(f)` returns true if the format of `f` is unknown.
"""
unknown(@nospecialize(f::Union{Formatted,Type})) = unknown(formatname(f)::Symbol)
unknown(name::Symbol) = name === :UNKNOWN

const unknown_df = DataFormat{:UNKNOWN}

"""
`info(fmt)` returns the magic bytes/extension information for
`fmt`.
"""
info(@nospecialize(f::Union{Formatted,Type})) = info(formatname(f)::Symbol)
info(sym::Symbol) = sym2info[sym]

"`magic(fmt)` returns the magic bytes of format `fmt`"
magic(@nospecialize(fmt::Type)) = magic(formatname(fmt)::Symbol)
magic(sym::Symbol) = info(sym)[1]

"""
`skipmagic(s::Stream)` sets the position of `s` to be just after the magic bytes.
For a plain IO object, you can use `skipmagic(io, fmt)`.
"""
skipmagic(@nospecialize(s::Stream)) = (skipmagic(stream(s), formatname(s)::Symbol); s)
skipmagic(io, @nospecialize(fmt::Type)) = skipmagic(io, formatname(fmt)::Symbol)
function skipmagic(io, sym::Symbol)
    magic, _ = sym2info[sym]
    skipmagic(io, magic)
    nothing
end
skipmagic(io, @nospecialize(magic::Function)) = nothing
skipmagic(io, magic::Vector{UInt8}) = seek(io, length(magic))
function skipmagic(io, magics::Vector{Vector{UInt8}})
    lengths = map(length, magics)
    l1 = lengths[1]
    all(isequal(l1), lengths) && return seek(io, l1) # it doesn't matter what magic bytes get skipped as they all have the same length
    len = getlength(io)
    tmp = read(io, min(len, maximum(lengths)))
    for m in reverse(magics)  # start with the longest since they are most specific
        if magic_equal(m, tmp)
            seek(io, length(m))
            return nothing
        end
    end
    error("tried to skip magic bytes of an IO that does not contain the magic bytes of the format. IO: $io")
end

function magic_equal(magic, buffer)
    length(magic) > length(buffer) && return false
    for (i,elem) in enumerate(magic)
        buffer[i] != elem && return false
    end
    true
end

function getlength(io, pos=position(io))
    seekend(io)
    len = position(io)
    seek(io, pos)
    return len
end

"""
    query(filename; checkfile=true)

Return a `File` object with information about the
format inferred from the file's extension and/or magic bytes.
If `filename` already exists, the file's magic bytes will take priority
unless `checkfile` is false.
"""
function query(filename; checkfile::Bool=true)
    filename = abspath(filename)
    sym = querysym(filename; checkfile=checkfile)
    return File{DataFormat{sym}}(filename)
end
query(@nospecialize(f::Formatted)) = f

# This is recommended for internal use because it returns Symbol (or errors)
function querysym(filename; checkfile::Bool=true)
    hasmagic(@nospecialize(magic)) = !(isa(magic, Vector{UInt8}) && isempty(magic))

    checkfile &= isfile(filename)
    _, ext = splitext(filename)
    if haskey(ext2sym, ext)
        sym = ext2sym[ext]
        if isa(sym, Symbol)               # there's only one format with this extension
            checkfile || return sym       # since we're not checking, we can return it immediately
            magic = sym2info[sym][1]
            hasmagic(magic) || return sym
            return open(filename) do io
                match(io, magic) && return sym
                # if it doesn't match, we prioritize the magic bytes over the guess based on extension
                return querysym_all(io)[1]
            end
        end
        # There are multiple formats consistent with this extension
        syms = sym::Vector{Symbol}
        checkfile || return syms[1]     # with !checkfile we default to the first. TODO?: change to an error?
        return open(filename) do io
            badmagic = false
            for sym in syms
                magic = sym2info[sym][1]
                if !hasmagic(magic)
                    badmagic = true
                    continue
                end
                match(io, magic) && return sym
            end
            badmagic && error("Some formats with extension ", ext, " have no magic bytes; use `File{format\"FMT\"}(filename)` to resolve the ambiguity.")
            return querysym_all(io)[1]
        end
    end
    !checkfile && return :UNKNOWN
    return open(filename) do io
        return querysym_all(io)[1]
    end
end

function match(io, magic::Vector{UInt8})
    len = getlength(io)
    len < length(magic) && return false
    return magic_equal(magic, read(io, length(magic)))
end

function match(io, magics::Vector{Vector{UInt8}})
    lengths = map(length, magics)
    len = getlength(io)
    tmp = read(io, min(len, maximum(lengths)))
    for m in reverse(magics)  # start with the longest since they are most specific
        if magic_equal(m, tmp)
            return true
        end
    end
    return false
end

function match(io, @nospecialize(magic::Function))
    seekstart(io)
    try
        magic(io)
    catch e
        @error("""There was an error in magic function $magic.
                  Please open an issue at FileIO.jl.""", exception=(e, catch_backtrace()))
        false
    end
end

# Returns sym, magic (the latter may be empty if a magic-function matched)
# Upon return the stream position is set to the end of magic.
function querysym_all(io)
    seekstart(io)
    len = getlength(io)
    lengths = map(magic_list) do p
        length(p.first)
    end
    tmp = read(io, min(len, maximum(lengths)))
    for (magic, sym) in reverse(magic_list)
        isempty(magic) && break
        if magic_equal(magic, tmp)
            seek(io, length(magic))
            return sym, magic
        end
    end
    for (magic, sym) in magic_func
        seekstart(io)
        match(io, magic) && return sym, empty_magic
    end
    seekstart(io)
    return :UNKNOWN, empty_magic
end

function querysym(io::IO)
    if seekable(io)
        sym, _ = querysym_all(io)
        seekstart(io)
        return sym
    end
    # When it's not seekable, we can only work our way upwards in length of magic bytes
    # We're essentially counting on the fact that one of them will match, otherwise the stream
    # is corrupted.
    buffer = UInt8[]
    for (magic, sym) in magic_list
        isempty(magic) && continue
        while length(buffer) < length(magic) && !eof(io)
            push!(buffer, read(io, UInt8))
        end
        if magic_equal(magic, buffer)
            return sym
        end
        eof(io) && break
    end
    return :UNKNOWN
end

"""
`query(io, [filename])` returns a `Stream` object with information about the
format inferred from the magic bytes.
"""
function query(io::IO, filename = nothing)
    sym = querysym(io)
    return Stream{DataFormat{sym}}(io, filename)
end
query(io::IO, @nospecialize(file::Formatted)) = Stream{DataFormat{formatname(file)::Symbol}}(io, filename(file))

seekable(io::IOBuffer) = io.seekable
seekable(::IOStream) = true
seekable(::Any) = false


text:
module Maxwell1D
export advance_fields!

function advance_fields!()
end

end


text:
__precompile__()
module Ogg
using FileIO
import Base: show, convert
export load, save, OggDecoder, OggEncoder

const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("libogg not properly installed. Please run Pkg.build(\"Ogg\")")
end

include("types.jl")
include("decoder.jl")
include("encoder.jl")

end # module


text:
using Graphs
using SimpleGraphs

function AlgConnectivity(L::Symmetric;vecs=false)
    # Calculates the algebraic connectivity of graph with laplacian L
    # Input: L a symmetric laplacian matrix
    # Output: value of algebraic connectivity λ_2
    e_vals = eigfact(L,2:2)
    λ = e_vals.values[1]
    ν = e_vals.vectors
    if vecs
        return λ, ν
    else
        return λ
    end
end
function AlgConnectivity(L::Array{Float64,2},vecs=false)
    Ls = Symmetric(L)
    results = AlgConnectivity(Ls,vecs=vecs)
    return results
end
function AlgConnectivity(L::Array{Int64,2},vecs=false)
    Ls = Symmetric(float(L))
    results = AlgConnectivity(Ls,vecs=vecs)
    return results
end
function AlgConnectivity(g::SimpleGraphs.SimpleGraph,sparse::Bool;vecs=false)
    L = laplace(g)
    results = AlgConnectivity(L,vecs)
    return results
end
function AlgConnectivity(g::GenericGraph,sparse::Bool;vecs=false)
    if sparse
        L = laplacian_matrix_sparse(g)
    else
        L = laplacian_matrix(g)
    end
    results = AlgConnectivity(L,vecs)
    return results
end

function AlgConnectivity(L::SparseMatrixCSC;vecs=false)
    # Calculates the algebraic connectivity of graph with laplacian L
    # Input: L a symmetric laplacian matrix
    # Output: value of algebraic connectivity λ_2
    s = svds(L,nsv=size(L)[1]-1)
    if length(s[2]) < size(L)[1]
        # there are some zero eigen values that weren't calculated
        λ = s[2][end]
    else
        λ = s[2][end-1]
    end

    if vecs
        ν = s[1][:,end-1]
        return λ, ν
    else
        return λ
    end
end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
jlcoverage/utils.jl defines utility functions to support code coverage analysis.
"""

# --- Exports

export display_coverage

# --- Imports

# Standard library
using Logging
using Printf

# External packages
using Coverage: Coverage

# --- Public Functions/Methods

"""
    display_coverage(coverage_data::Vector; startpath::AbstractString)

Display coverage results provided in `coverage_data`. File names are displayed relative
to `startpath`. To display absolute paths, set `startpath` to an empty string. By default,
`startpath` is set to the current working directory.
"""
function display_coverage(coverage_data::Vector; startpath::AbstractString=pwd())
    # --- Check arguments

    # Ensure that `startpath` is an absolute path
    if !isempty(startpath)
        startpath = abspath(startpath)
    end

    # --- Generate coverage report

    # Line formats
    header_line_format = "%-35s %15s %10s %10s\n"
    results_line_format = "%-35s %15d %10d %10s\n"
    horizontal_rule = "-"^79

    # Print header line
    println(horizontal_rule)
    printf(header_line_format, "File", "Lines of Code", "Missed", "Coverage")
    println(horizontal_rule)

    # Initialize line counters
    total_lines_of_code = 0
    total_covered_lines_of_code = 0

    # Print coverage for individual files
    for file_coverage in coverage_data
        if !isempty(startpath)
            filename = relpath(file_coverage.filename, startpath)
        else
            filename = file_coverage.filename
        end

        covered_lines_of_code, lines_of_code = Coverage.get_summary(file_coverage)
        missed_lines_of_code = lines_of_code - covered_lines_of_code
        coverage_pct = 100 * covered_lines_of_code / lines_of_code
        coverage_pct_str = isnan(coverage_pct) ? "N/A" : @sprintf "%9.1f%%" coverage_pct

        printf(
            results_line_format,
            filename,
            lines_of_code,
            missed_lines_of_code,
            coverage_pct_str,
        )

        # Increment line counteres
        total_lines_of_code += lines_of_code
        total_covered_lines_of_code += covered_lines_of_code
    end

    # Print coverage summary
    total_missed_lines_of_code = total_lines_of_code - total_covered_lines_of_code
    coverage_pct = 100 * total_covered_lines_of_code / total_lines_of_code
    coverage_pct_str = isnan(coverage_pct) ? "N/A" : @sprintf "%9.1f%%" coverage_pct

    println(horizontal_rule)
    printf(
        results_line_format,
        "TOTAL",
        total_lines_of_code,
        total_missed_lines_of_code,
        coverage_pct_str,
    )

    # TODO: add count of tests passed, skipped, failed.
    # TODO: add test runtime

    return nothing
end

# --- Private Methods

"""
    printf(fmt, args...)

Print formatted text using format defined by `fmt` and values provided in `args`.
"""
printf(fmt, args...) = @eval @printf($fmt, $(args...))


text:
"""
Threaded replicate simulations. If `nreplicates` is set the data object
will be a vector of replicate data, so we loop over it with threads.
"""
sequencerules!(data::AbstractVector{<:AbstractSimData}, rules) = begin
    Threads.@threads for i in 1:length(data)
        sequencerules!(data[i], rules)
    end
    data
end

"""
Iterate over all rules recursively, swapping source and dest arrays.
Returns the data object with source and dest arrays ready for the next iteration.
"""
sequencerules!(data::SimData) = sequencerules!(data, rules(data))
sequencerules!(data::SimData, rules::Tuple) = begin
    # Run the first rule for the whole frame
    maprule!(data, rules[1])
    # Swap the source and dest arrays
    data = swapsource(data)
    # Run the rest of the rules, recursively
    sequencerules!(data, tail(rules))
end
@inline sequencerules!(data::SimData, rules::Tuple{}) = data

sequencerules!(data::MultiSimData) = begin
    @set! data.data = map(sequencerules!, data.data)
    sequenceinteractions!(data)
end

sequenceinteractions!(multidata::MultiSimData) =
    sequenceinteractions!(multidata, interactions(multidata))
sequenceinteractions!(multidata::MultiSimData, interactions::Tuple) = begin
    # Run the first rule for the whole frame
    mapinteraction!(multidata, interactions[1])
    @set! multidata.data = map(swapsource, data(multidata))

    # Run the rest of the interactions, recursively
    sequenceinteractions!(multidata, tail(interactions))
end
sequenceinteractions!(multidata::MultiSimData, interactions::Tuple{}) = multidata


text:
# # Univariate Solver

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Extension/univariate_solver.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Extension/univariate_solver.ipynb)
# **Contributed by**: Benoît Legat

# When using an SDP solver, the Sum-of-Squares constraint is bridged
# into a semidefinite constraint. This reformulation is done only if the
# solver does not support the Sum-of-Squares constraint. We show in this tutorial how to define a solver that supports such constraint.
# The *same* model with be then solved with or without reformulation depending on the solver.

# We consider a solver that would support finding the SOS decomposition
# of nonnegative univariate polynomials.

using Test #src

module MyUnivariateSolver

import LinearAlgebra
import MathOptInterface
const MOI = MathOptInterface
import MultivariatePolynomials
const MP = MultivariatePolynomials
import SumOfSquares
const SOS = SumOfSquares

function decompose(p::MP.AbstractPolynomial, tol=1e-6)
    vars = MP.effective_variables(p)
    if length(vars) != 1
        error("$p is not univariate")
    end
    x = first(vars)
    lead = MP.leadingcoefficient(p)
    if !isone(lead)
        p = p / lead
    end
    deg = MP.maxdegree(p)
    if isodd(deg)
        return
    end
    d = div(deg, 2)
    companion = zeros(2d, 2d)
    for i in 0:(2d-1)
        if i > 0
            companion[i + 1, i] = 1.0
        end
        companion[i + 1, end] = -MP.coefficient(p, x^i)
    end
    F = LinearAlgebra.schur(complex(companion))
    q = one(p)
    i = 1
    while i <= length(F.values)
        root = F.values[i]
        q *= (x - root)
        if !isapprox(real(root), real(F.values[i+1]), rtol=tol, atol=tol)
            return # Cannot happen for complex conjugate root so it means that we have a root which does not have an even multiplicity This means that the polynomial is not nonnegative
        end
        i += 2
    end
    q1 = MP.mapcoefficientsnz(real, q)
    q2 = MP.mapcoefficientsnz(imag, q)
    return SOS.SOSDecomposition([q1, q2])
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    p::Union{Nothing,MP.AbstractPolynomial}
    decomposition::Union{Nothing,SOS.SOSDecomposition}
    tol::Float64
    function Optimizer()
        return new(nothing, nothing, 1e-6)
    end
end

MOI.is_empty(optimizer::Optimizer) = optimizer.p === nothing
function MOI.empty!(optimizer::Optimizer)
    optimizer.p = nothing
    return
end

function MOI.supports_constraint(::Optimizer, ::Type{<:MOI.VectorAffineFunction}, ::Type{<:SOS.SOSPolynomialSet{SOS.FullSpace}})
    return true
end
function MOI.add_constraint(optimizer::Optimizer, func::MOI.VectorAffineFunction, set::SOS.SOSPolynomialSet{SOS.FullSpace})
    if optimizer.p !== nothing
        error("Only one constraint is supported")
    end
    if !isempty(func.terms)
        error("Only supports constant polynomials")
    end
    optimizer.p = MP.polynomial(func.constants, set.monomials)
    return MOI.ConstraintIndex{typeof(func),typeof(set)}(1) # There will be only ever one constraint so the index does not matter.
end

function MOI.optimize!(optimizer::Optimizer)
    optimizer.decomposition = decompose(optimizer.p, optimizer.tol)
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    if optimizer.decomposition === nothing
        return MOI.INFEASIBLE
    else
        return MOI.OPTIMAL
    end
end

function MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus)
    if optimizer.decomposition === nothing
        return MOI.NO_SOLUTION
    else
        return MOI.FEASIBLE_POINT
    end
end

function MOI.get(optimizer::Optimizer, ::SOS.SOSDecompositionAttribute, ::MOI.ConstraintIndex)
    return optimizer.decomposition
end

end

# We define the same function both for our new solver and SDP solvers.

using SumOfSquares
function decompose(p, solver)
    model = Model(solver)
    con = @constraint(model, p in SOSCone())
    optimize!(model)
    @assert primal_status(model) == MOI.FEASIBLE_POINT
    return sos_decomposition(con)
end

# We consider the following univariate polynomial from
# Example 3.35 of [BPT12].
#
# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

# Using our solver we find the following decomposition in two squares.

using DynamicPolynomials
@polyvar x
p = x^4 + 4x^3 + 6x^2 + 4x + 5
dec = decompose(p, MyUnivariateSolver.Optimizer)

# We can verify as follows that it gives a correct decomposition:

@test polynomial(dec) ≈ p #src
polynomial(dec)

# We can also use a semidefinite solver:

import CSDP
dec = decompose(p, CSDP.Optimizer)
@test length(dec.ps) == 3 #src
@test polynomial(dec) ≈ p #src

# The decomposition is different, it is the sum of 3 squares.
# However, it is also valid:

polynomial(dec)


text:
# Read the image format written by RADMC3D.

# Read the ascii text in `image.out` and parse things into a 3 dimensional matrix (x, y, lambda)

# The first four lines are format information
# iformat # = 1 (2 is local observer)
# im_nx   im_ny #number of pixels in x and y directions
# nlam # number of images at different wavelengths
# pixsize_x  pixsize_y # size of the pixels in cm
# lambda[1] ... lambda[nlam + 1] # wavelengths (um) correspending to images
# pixels, ordered from left to right (increasing x) in the inner loop, and from bottom to top (increasing y) in the outer loop. And wavelength is the outermost loop.

"The image module contains various data types for reading and holding images produced by the radiative transfer programs (via `RADMC-3D`), as well as routines for processing these images."
module image

export imread, imToSky, imToSpec, SkyImage, -
export taureadImg, taureadPos

using ..constants
# import Images # The Images.jl package, not affiliated w/ DiskJockey
# using Dierckx

import Base.- # extend this for Image

# Define an image type, which can store the data as well as pixel spacing

abstract type Image end

"
    RawImage(data, pixsize_x, pixsize_y, lams)

Hold the raw output from `RADMC-3D` in a 3D array (npix_y, npix_x, nlam).
RawImage reflects the RADMC convention that both x and y are increasing
with array index. This means that to display the image as RADMC intends it,
you must set the first array element to the lower left corner."
mutable struct RawImage <: Image
    data::Array{Float64,3} # [ergs/s/cm^2/Hz/ster]
    pixsize_x::Float64 # [cm]
    pixsize_y::Float64 # [cm]
    lams::Vector{Float64} # [μm]
end

"SkyImage is a holder that has both RA and DEC increasing with array index
This convention is necessary for the FFT step
However, to display this image in the traditional sky convention (North up,
East to the left), you must set the first array element to the lower left
corner *and* flip the array along the RA axis: `fliplr(data)` or flipdim(data, 2)"
mutable struct SkyImage <: Image
    data::Array{Float64,3} # [Jy/pixel]
    ra::Vector{Float64} # [arcsec]
    dec::Vector{Float64} # [arcsec]
    lams::Vector{Float64} # [μm]
end

"TausurfImage is designed to hold the results of the RADMC-3D `tausurf` operation. This is
The distance in cm above or below the plane tangent to the observer, which intersects the origin of the model.

From the RADMC3D manual:
The image output file image.out will now contain, for each pixel, the position along the ray in centimeters where τ = τs. The zero point is the surface perpendicular to the direction of observation, going through the pointing position (which is, by default the origin (0, 0, 0)). Positive values mean that the surface is closer to the observer than the plane, while negative values mean that the surface is behind the plane.
So for this datastructure, it's the same thing as RawImage, just instead of intensity, we have distance above/below plane."
mutable struct TausurfImg <: Image
    data::Array{Float64,3} # cm above/behind central projected plane of disk
    pixsize_x::Float64 # [cm]
    pixsize_y::Float64 # [cm]
    lams::Vector{Float64} # [μm]
end

"Encapsulates the 3D position of the pixels representing the tau=1 surface, in the same datashape as the image.
For each pixel, this is the x, y, or z position."
mutable struct TausurfPos
    data_x::Array{Float64,3} # [cm]
    data_y::Array{Float64,3} # [cm]
    data_z::Array{Float64,3} # [cm]
    lams::Vector{Float64} # [μm]
end

"Subtraction for `SkyImage`s"
function -(img1::SkyImage, img2::SkyImage)
    # @assert img1.lams == img2.lams "Images must have the same wavelengths."
    @assert isapprox(img1.ra, img2.ra) "Images must have same RA coordinates."
    @assert isapprox(img1.dec, img2.dec) "Images must have same DEC coordinates."
    data = img1.data - img2.data
    return SkyImage(data, img1.ra, img1.dec, img1.lams)
end

"SkyImage constructor for just a single frame"
SkyImage(data::Matrix{Float64}, ra::Vector{Float64}, dec::Vector{Float64}, lam::Float64) =
SkyImage(reshape(data, tuple(size(data)..., 1)), ra, dec, [lam])

"
    imread(file=\"image.out\")

Read the image file (default=image.out) and return it as an Image object, which contains the fluxes in Jy/pixel,
the sizes and locations of the pixels in arcseconds, and the wavelengths (in microns) corresponding to the images"
function imread(file = "image.out")

    fim = open(file, "r")
    iformat = parse(Int, readline(fim))
    im_nx, im_ny = split(readline(fim))
    im_nx = parse(Int, im_nx)
    im_ny = parse(Int, im_ny)
    nlam = parse(Int, readline(fim))
    pixsize_x, pixsize_y = split(readline(fim))
    pixsize_x = parse(Float64, pixsize_x)
    pixsize_y = parse(Float64, pixsize_y)

    # Read the wavelength array
    lams = Array{Float64}(undef, nlam)
    for i = 1:nlam
        lams[i] = parse(Float64, readline(fim))
    end

    # Create an array with the proper size, and then read the file into it
    data = Array{Float64}(undef, im_ny, im_nx, nlam)

    # According to the RADMC manual, section A.15, the pixels are ordered
    # left to right (increasing x) in the inner loop, and from bottom to top
    # (increasing y) in the outer loop.
    # Basically, pack the array in order but display with origin=lower.

    # Because of the way an image is stored as a matrix, we actually pack the
    # array indices as data[y, x, lam]
    # radmc3dPy achieves something similar by keeping indices the x,y but
    # swaping loop order (radmcPy/image.py:line 675)
    for k = 1:nlam
        readline(fim) # Junk space
        for j = 1:im_ny
            for i = 1:im_nx
                data[j,i,k] = parse(Float64, readline(fim))
            end
        end
    end

    close(fim)

    # According to the RADMC3D manual, the units are *intensity* [erg cm−2 s−1 Hz−1 ster−1]

    return RawImage(data, pixsize_x, pixsize_y, lams)
end

"
    taureadImg(file=\"image_tausurf.out\")

Like imread, but for tausurf. Pixels that have no ``\\tau`` surface are set to `NaN`."
function taureadImg(file = "image_tausurf.out")
    fim = open(file, "r")
    iformat = parse(Int, readline(fim))
    im_nx, im_ny = split(readline(fim))
    im_nx = parse(Int, im_nx)
    im_ny = parse(Int, im_ny)
    nlam = parse(Int, readline(fim))
    pixsize_x, pixsize_y = split(readline(fim))
    pixsize_x = parse(Float64, pixsize_x)
    pixsize_y = parse(Float64, pixsize_y)

    # Read the wavelength array
    lams = Array{Float64}(nlam)
    for i = 1:nlam
        lams[i] = parse(Float64, readline(fim))
    end

    # Create an array with the proper size, and then read the file into it
    data = Array{Float64}(im_ny, im_nx, nlam)

    for k = 1:nlam
        readline(fim) # Junk space
        for j = 1:im_ny
            for i = 1:im_nx
                val = parse(Float64, readline(fim))

                data[j,i,k] = val

                # If RADMC3D signaled that there is no tau=1 surface here, set height to NaN
                if isapprox(val, -1e91)
                    data[j,i,k] = NaN #
                else
                    data[j,i,k] = val
                end
            end
        end
    end

    close(fim)

    # According to the RADMC3D manual, the units are *intensity* [erg cm−2 s−1 Hz−1 ster−1]

    return TausurfImg(data, pixsize_x, pixsize_y, lams)

end

"Read the (x,y,z) positions of the ``\\tau=1`` pixels."
function taureadPos(file = "tausurface_3d.out")
    fim = open(file, "r")
    iformat = parse(Int, readline(fim))
    im_nx, im_ny = split(readline(fim))
    im_nx = parse(Int, im_nx)
    im_ny = parse(Int, im_ny)
    nlam = parse(Int, readline(fim))
    # pixsize_x, pixsize_y = split(readline(fim))
    # pixsize_x = parse(Float64, pixsize_x)
    # pixsize_y = parse(Float64, pixsize_y)

    # Read the wavelength array
    lams = Array{Float64}(nlam)
    for i = 1:nlam
        lams[i] = parse(Float64, readline(fim))
    end

    # Create an array with the proper size, and then read the file into it
    data_x = Array{Float64}(im_ny, im_nx, nlam)
    data_y = Array{Float64}(im_ny, im_nx, nlam)
    data_z = Array{Float64}(im_ny, im_nx, nlam)

    # In contrast to the other image formats, apparently there is only a space before the lams, not inbetween lams.
    readline(fim) # Junk space

    for k = 1:nlam
        for j = 1:im_ny
            for i = 1:im_nx
                val_x, val_y, val_z = split(readline(fim))
                val_x = parse(Float64, val_x)
                val_y = parse(Float64, val_y)
                val_z = parse(Float64, val_z)

                # If RADMC3D signaled that there is no tau=1 surface here (any of the xyz values are -1e91)
                # then set the height to NaN
                if isapprox(val_x, -1e91)
                    data_x[j,i,k] = NaN
                    data_y[j,i,k] = NaN
                    data_z[j,i,k] = NaN
                else
                    data_x[j,i,k] = val_x
                    data_y[j,i,k] = val_y
                    data_z[j,i,k] = val_z
                end
            end
        end
    end

    close(fim)

    return TausurfPos(data_x, data_y, data_z, lams)

end

"
    imToSky(img::RawImage, dpc::Float64)

Convert a `RawImage` to a `SkyImage`. Assumes dpc is parsecs."
function imToSky(img::RawImage, dpc::Float64)

    # The RawImage is oriented with North up and East increasing to the left.
    # this means that for the RawImage, the delta RA array goes from + to -

    # However, the SkyImage actually requires RA (ll) in increasing form.
    # Therefore we flip along the RA axis, fliplr(data) or flipdim(data, 2)

    # println("Min and max intensity ", minimum(img.data), " ", maximum(img.data))
    # println("Pixel size ", img.pixsize_x)
    # println("Steradians subtended by each pixel ",  img.pixsize_x * img.pixsize_y / (dpc * pc)^2)

    # convert from ergs/s/cm^2/Hz/ster to to Jy/ster
    conv = 1e23 # [Jy/ster]

    # Conversion from erg/s/cm^2/Hz/ster to Jy/pixel at 1 pc distance.
    # conv = 1e23 * img.pixsize_x * img.pixsize_y / (dpc * pc)^2

    # Flip across RA dimension
    # dataJy = flipdim(img.data, 2) .* conv
    dataJy = reverse(img.data, dims = 2) .* conv

    (im_ny, im_nx) = size(dataJy)[1:2] # y and x dimensions of the image

    # The locations of pixel centers in cm
    # if n_x = 16, goes [-7.5, -6.5, ..., -0.5, 0.5, ..., 6.5, 7.5] * pixsize
    xx = ((Float64[i for i = 0:im_nx - 1] .+ 0.5) .- im_nx / 2.) * img.pixsize_x
    yy = ((Float64[i for i = 0:im_ny - 1] .+ 0.5) .- im_ny / 2.) * img.pixsize_y

    # The locations of the pixel centers in relative arcseconds
    # Note both RA and DEC increase with array index.
    ra = xx ./ (AU * dpc)
    dec = yy ./ (AU * dpc)

    return SkyImage(dataJy, ra, dec, img.lams)

end

"Take an image and integrate all the frames to create a spatially-integrated spectrum"
function imToSpec(img::SkyImage)

    # pixels in SkyImage are Jy/ster

    # convert from Jy/str to Jy/pixel using str/pixel
    dRA = abs(img.ra[2] - img.ra[1]) * arcsec
    dDEC = abs(img.dec[2] - img.dec[1]) * arcsec

    # Add up all the flux in the pixels to create the spectrum
    flux = dropdims(sum(img.data .* dRA .* dDEC, dims = (1, 2)), dims = (1, 2))
    spec = hcat(img.lams, flux) # First column is wl, second is flux

    return spec
end

"Calculate the integrated line flux"
function integrateSpec(spec::Matrix{Float64}, lam0::Float64)
    # First column is wl, second is flux

    wl = spec[:,1]
    fl = spec[:,2]

    # Convert wl to kms
    vs = c_kms * (spec[:,1] .- lam0) / lam0

    if vs[2] - vs[1] < 0
        reverse!(vs)
        reverse!(fl)
    end

    # println(vs)
    # println(fl)

    spl = Spline1D(vs, fl)
    tot = integrate(spl, vs[1], vs[end])

    return tot
end

"Storage for zeroth moment map"
mutable struct ZerothMoment <: Image
   data::Array{Float64,2} # [Jy · Hz / pixel]
   ra::Vector{Float64} # [arcsec]
   dec::Vector{Float64} # [arcsec]
end

"
    convert(::Type{ZerothMoment}, img::SkyImage)

Convert a `SkyImage` to a `ZerothMoment` map."
function convert(::Type{ZerothMoment}, img::SkyImage)
    # Sum along the frequency axis
    data = squeeze(sum(img.data, (3)), 3)
    return ZerothMoment(data, img.ra, img.dec)
end

end # model


text:
using Latexify
using Markdown
using Test

#inline
@test latexify(:(x * y); env=:inline, cdot=false) == raw"$x y$"

@test latexify(:(x * y); env=:inline, cdot=true) == raw"$x \cdot y$"

@test latexify(:(x*(y+z)*y*(z+a)*(z+b)); env=:inline, cdot=false) ==
raw"$x \left( y + z \right) y \left( z + a \right) \left( z + b \right)$"

@test latexify(:(x*(y+z)*y*(z+a)*(z+b)); env=:inline, cdot=true) ==
raw"$x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right)$"

# raw
@test latexify(:(x * y); env=:raw, cdot=false) == raw"x y"

@test latexify(:(x * y); env=:raw, cdot=true) == raw"x \cdot y"

@test latexify(:(x * (y + z) * y * (z + a) * (z + b)); env=:raw, cdot=false) ==
raw"x \left( y + z \right) y \left( z + a \right) \left( z + b \right)"

@test latexify(:(x * (y + z) * y * (z + a) * (z + b)); env=:raw, cdot=true) ==
raw"x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right)"

# array
@test latexify( [:(x*y), :(x*(y+z)*y*(z+a)*(z+b))]; env=:equation, transpose=true, cdot=false) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
x y & x \left( y + z \right) y \left( z + a \right) \left( z + b \right) \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify( [:(x*y), :(x*(y+z)*y*(z+a)*(z+b))]; env=:equation, transpose=true, cdot=true) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
x \cdot y & x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right) \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

# mdtable
arr = ["x*(y-1)", 1.0, 3*2, :(x-2y), :symb]

@test latexify(arr; env=:mdtable, cdot=false) ==
Markdown.md"| $x \left( y - 1 \right)$ |
| ------------------------:|
|                    $1.0$ |
|                      $6$ |
|                $x - 2 y$ |
|                   $symb$ |
"

@test latexify(arr; env=:mdtable, cdot=true) ==
Markdown.md"| $x \cdot \left( y - 1 \right)$ |
| ------------------------------:|
|                          $1.0$ |
|                            $6$ |
|                $x - 2 \cdot y$ |
|                         $symb$ |
"


text:
struct OptResults
    LB::Float64                         # Lower Bound on the objective
    UB::Float64                         # Upper Bound on the objective
    status::String                      # Status
    solvetime::Float64                  # Solve time in seconds
    iterations::Union{Nothing,Int}      # Number of iterations (for local search)
    nodecount::Union{Nothing,Int}       # Node count of the Branch-and-Bound tree
    lazycount::Union{Nothing,Int}       # Number of lazy constraints added during solving
    num_constraints::Union{Nothing,Int} # Number of constraints in the original model (before solving)

    function OptResults(lb::Float64, ub::Float64, status::String, solvetime::Float64,
                        iterations::Union{Nothing,Int}, nodecount::Union{Nothing,Int},
                        lazycount::Union{Nothing,Int}, num_constraints::Union{Nothing,Int})
        if ~(lb <= ub) && ~(isapprox(lb, ub, rtol=1e-4))
            throw(ArgumentError("Lower Bound value must be less than or equal to Upper Bound. Got values LB=$lb and UB=$ub"))
        end
        if (status == "OPTIMAL") && !(isapprox(lb, ub, rtol=1e-4))
            throw(ArgumentError("An optimal solution implies that the Lower Bound is equal to the Upper Bound. Got values LB=$lb and UB=$ub"))
        end
        if solvetime < 0
            throw(DomainError("Solve time must be a non-negative number (in seconds)."))
        end
        if (!isnothing(iterations)) && (iterations < 0)
            throw(DomainError("Number of iterations must be a positive integer number."))
        end
        if (!isnothing(nodecount)) && (nodecount < 0)
            throw(DomainError("Number of nodes in B&B must be a non-negative integer number."))
        end
        if (!isnothing(lazycount)) && (lazycount < 0)
            throw(DomainError("Number of lazy constraints added must be a non-negative integer number."))
        end
        if (!isnothing(num_constraints)) && (num_constraints < 0)
            throw(DomainError("Number of constraints must be a non-negative integer number."))
        end

        return new(lb, ub, status, solvetime, iterations, nodecount, lazycount, num_constraints)
    end
end

function displayResults(opt_results::OptResults)
    """
    Method to display optimization results.

    Parameters
    ----------
    opt_results : OptResults
        Optimization results object.
    """
    print_string = string("\n",
          "--------------- Opt Results ---------------\n",
          "Obj. LB = ",round(opt_results.LB, digits=5),"\n",
          "Obj. UB = ",round(opt_results.UB, digits=5),"\n",
          "Status: $(opt_results.status)\n",
          "Solve time: ",round(opt_results.solvetime, digits=5)," seconds\n")
    if !isnothing(opt_results.iterations)
        print_string = string(print_string, "Iterations: $(opt_results.iterations)\n")
    end
    if !isnothing(opt_results.nodecount)
        print_string = string(print_string, "Number of nodes: $(opt_results.nodecount)\n")
    end
    if !isnothing(opt_results.num_constraints)
        print_string = string(print_string, "Number of constraints: $(opt_results.num_constraints)\n")
    end
    if !isnothing(opt_results.lazycount)
        print_string = string(print_string, "Number of lazy constraints: $(opt_results.lazycount)\n")
    end
    print(print_string)
end

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ cb162dc0-6713-11eb-2681-ffafd6baf172
using GMT

# ╔═╡ f9bdec00-6720-11eb-326e-51e2f53b642b
md"## Line colors with custom CPT and auto color wrap

In this example the lines color is set using a custom CPT. Pen thickness is assigned automatically.

The custom CPT is used by setting the plot command’s cmap argument to ``true``. This workas because we previously computed the CPT and it will remain in memory until it's consumed when we finish the plot. The **level** argument sets the color to be used from the custom CPT.

In fact, in this case with a CPT already in memory, the **level** option alone would have triggered the line coloring and the **cmap** option could have been droped.
"

# ╔═╡ 4bb76890-6714-11eb-3573-73aa8df83b62
C = makecpt(range=(0,10,1));

# ╔═╡ ff01bd70-6754-11eb-231c-43976dec9ae5
md"Normally we don't need to start a figure with a call to *basemap* because the *plot* function takes care of guessing reasonable  defaults, but in this case we start with a curve with small amplitude and we grow the figure by adding more lines. So if we leave it to automatic guessing one would have to start by the largest amplitude curve."

# ╔═╡ 8a0fcbf0-6714-11eb-12ea-51e244aa662d
basemap(region=(20,30,-10,10), figsize=(12,8))

# ╔═╡ dc6295a0-6713-11eb-330d-0f4f8e3a60e1
x = 20:0.1:30;

# ╔═╡ 8e7f8220-6714-11eb-24fa-c746e71c6f35
for amp=0:10
	y = amp .* sin.(x)
	plot!(x, y, cmap=true, level=amp)
end

# ╔═╡ 3391f530-6720-11eb-2ca9-137fb5325fea
colorbar!(show=true)

# ╔═╡ 7d0fc080-6724-11eb-028e-958fc8284904
md"## Line colors with the automatic color scheme

Here we are showing how to plot several lines at once and color them according to a circular color scheme comprised of 7 distinct colors. We start by generating a dummy matrix 8x5, where rows represent the vertex and the columns hold the lines. To tell the program that first column contains the coordinates and the remaining are all lines to be plotted we use the option **multicol**=*true*
"

# ╔═╡ 5eea7020-6727-11eb-3514-9307a1afbee7
mat = GMT.fakedata(8, 5);

# ╔═╡ b4f7f1b0-6743-11eb-197b-372a05795fdc
lines(mat, multicol=true, name="clines2.png", show=true)

# ╔═╡ 0624c710-6745-11eb-2f33-57f86d0b7c17
md"But if we want choose the colors ourselves, it is also easy though we need to go a bit lower in the data preparation.

The basic data type to transfer tabular data to GMT is the *GMTdataset* and the above command has converted the matrix into a *GMTdataset* under the hood but now we need to create one ourselves and fine control more details, like the colors and line thickness of the individual lines. Not that we have 5 lines but will provide 3 colors and 3 lines thicknesses. When we do this those properties are wrapped modulo its number."

# ╔═╡ c2be1c30-6747-11eb-3650-41893fcf806f
D = mat2ds(mat, color=["brown", "green", "blue"], linethick=[2, 1.0, 0.5, 0.25], multi=true);

# ╔═╡ a6175150-6751-11eb-17dc-b39fcfa29302
md"And now we just call *lines* (but using *plot* would have been the same) with the **D** argument."

# ╔═╡ d7d2cb20-6751-11eb-3e47-d7ac04cfbaac
lines(D, show=true)

# ╔═╡ Cell order:
# ╟─f9bdec00-6720-11eb-326e-51e2f53b642b
# ╠═cb162dc0-6713-11eb-2681-ffafd6baf172
# ╠═4bb76890-6714-11eb-3573-73aa8df83b62
# ╟─ff01bd70-6754-11eb-231c-43976dec9ae5
# ╠═8a0fcbf0-6714-11eb-12ea-51e244aa662d
# ╠═dc6295a0-6713-11eb-330d-0f4f8e3a60e1
# ╠═8e7f8220-6714-11eb-24fa-c746e71c6f35
# ╠═3391f530-6720-11eb-2ca9-137fb5325fea
# ╟─7d0fc080-6724-11eb-028e-958fc8284904
# ╠═5eea7020-6727-11eb-3514-9307a1afbee7
# ╠═b4f7f1b0-6743-11eb-197b-372a05795fdc
# ╟─0624c710-6745-11eb-2f33-57f86d0b7c17
# ╠═c2be1c30-6747-11eb-3650-41893fcf806f
# ╟─a6175150-6751-11eb-17dc-b39fcfa29302
# ╠═d7d2cb20-6751-11eb-3e47-d7ac04cfbaac


text:
"""
    MetaUndirectedGraph

A `MetaGraph` whose underlying graph is of type `Graphs.SimpleGraph`.
"""
const MetaUndirectedGraph = MetaGraph{<:Any,<:Any,<:SimpleGraph}

Graphs.SimpleGraph(g::MetaUndirectedGraph) = g.graph

Graphs.is_directed(::Type{<:MetaUndirectedGraph}) = false
Graphs.is_directed(::MetaUndirectedGraph) = false

function arrange(
    ::MetaUndirectedGraph{<:Any,Label},
    label_1::Label,
    label_2::Label,
    v1::Integer,
    v2::Integer,
) where {T,Label}
    if v1 < v2
        return (label_1, label_2)
    else
        return (label_2, label_1)
    end
end

function arrange(
    g::MetaUndirectedGraph{<:Any,Label}, label_1::Label, label_2::Label
) where {Label}
    v1, v2 = code_for(g, label_1), code_for(g, label_2)
    return arrange(g, label_1, label_2, v1, v2)
end