text:
"""
    cmfunc(match::Symbol)

Returns the full multispectral matching functions based on the CMF data in `cmf_tables.jl`

Possible `match` symbols:

* `:cie31`, CIE1931 2° observer
* `:cie64`, CIE1964 10° observer
* `:cie31_j`, CIE1931 observer with corrections by Judd
* `:cie31_jv`, CIE1964 observer wtih correction by Judd and Vos
* `:cie12_2`, CIE2012 2° observer (current CIE 2° standard observer)
* `:cie12_10`, CIE2012 10° observer (current CIE 10° standard observer)

# Example

```jldoctest
julia> cmfunc(:cie12_2)
CIE12_2(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.003769647, 0.004532416, 0.005446553, 0.006538868, 0.007839699, 0.009382967, 0.01120608, 0.01334965, 0.0158569, 0.01877286  …  2.986206e-6, 2.814999e-6, 2.653663e-6, 2.501725e-6, 2.358723e-6, 2.224206e-6, 2.097737e-6, 1.978894e-6, 1.867268e-6, 1.762465e-6], Real[0.0004146161, 0.0005028333, 0.0006084991, 0.0007344436, 0.0008837389, 0.001059646, 0.001265532, 0.001504753, 0.001780493, 0.002095572  …  1.190771e-6, 1.123031e-6, 1.059151e-6, 9.989507e-7, 9.422514e-7, 8.888804e-7, 8.38669e-7, 7.914539e-7, 7.47077e-7, 7.05386e-7], Real[0.0184726, 0.02221101, 0.02669819, 0.03206937, 0.03847832, 0.04609784, 0.05511953, 0.06575257, 0.07822113, 0.09276013  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function cmfunc(match::Symbol)
    if match in (:cie31, :cie64, :cie12_2, :cie12_10, :cie12_2, :cie12_10)
        if match == :cie31
            CIE31(CMF1931[:, 1], CMF1931[:, 2], CMF1931[:, 3], CMF1931[:, 4])
        elseif match == :cie64
            CIE64(CMF1964[:, 1], CMF1964[:, 2], CMF1964[:, 3], CMF1964[:, 4])
        elseif match == :cie31_j
            CIE31_J(CMF1931_J[:, 1], CMF1931_J[:, 2], CMF1931_J[:, 3], CMF1931_J[:, 4])
        elseif match == :cie31_jv
            CIE31_JV(CMF1931_JV[:, 1], CMF1931_JV[:, 2], CMF1931_JV[:, 3], CMF1931_JV[:, 4])
        elseif match == :cie12_2
            CIE12_2(CMF2012_2[:, 1], CMF2012_2[:, 2],CMF2012_2[:, 3], CMF2012_2[:, 4])
        elseif match == :cie12_10
            CIE12_10(CMF2012_10[:, 1], CMF2012_10[:, 2], CMF2012_10[:, 3], CMF2012_10[:, 4])
        end
    else
        DomainError(match, "Expecting :cie31, :cie64, :cie12_2, or :cie12_10")
    end
end

"""
    cmfunc(conefund::ConeFund)

Returns the full multispectral color matching functions based on the cone fundamental data in `cmf_tables.jl`.
Useful to convert (altered) cone fundamental functions back to color matching functions, e.g. for simulations of color anomalous observers with missing or shifted cone responses.

Possible `conefund` types:

* `LMS31`, CIE1931 2° observer cone fundamentals
* `LMS64`, CIE1964 10° observer cone fundamentals
* `LMS06_2`, CIE2012 2° observer (current CIE 2° standard observer) cone fundamentals
* `LMS06_10`, CIE2012 10° observer (current CIE 10° standard observer) cone fundamentals

# Example

```jldoctest
julia> cmfunc(conefund(:cie12_10))
CIE12_10(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.0001205700118214772, 0.00015063103540534432, 0.00018705601478443535, 0.0002308215774282965, 0.00028292144443875793, [...], 0.010780088746948505, 0.012937395576635658, 0.015480114358540253  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
```
"""
function cmfunc(conefund::ConeFund)
    if typeof(conefund) in (LMS31, LMS64, LMS06_2, LMS06_10)
        if isa(conefund, LMS31)
            ncmf = convmat(:cie31)^-1 * [conefund.l, conefund.m, conefund.s]
            CIE31(conefund.λ, ncmf[1], ncmf[2], ncmf[3])
        elseif isa(conefund, LMS64)
            ncmf = convmat(:cie64)^-1 * [conefund.l, conefund.m, conefund.s]
            CIE64(conefund.λ, ncmf[1], ncmf[2], ncmf[3])
        elseif isa(conefund, LMS06_2)
            ncmf = convmat(:cie12_2)^-1 * [conefund.l, conefund.m, conefund.s]
            CIE12_2(conefund.λ, ncmf[1], ncmf[2], ncmf[3])
        elseif isa(conefund, LMS06_10)
            ncmf = convmat(:cie12_10)^-1 * [conefund.l, conefund.m, conefund.s]
            CIE12_10(conefund.λ, ncmf[1], ncmf[2], ncmf[3])
        end
    else
        DomainError(conefund, "Expecting cone fundamentals of type `LMS31`, `LMS64`, `LMS06_2`, or `LMS06_10`")
    end
end

"""
    conefund(g::Symbol)

Returns the CIE cone fundamentals according to the requested lms function.

Possible symbols for `g`:

* `:lms31`, CIE1931 2° observer
* `:lms64`, CIE1964 10° observer
* `:lms06_2`, CIE2012 2° observer (current CIE 2° standard observer)
* `:lms06_10`, CIE2012 10° observer (current CIE 10° standard observer)

# Example

```jldoctest
julia> conefund(:cie64)
LMS64(Real[360, 361, 362, 363, 364, 365, 366, 367, 368, 369  …  821, 822, 823, 824, 825, 826, 827, 828, 829, 830], Real[1.5474770176000004e-8, 2.3435133003999998e-8, 3.5282213430000004e-8, [...], 1.4321467274333274e-9, 1.384067146347377e-9, 1.3377105658218067e-9, 1.2930114805124651e-9, 1.2499070138616995e-9, 1.208336807203753e-9])
```
"""
function conefund(g::Symbol)
    if g == :cie64
        lms = convmat(g) * [CMF1964[:,2], CMF1964[:,3], CMF1964[:,4]]
        @inbounds for l in 161:471
            lms[3][l] = 10^(10402.1 / (371 + l) - 21.27185)
        end
        return LMS64(collect(360:830),lms[1],lms[2],lms[3])
    elseif g == :cie31
        lms = convmat(g) * [CMF1931[:,2], CMF1931[:,3], CMF1931[:,4]]
        return LMS31(collect(360:830), lms[1], lms[2], lms[3])
    elseif g == :cie12_2
        return LMS06_2(LMS2006_2[:, 1], LMS2006_2[:, 2], LMS2006_2[:, 3], LMS2006_2[:, 4])
    elseif g == :cie12_10
        return LMS06_10(LMS2006_10[:, 1], LMS2006_10[:, 2], LMS2006_10[:, 3], LMS2006_10[:, 4])
    else
        throw(DomainError(g, "Color Matching function does not exist."))
    end
end

"""
    conefund(g::CMatch)

Returns the CIE cone fundamentals associated with a CMF.

Possible types for `g`:

* `CIE31`, CIE1931 2° observer
* `CIE64`, CIE1964 10° observer
* `CIE12_2`, CIE2012 2° observer (current CIE 2° standard observer)
* `CIE12_10`, CIE2012 10° observer (current CIE 10° standard observer)

# Example

```jldoctest
julia> conefund(cmfunc(:cie31))
LMS31(Real[360, 361, 362, 363, 364, 365, 366, 367, 368, 369  …  821, 822, 823, 824, 825, 826, 827, 828, 829, 830], Real[6.209077859e-6, 6.950775712086998e-6, 7.782164998607999e-6, [...], 5.5398906285000006e-8, 5.1646151529000004e-8, 4.814798721600001e-8, 4.488794319300001e-8, 4.1849549307000004e-8])
```
"""
function conefund(g::CMatch)
    if g isa CIE31
        conefund(:cie31)
    elseif g isa CIE64
        conefund(:cie64)
    elseif g isa CIE12_2
        c = convmat(:cie12_2) * [g.x, g.y, g.z]
        LMS06_2(g.λ, c[1], c[2], c[3])
    elseif g isa CIE12_10
        c = convmat(:cie12_10) * [g.x, g.y, g.z]
        LMS06_10(g.λ, c[1], c[2], c[3])
    else
        throw(DomainError(g, "Cone fundamental for this type does not exist."))
    end
end

"""
    convmat(g::Symbol)

Returns the conversion matrix to convert CMFs to their associated cone fundamentals (LMS).

Result: 3 × 3 matrix

# Example

`convmat(:cie12_10)`

```3×3 Matrix{Float64}:
 1.93986   -1.34664   0.430449
 0.692839   0.349676  0.0
 0.0        0.0       2.14688```
"""
function convmat(g::Symbol)
    if g == :cie12_2
               [1.94735469 -1.41445123 0.36476327;
                0.68990272  0.34832189 0.00000000;
                0.00000000  0.00000000 1.93485343]
    elseif g == :cie12_10
               [1.93986443 -1.34664359 0.43044935;
                0.69283932  0.34967567 0.00000000;
                0.00000000  0.00000000 2.14687945]
    elseif g == :cie64
               [0.236157  0.826427  -0.045710;
               -0.431117  1.206922   0.090020;
                0.040557 -0.019683   0.48195]
    elseif g == :cie31
               [0.236157  0.826427  -0.045710;
               -0.431117  1.206922   0.090020;
                0.040557 -0.019683   0.48195]
    else throw(DomainError(g, "Color matching function can’t be converted"))
    end
end

function invconvmat(g::Symbol)
    g == :lms06_2 ? convmat(:cie12_2)^-1 :
    g == :lms06_10 ? convmat(:cie12_10)^-1 :
    g == :lms64 ? convmat(:cie64)^-1 :
    g == :lms31 ? convmat(:cie31)^-1 : throw(DomainError(g, "Cone fundamentals don’t exist."))
end

"""
    colorvisiondeficiency(deftype::Symbol, Δλ::Real = 0.0)

Returns the cone fundamentals for a color deficient observer of a given type.

Allowed types:

* `:protanope`
* `:deuteranope`
* `:tritanope`
* `:protanomalous`
* `:deuteranomalous`
* `:tritanomalous`
"""
function colorvisiondeficiency(model::Symbol, deftype::Symbol, shift::Real = 0.0)
    # color primaries SPDs of RBB color space

    # transform RGB triples to opponent colors
    WSR = ρWS * sum(φR .* WS .* SPECENV.Δλ)
    WSG = ρWS * sum(φG .* WS .* SPECENV.Δλ)
    WSB = ρWS * sum(φB .* WS .* SPECENV.Δλ)
    YBR = ρYB * sum(φR .* YB .* SPECENV.Δλ)
    YBG = ρYB * sum(φG .* YB .* SPECENV.Δλ)
    YBB = ρYB * sum(φB .* YB .* SPECENV.Δλ)
    RGR = ρRG * sum(φR .* RG .* SPECENV.Δλ)
    RGG = ρRG * sum(φG .* RG .* SPECENV.Δλ)
    RGB = ρRG * sum(φB .* RG .* SPECENV.Δλ)
    # map RGB to opponent-color space for anomalous observer
    Γ = [WSR WSG WSB;
         YBR YBG YBB;
         RGR RGG RGB]
    # map RGB to opponent-color space for normal observer
    Γn= []
end

"""
    colorvisiondeficiency(model::Symbol, deftype::Symbol, c::Spectrum, illum::LSpec = D65_illuminant() , shift::Real = 0.0)

color vision deficiency simulation using Jozef Cohen’s matrix-R method.
"""
function colorvisiondeficiency(model::Symbol, deftype::Symbol, c::Spectrum, illum::LSpec = D65_illuminant() , shift::Real = 0.0)
    cflms = conefund(SPECENV.cmf)
    if deftype == :protanopia
        Acvd = conefund(match)
        Acvd.l = zeros(length(Acvd.l))
    elseif deftype == :deuteranopia
        Acvd = conefund(match)
        Acvd.m = zeros(length(Acvd.m))
    elseif deftype == :tritanopia
        Acvd = conefund(match)
        Acvd.s = zeros(length(Acvd.s))
    elseif deftype == :protanomaly
        Acvd = shift_colormatch(conefund(match),:l,shift)
    elseif deftype == :deuteranomaly
        Acvd = shift_colormatch(conefund(match),:m,shift)
    elseif deftype == :tritanomaly
        Acvd = shift_colormatch(conefund(match),:s,shift)
    else
        DomainError(deftype, "Color anomaly type does not exist.")
    end

    Clms = fundamental_metamer(illum, c, cflms)
    Ccvd = fundamental_metamer(illum, c, Acvd)
    # B = metameric_black(illum, c, cf)
    # Rlms = matrix_R(cf)
    # Rcvd = matrix_R(Acvd)
    diff = tristimulus(illum::LSpec, rad::Spectrum, match = conefund(SPECENV.cmf))
    env = SPECENV
    λc = collect(env.λmin:env.Δλ:env.λmax)
    Wee = luminance_spec(λc, ones(Float64,length(λc)))
    eeClms = matrix_R(cflms)
    Alms = matrix_A(cflms)
    Pinv = Alms * (Alms' * Alms)^-1
    ycrcb = mlms_ycc() * Alms * eeClms

end

cvd = colorvisiondeficiency

"""
    vybrg(conefund::ConeFund = conefund(SPECENV.cmf))

Returns the Ingling-Tsou opponent color functions based on the cone fundamental functios.

# Example

```jldoctest
julia> c = conefund(:cie12_2); vybrg(c)
VYBRG(Real[390.0, 391.0, 392.0, 393.0, 394.0, 395.0, 396.0, 397.0, 398.0, 399.0  …  821.0, 822.0, 823.0, 824.0, 825.0, 826.0, 827.0, 828.0, 829.0, 830.0], Real[0.0003963417152, 0.0004807961832, 0.0005820062634, [...], 1.4450193672e-6, 1.3620748624e-6, 1.2840832415999998e-6,] 1.2107868407999998e-6, 1.1419364999999999e-6, 1.0772919015999998e-6, 1.0166218466399999e-6])
```
"""
function vybrg(conefund::ConeFund = conefund(SPECENV.cmf))
    inglingtsou_opponent(conefund)
end


text:
using TargetedLearning.LReg
using StatsFuns

facts("Testing LReg") do

    n, p, newn = 50, 3, 20
    srand(10)
    x = rand(n, p)
    y = round(rand(n))
    newx = rand(newn, p)

    context("without offest or weights") do
        lrfit = lreg(x, y)
        @fact lrfit --> is_a(LR)
        @fact predict(lrfit, newx) --> roughly(logistic(newx * lrfit.β))

        @fact_throws predict(lrfit, newx, offset=rand(newn))

        pred = logistic(x * lrfit.β)
        total_loss = 0.0
        lp=linpred(lrfit, x)
        for i in 1:length(y) #NumericExtensions
          total_loss += LReg.loss(y[i], lp[i])
        end
        @fact total_loss/length(y) --> roughly(mean(- y .* log(pred) .- (1-y) .* log(1 .- pred)))
    end

    context("with offset") do
        off = rand(n)
        newoff = rand(newn)
        lrfitoff = lreg(x, y, offset=off)
        @fact predict(lrfitoff, newx, offset=newoff) -->  roughly(logistic(newx * lrfitoff.β .+ newoff))
        @fact_throws predict(lrfitoff, newx)
    end

    context("with weights") do
        lrfitwts = lreg(x, y, wts=ones(50))
        @fact lrfitwts.β --> roughly(lreg(x,y).β)
        rwts = rand(50)
        lrfitwts2 = lreg(x, y, wts=rwts)
        @fact lrfitwts2.β --> not(roughly(lrfitwts.β))
    end

    context("fitting on a subset of columns") do
        sslrfit = lreg(x, y, subset=[3,1])
        @fact sslrfit.β[[1,3]] --> roughly(lreg(x[:, [1,3]], y).β)
        @fact sslrfit.β[2] --> 0.0
        @fact sslrfit.idx --> [1,3]
        @fact sslrfit.β --> lreg(x, y, subset=[1,3]).β
    end

    context("works when x is a vector") do
        lrvec = lreg(rand(10), round(rand(10)))
        @fact lrvec --> is_a(LR)
        @fact predict(lrvec, rand(10)) --> is_a(Vector)
    end

    context("intercept only") do
        @pending "without weights" --> nothing
        @pending "with weights" --> nothing
    end

end

ext: jl


text:
# ---
# title: 1021. Remove Outermost Parentheses
# id: problem1021
# author: Indigo
# date: 2022-03-19
# difficulty: Easy
# categories: Stack
# link: <https://leetcode.com/problems/remove-outermost-parentheses/description/>
# hidden: true
# ---
#
# A valid parentheses string is either empty `("")`, `"(" + A + ")"`, or `A +
# B`, where `A` and `B` are valid parentheses strings, and `+` represents string
# concatenation.  For example, `""`, `"()"`, `"(())()"`, and `"(()(()))"` are
# all valid parentheses strings.
#
# A valid parentheses string `S` is **primitive** if it is nonempty, and there
# does not exist a way to split it into `S = A+B`, with `A` and `B` nonempty
# valid parentheses strings.
#
# Given a valid parentheses string `S`, consider its primitive decomposition: `S
# = P_1 + P_2 + ... + P_k`, where `P_i` are primitive valid parentheses strings.
#
# Return `S` after removing the outermost parentheses of every primitive string
# in the primitive decomposition of `S`.
#
#
#
# **Example 1:**
#
#
#
#     Input: "(()())(())"
#     Output: "()()()"
#     Explanation:
#     The input string is "(()())(())", with primitive decomposition "(()())" + "(())".
#     After removing outer parentheses of each part, this is "()()" + "()" = "()()()".
#
#
# **Example 2:**
#
#
#
#     Input: "(()())(())(()(()))"
#     Output: "()()()()(())"
#     Explanation:
#     The input string is "(()())(())(()(()))", with primitive decomposition "(()())" + "(())" + "(()(()))".
#     After removing outer parentheses of each part, this is "()()" + "()" + "()(())" = "()()()()(())".
#
#
# **Example 3:**
#
#
#
#     Input: "()()"
#     Output: ""
#     Explanation:
#     The input string is "()()", with primitive decomposition "()" + "()".
#     After removing outer parentheses of each part, this is "" + "" = "".
#
#
#
#
# **Note:**
#
#   1. `S.length <= 10000`
#   2. `S[i]` is `"("` or `")"`
#   3. `S` is a valid parentheses string
#
#
#
#
## @lc code=start
using LeetCode

function remove_outer_parentheses(s::String)
    res = ""
    cnt = 0
    for ch in s
        flg = cnt < 0
        cnt += ch == '(' ? -1 : 1
        (flg && cnt < 0) && (res *= ch)
    end
    return res
end
## @lc code=end


text:
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

mutable struct MOISolution{R}
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    dual_status::MOI.ResultStatusCode
    primal::Vector{R}
    dual::Vector{R}
    slack::Vector{R}
    objval::R
    solve_time::R
    # other params
    objective_bound::R
    relative_mip_gap::R
    iteration_count::Int
    barrier_iterations::Int
    node_count::Int
end
MOISolution(R::DataType) = MOISolution{R}(MOI.OPTIMIZE_NOT_CALLED,
                            MOI.NO_SOLUTION,
                            MOI.NO_SOLUTION,
                            R[], R[], R[], NaN, 0.0,
                            #other
                            NaN, NaN, 1, 1, 0)

# Used to build the data with allocate-load during `copy_to`.
# When `optimize!` is called, a the data is passed to SCS
# using `SCS_solve` and the `ModelData` struct is discarded
mutable struct ModelData
    # problem dimensions
    m::Int # Number of rows/constraints
    n::Int # Number of cols/variables
    # A matrix in sparse form
    I::Vector{Int} # List of rows
    J::Vector{Int} # List of cols
    V::Vector{Float64} # List of coefficients
    # constraint bounds
    b::Vector{Float64} # constants
    sense::Vector{ConstraintSense} # constants
    # objective
    objconstant::Float64 # The objective is min c'x + objconstant
    c::Vector{Float64}
    # variable type
    # v_type::Vector{}
    # function ModelData()
    #     new(0,
    #     0,
    #     Int[],
    #     Int[],
    #     Float64[],
    #     Float64[],
    #     ConstraintSense[],
    #     NaN,
    #     Float64[]
    #     )
    # end
end

# This is tied to SCS's internal representation
mutable struct ConeData
    nrows::Dict{Int,Int}
    eq::Int # number of == constraints
    gt::Int # number of >= constraints
    # lt::Int # number of <= equality constraints
    function ConeData()
        new(Dict{Int,Int}(), 0, 0)#, 0)
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    solve_function::Function
    cone::ConeData
    maxsense::Bool
    data::Union{Nothing, ModelData} # only non-Void between MOI.copy_to and MOI.optimize!
    sol::MOISolution
    options
    function Optimizer(solve_function = no_solve_function; options...)
        new(solve_function, ConeData(), false, nothing, MOISolution(Float64), options)
    end
end

function no_solve_function(optimizer::Optimizer)
    println("No solver set.")
    MOISolution(Float64)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Matrix optimizer"

function MOI.is_empty(optimizer::Optimizer)
    !optimizer.maxsense && optimizer.data === nothing
end
function MOI.empty!(optimizer::Optimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    # optimizer.sol.ret_val = 0
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.SingleVariable},
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return true
end

function MOI.supports(::Optimizer,
                      ::Union{MOI.ConstraintPrimalStart,
                              MOI.ConstraintDualStart},
                      ::Type{<:MOI.ConstraintIndex})
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{<:MOI.VectorAffineFunction{Float64}},
    ::Type{<:Union{MOI.Zeros, MOI.Nonnegatives}})
    # ::Type{MOI.Nonnegatives})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

using Compat.SparseArrays

# Computes cone dimensions
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    return ci.value
end
#_allocate_constraint: Allocate indices for the constraint `f`-in-`s`
# using information in `cone` and then update `cone`
function _allocate_constraint(cone::ConeData, f, s::MOI.Zeros)
    ci = cone.f
    cone.eq += MOI.dimension(s)
    return ci
end
function constroffset(cone::ConeData,
                      ci::CI{<:MOI.AbstractFunction, MOI.Nonnegatives})
    return cone.eq + ci.value
end
function _allocate_constraint(cone::ConeData, f, s::MOI.Nonnegatives)
    ci = cone.gt
    cone.gt += MOI.dimension(s)
    return ci
end
function constroffset(optimizer::Optimizer, ci::CI)
    return constroffset(optimizer.cone, ci::CI)
end
function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    return CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
end

output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
# constrrows: Recover the number of rows used by each constraint.
# When, the set is available, simply use MOI.dimension
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
# When only the index is available, use the `optimizer.ncone.nrows` field
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.nrows[constroffset(optimizer, ci)]

constraint_senses(s::MOI.Zeros) = fill(EQUAL_TO, MOI.dimension(s))
constraint_senses(s::MOI.Nonnegatives) = fill(GREATER_THAN, MOI.dimension(s))
function MOIU.load_constraint(optimizer::Optimizer, ci, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    I, J, V = findnz(A)
    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    optimizer.cone.nrows[offset] = length(rows)
    i = offset .+ rows
    b = f.constants
    # Input format is Ax + b >= 0
    # Ouput format is A2x <= b -> -Ax <= b
    # The SCS format is b - Ax ∈ cone
    optimizer.data.b[i] = b
    optimizer.data.sense[i] = constraint_senses(s)
    append!(optimizer.data.I, offset .+ I)
    append!(optimizer.data.J, J)
    append!(optimizer.data.V, -V)
end

function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
    optimizer.cone = ConeData()
    VI.(1:nvars)
end

function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.eq + cone.gt
    I = Int[]
    J = Int[]
    V = Float64[]
    b = zeros(m)
    sense = fill(EQUAL_TO, m)
    c = zeros(nvars)
    optimizer.data = ModelData(m, nvars, I, J, V, b, sense, 0., c)
    # `optimizer.sol` contains the result of the previous optimization.
    # It is used as a warm start if its length is the same, e.g.
    # probably because no variable and/or constraint has been added.
    optimizer.sol.primal = zeros(nvars)
    optimizer.sol.dual = zeros(m)
    optimizer.sol.slack = zeros(m)
    nothing
end

function MOIU.allocate(::Optimizer, ::MOI.VariablePrimalStart,
                       ::MOI.VariableIndex, ::Union{Nothing, Float64})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintPrimalStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(::Optimizer, ::MOI.ConstraintDualStart,
                       ::MOI.ConstraintIndex,
                       ::Union{Nothing, AbstractVector{Float64}})
end
function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MAX_SENSE
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.Union{MOI.SingleVariable,
                                   MOI.ScalarAffineFunction{Float64}})
end

function MOIU.load(::Optimizer, ::MOI.VariablePrimalStart,
                   ::MOI.VariableIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.VariablePrimalStart,
                   vi::MOI.VariableIndex, value::Float64)
    optimizer.sol.primal[vi.value] = value
end
function MOIU.load(::Optimizer, ::MOI.ConstraintPrimalStart,
                   ::MOI.ConstraintIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ConstraintPrimalStart,
                   ci::MOI.ConstraintIndex, value)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    optimizer.sol.slack[offset .+ rows] .= value
end
function MOIU.load(::Optimizer, ::MOI.ConstraintDualStart,
                   ::MOI.ConstraintIndex, ::Nothing)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ConstraintDualStart,
                   ci::MOI.ConstraintIndex, value)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    optimizer.sol.dual[offset .+ rows] .= value
end
function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.SingleVariable)
    MOIU.load(optimizer,
              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
              MOI.ScalarAffineFunction{Float64}(f))
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          optimizer.data.n))
    optimizer.data.objconstant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
    return nothing
end

function MOI.optimize!(optimizer::Optimizer)
    sol = optimizer.solve_function(optimizer)
    # sol.objval = (optimizer.maxsense ? -1 : 1) *
    #                 dot(optimizer.data.c, sol.primal) +
    #                 optimizer.objconstant
    optimizer.sol = sol
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
    return optimizer.sol.solve_time
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    return optimizer.sol.termination_status
end

MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue) = optimizer.sol.objval

function MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus)
    return optimizer.sol.primal_status
end
function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, vi::VI)
    optimizer.sol.primal[vi.value]
end
function MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI})
    return MOI.get.(optimizer, a, vi)
end
function MOI.get(optimizer::Optimizer, ::MOI.ConstraintPrimal,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    primal = optimizer.sol.slack[offset .+ rows]
    return primal
end

function MOI.get(optimizer::Optimizer, ::MOI.DualStatus)
    return optimizer.sol.dual_status
end
function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual,
                 ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    dual = optimizer.sol.dual[offset .+ rows]
    return dual
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1

#=
    query model data
=#

A_matrix_triplets(data::ModelData) = (data.I, data.J, data.V)
A_matrix_sparse(data::ModelData) = sparse(A_matrix_triplets(data)..., data.m, data.n)
A_matrix_dense(data::ModelData) = Matrix(A_matrix_sparse(data))

function LPSolverForm(optimizer)
    data = optimizer.data
    raw_lp = LPSolverForm{Float64}(
        optimizer.maxsense ? MOI.MAX_SENSE : MOI.MIN_SENSE,
        data.c,
        A_matrix_dense(data),
        data.b,
        data.sense,
        fill(-Inf, data.n),
        fill(Inf, data.n)
    )
end


text:
module ViennaRNA

import Base
using Unitful: @u_str, Quantity, uconvert, ustrip
using UnsafePointers: UnsafePtr

# types
export FoldCompound, Pairtable
# wrapper functions
export bp_distance, bpp, centroid, energy, ensemble_defect,
    heat_capacity, inverse_fold, inverse_pf_fold, mea,
    mean_bp_distance, mfe, neighbors, partfn, pbacktrack, plot_coords,
    prob_of_structure, subopt, subopt_zuker, tree_edit_dist
# additional utility functions
export basepairs, plot_structure, prob_of_basepairs

include("../lib/LibRNA.jl")
import .LibRNA

module LibRNA_Helper
function free_subopt_solutions(ptr::Ptr)
    ptr == C_NULL && return
    i = 1
    while true
        sol = unsafe_load(ptr, i)
        sol.structure == C_NULL && break
        Libc.free(sol.structure)
        i += 1
    end
    Libc.free(ptr)
end

function free_structure_list(ptr::Ptr, num::Integer)
    for i = 1:num
        Libc.free(unsafe_load(ptr, i))
    end
    Libc.free(ptr)
end
end
import .LibRNA_Helper

# TODO: rename en_unit to unit_energy, en_int_unit to unit_energy_int
const en_unit = 1.0u"kcal/mol"
const en_int_unit = 0.01u"kcal/mol"
const unit_temperature = u"°C"

# TODO: use $(LibRNA.VRNA_MODEL_DEFAULT_GQUAD) string interpolation
# here instead of hardcoding defaults
"""
    FoldCompound(seq::AbstractString; [params, temperature], [model_details...])
    FoldCompound(msa::Vector{<:AbstractString}; [params, temperature], [model_details...])

A `FoldCompound` encapsulates nucleotide sequences, energy
parameters, and model details.

Input arguments:
- `seq`: nucleotide sequence, multiple strands are separated by an
  '&' character
- `msa`: multiple sequence alignment, for comparative folding
   (alifold). A vector of sequences which may contain multiple
   strands, denoted by '&', and gap '-' characters
- `params`: energy parameter set, legal values are `:RNA_Turner1999`,
  `:RNA_Turner2004`, `:RNA_Andronescu2007`, `:RNA_Langdon2018`.
  Default is `:RNA_Turner2004`.
- `temperature`: the temperature at which calculations are performed.
  Default is `37u"°C"`.

Model details (additional keyword arguments):
- `circular`: determines if the RNA strand is circular, i.e. the
  5'-end and 3'-end are covalently bonded. Default is `false`.
- `dangles`: how to treat dangling base pairs in multiloops and the
  exterior loop. Can be 0, 1, 2, or 3. See ViennaRNA docs for
  details. Default is `2`.
- `gquadruplex`: allow G-quadruplexes in predictions. Default is
  `false`.
- `log_ML`: use logarithmic energy model for multiloops. Default is
  `false`.
- `min_loop_length`: the minimum size of a loop (without the closing
   base pair). Default is `3`.
- `no_GU_basepairs`: disallow G-U basepairs. Default is `false`.
- `no_GU_closure`: disallow G-U basepairs as closing pairs for
  loops. Default is `false`.
- `no_lonely_pairs`: disallow isolated base pairs. Default is `false`.
- `special_hairpins`: use special hairpin energies for certain tri-,
  tetra- and hexloops. Default is `true`.
- `uniq_ML`: use unique decomposition for multiloops, needed for
  `pbacktrack` and `subopt`. Default is `false`.
"""
mutable struct FoldCompound
    ptr         :: Ptr{LibRNA.vrna_fc_s}
    uptr        :: UnsafePtr{LibRNA.vrna_fc_s}
    msa         :: Vector{String}
    msa_strands :: Vector{Vector{SubString{String}}}

    FoldCompound(seq::AbstractString; kwargs...) = FoldCompound([seq]; kwargs...)
    function FoldCompound(msa::Vector{<:AbstractString};
                          model_details::Ptr{LibRNA.vrna_md_s}=Ptr{LibRNA.vrna_md_s}(C_NULL),
                          options::Unsigned=LibRNA.VRNA_OPTION_DEFAULT,
                          params::Symbol=:RNA_Turner2004,
                          temperature::Quantity=37.0u"°C",
                          # model_details options
                          circular::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_CIRC),
                          dangles::Int=Int(LibRNA.VRNA_MODEL_DEFAULT_DANGLES),
                          gquadruplex::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_GQUAD),
                          log_ML::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_LOG_ML),
                          min_loop_size::Int=Int(LibRNA.TURN),
                          no_GU_basepairs::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_NO_GU),
                          no_GU_closure::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_NO_GU_CLOSURE),
                          no_lonely_pairs::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_NO_LP),
                          special_hairpins::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_SPECIAL_HP),
                          uniq_ML::Bool=Bool(LibRNA.VRNA_MODEL_DEFAULT_UNIQ_ML))
        if dangles < 0 || dangles > 3
            throw(ArgumentError("dangles must be 0, 1, 2, or 3"))
        end
        if min_loop_size < 0
            throw(ArgumentError("min_loop_size must be >= 0"))
        end
        if length(msa) > 1
            firstlen = length(first(msa))
            if any(s -> length(s) != firstlen, msa)
                throw(ArgumentError("all sequences in msa must have same length"))
            end
            if any(s -> contains(s, '&'), msa)
                throw(ArgumentError("multiple strands not supported in comparative mode (alifold)"))
            end
        end
        msa_strands = split.(msa, '&')
        if params == :RNA_Turner1999
            err = LibRNA.vrna_params_load_RNA_Turner1999()
        elseif params == :RNA_Turner2004
            err = LibRNA.vrna_params_load_RNA_Turner2004()
        elseif params == :RNA_Andronescu2007
            err = LibRNA.vrna_params_load_RNA_Andronescu2007()
        elseif params == :RNA_Langdon2018
            err = LibRNA.vrna_params_load_RNA_Langdon2018()
        else
            throw(ArgumentError("unknown energy parameters: $(params)"))
        end
        if err == 0
            error("Failed to load energy parameters $params_name")
        end
        temperature_nounit = ustrip(uconvert(u"°C", temperature))
        LibRNA.vrna_md_defaults_temperature(temperature_nounit)
        LibRNA.vrna_md_defaults_circ(Int(circular))
        LibRNA.vrna_md_defaults_dangles(dangles)
        LibRNA.vrna_md_defaults_gquad(Int(gquadruplex))
        LibRNA.vrna_md_defaults_logML(Int(log_ML))
        LibRNA.vrna_md_defaults_min_loop_size(min_loop_size)
        LibRNA.vrna_md_defaults_noGU(Int(no_GU_basepairs))
        LibRNA.vrna_md_defaults_noGUclosure(Int(no_GU_closure))
        LibRNA.vrna_md_defaults_noLP(Int(no_lonely_pairs))
        LibRNA.vrna_md_defaults_special_hp(Int(special_hairpins))
        LibRNA.vrna_md_defaults_uniq_ML(Int(uniq_ML))

        ptr = if length(msa) == 1
            LibRNA.vrna_fold_compound(first(msa), model_details, options)
        else
            LibRNA.vrna_fold_compound_comparative(msa, model_details, options)
        end
        ptr != C_NULL || error("vrna_fold_compound returned pointer == C_NULL")

        fc = new(ptr, UnsafePtr(ptr), msa, msa_strands)
        if fc.min_loop_size != min_loop_size
            @warn "min_loop_size was specified as $min_loop_size, but got set to $(fc.min_loop_size) by ViennaRNA"
        end
        finalizer(fc) do x
            # TODO: do we have to call vrna_mx_mfe_free or
            #       vrna_mx_pf_free here ourselves?
            LibRNA.vrna_fold_compound_free(x.ptr)
        end
    end
end

# TODO: fc.uptr.length[] |> Int
Base.length(fc::FoldCompound) = sum(size(fc))

Base.size(fc::FoldCompound) = ntuple(i -> length(first(fc.msa_strands)[i]), fc.nstrands)

function Base.getproperty(fc::FoldCompound, sym::Symbol)
    if sym == :circular
        return Bool(fc.uptr.params[].model_details.circ[])
    elseif sym == :dangles
        return Int(fc.uptr.params[].model_details.dangles[])
    elseif sym == :gquadruplex
        return Bool(fc.uptr.params[].model_details.gquad[])
    elseif sym == :has_exp_matrices
        fc.uptr.exp_matrices[] != C_NULL
    elseif sym == :log_ML
        return Bool(fc.uptr.params[].model_details.logML[])
    elseif sym == :min_loop_size
        return Int(fc.uptr.params[].model_details.min_loop_size[])
    elseif sym == :no_GU_basepairs
        return Bool(fc.uptr.params[].model_details.noGU[])
    elseif sym == :no_GU_closure
        return Bool(fc.uptr.params[].model_details.noGUclosure[])
    elseif sym == :no_lonely_pairs
        return Bool(fc.uptr.params[].model_details.noLP[])
    elseif sym == :nstrands
        return Int(fc.uptr.strands[])
    elseif sym == :params_name
        return unsafe_string(reinterpret(Ptr{UInt8}, pointer(fc.uptr.params[].param_file)))
    elseif sym == :special_hairpins
        return Bool(fc.uptr.params[].model_details.special_hp[])
    elseif sym == :temperature
        par_temperature = fc.uptr.params[].temperature[]
        md_temperature = fc.uptr.params[].model_details.temperature[]
        if par_temperature != md_temperature
            error("params temperature and model_details temperature don't agree")
        end
        return par_temperature * unit_temperature
    elseif sym == :uniq_ML
        return Bool(fc.uptr.params[].model_details.uniq_ML[])
    else
        # fallback
        getfield(fc, sym)
    end
end

Base.propertynames(fc::FoldCompound) =
    (fieldnames(typeof(fc))...,
     :circular, :dangles, :gquadruplex, :has_exp_matrices, :log_ML,
     :min_loop_size, :no_GU_basepairs, :no_GU_closure,
     :no_lonely_pairs, :nstrands, :params_name, :special_hairpins,
     :temperature, :uniq_ML)

function Base.show(io::IO, mime::MIME"text/plain", fc::FoldCompound)
    strand = "$(fc.nstrands) strand" * (fc.nstrands > 1 ? "s" : "")
    nt = "$(length(fc)) nt$(fc.nstrands > 1 ? " total" : "")"
    circ = fc.circular ? " (circular)" : ""
    comparative = length(fc.msa) > 1 ? " [comparative]" : ""
    println(io, "FoldCompound, $strand, $nt$circ$comparative")
    println(io, "  params      : $(fc.params_name)")
    println(io, "  temperature : $(fc.temperature)")
    println(io, "  options     : circular=$(fc.circular), dangles=$(fc.dangles), gquadruplex=$(fc.gquadruplex), log_ML=$(fc.log_ML),")
    println(io, "                min_loop_size=$(fc.min_loop_size), no_GU_basepairs=$(fc.no_GU_basepairs), no_GU_closure=$(fc.no_GU_closure),")
    println(io, "                no_lonely_pairs=$(fc.no_lonely_pairs), special_hairpins=$(fc.special_hairpins), uniq_ML=$(fc.uniq_ML)")
    if length(fc.msa) == 1
        for (i,s) in enumerate(first(fc.msa_strands))
            println(io,   "  strand $i    : $(s)")
        end
    else
        println(io, "  MSA")
        for (i,strands) in enumerate(fc.msa_strands)
            println(io, "      ", join(strands, " & "))
        end
    end
end

"""
    Pairtable(structure)

A `Pairtable` represents a secondary structure given in dot-bracket
notation, e.g. `(((...)))`.  In dot-bracket notation, unpaired bases
are denoted by a dot `.` and base-pairs are denoted by matching
brackets `(` and `)`.
"""
mutable struct Pairtable
    # Note: ptr[1] contains the number of elements, ptr[i+1] is the
    #       i-th element this means that functions like getindex or
    #       setindex! must use ptr[i+1] to access the i-th element
    ptr :: Ptr{Cshort}
    function Pairtable(structure::AbstractString)
        ptr = LibRNA.vrna_ptable_from_string(structure, LibRNA.VRNA_BRACKETS_ANY)
        pt = new(ptr)
        finalizer(pt) do x
            Libc.free(x.ptr)
        end
    end
end

Base.length(pt::Pairtable) = Int(unsafe_load(pt.ptr, 1))

basepairs(pt::Pairtable) = [(i, pt[i]) for i = 1:length(pt) if i < pt[i]]

# Notes
# - only works for unpseudoknotted structures
# - vrna_db_from_ptable does the same, but we would have to incur an
#   extra malloc/free to move the string to Julia, so we do it here
#   ourselves
Base.String(pt::Pairtable) =
    join(pt[i] == 0 ? '.' : pt[i] > i ? '(' : ')' for i = 1:length(pt))

function Base.show(io::IO, mime::MIME"text/plain", pt::Pairtable)
    print(io, String(pt))
end

function Base.getindex(pt::Pairtable, i)
    @boundscheck if i < 1 || i > length(pt)
        throw(BoundsError(pt, i))
    end
    unsafe_load(pt.ptr, i + 1)
end

function Base.setindex!(pt::Pairtable, val, i)
    @boundscheck begin
        len = length(pt)
        if i < 1 || val < 0 || i > len || val > len || i == val
            throw(ArgumentError("illegal basepair set: pt[$i] = $val"))
        end
    end
    if pt[i] != 0
        unsafe_store!(pt.ptr, 0, pt[i] + 1)
    end
    if val != 0
        if pt[val] != 0
            unsafe_store!(pt.ptr, 0, pt[val] + 1)
        end
        unsafe_store!(pt.ptr, i, val + 1)
    end
    unsafe_store!(pt.ptr, val, i + 1)
end

# secondary structure distance measures

"""
    bp_distance(structure1, structure2)

Base-pair distance between two secondary structures `structure1` and
`structure2`.
"""
function bp_distance(structure1::AbstractString, structure2::AbstractString)
    length(structure1) == length(structure2) ||
        throw(ArgumentError("input structures have different lengths"))
    return LibRNA.vrna_bp_distance(structure1, structure2)
end

"""
    tree_edit_dist(structure1, structure2; [hit_type])

Tree edit distance between two secondary structures `structure1` and
`structure2`.
"""
function tree_edit_dist(structure1::AbstractString, structure2::AbstractString;
                        hit_type::Integer = LibRNA.VRNA_STRUCTURE_TREE_HIT)
    # vrna_db_to_tree_string prints result to stdout, which we discard
    ts1, ts2 = redirect_stdout(devnull) do
        return LibRNA.vrna_db_to_tree_string(structure1, hit_type),
               LibRNA.vrna_db_to_tree_string(structure2, hit_type)
    end
    t1 = LibRNA.make_tree(ts1)
    t2 = LibRNA.make_tree(ts2)
    d = LibRNA.tree_edit_distance(t1, t2)
    LibRNA.free_tree(t1)
    LibRNA.free_tree(t2)
    Libc.free(ts1)
    Libc.free(ts2)
    return d
end

# equilibrium properties

# mean basepair distance of ensemble to itself

"""
    mean_bp_distance(fc)
    mean_bp_distance(sequence)

Base-pair distance of all possible pairs of secondary structures
weighted by their Boltzmann probabilities.
"""
function mean_bp_distance(fc::FoldCompound)
    fc.has_exp_matrices ||
        throw(ArgumentError("partfn(fc::FoldCompound) must be run first"))
    return LibRNA.vrna_mean_bp_distance(fc.ptr)
end

function mean_bp_distance(sequence::AbstractString)
    fc = FoldCompound(sequence)
    partfn(fc)
    return mean_bp_distance(fc)
end

# ensemble defect: mean distance of ensemble to a target structure
# TODO: check if this is exactly expected basepair distance to target struct

"""
    ensemble_defect(fc, pt)
    ensemble_defect(sequence, pt)
    ensemble_defect(fc, structure)
    ensemble_defect(sequence, structure)

Ensemble defect between a sequence given as a `FoldCompound` or
`AbstractString` in dot-bracket notation and a secondary structure
given as a `Pairtable` or as an `AbstractString`.
"""
function ensemble_defect(fc::FoldCompound, pt::Pairtable)
    fc.has_exp_matrices ||
        throw(ArgumentError("partfn(::FoldCompound) must be run first"))
    length(fc) == length(pt) ||
        throw(ArgumentError("FoldCompound and Pairtable must have equal length"))
    return LibRNA.vrna_ensemble_defect_pt(fc.ptr, pt.ptr)
end

function ensemble_defect(sequence::AbstractString, pt::Pairtable)
    fc = FoldCompound(sequence)
    partfn(fc)
    return ensemble_defect(fc, pt)
end

ensemble_defect(fc::FoldCompound, structure::AbstractString) =
    ensemble_defect(fc, Pairtable(structure))

function ensemble_defect(sequence::AbstractString, structure::AbstractString)
    fc = FoldCompound(sequence)
    partfn(fc)
    return ensemble_defect(fc, Pairtable(structure))
end

# probability of a structure in the ensemble

"""
    prob_of_structure(fc, structure)
    prob_of_structure(sequence, structure)

Boltzmann probability of a structure for a sequence.
"""
function prob_of_structure(fc::FoldCompound, structure::AbstractString)
    fc.has_exp_matrices ||
        throw(ArgumentError("must first call ViennaRNA.partfn(::FoldCompound)"))
    length(fc) == length(structure) ||
        throw(ArgumentError("FoldCompound and structure must have equal length"))
    return LibRNA.vrna_pr_structure(fc.ptr, structure)
end

function prob_of_structure(sequence::AbstractString, structure::AbstractString)
    fc = FoldCompound(sequence)
    partfn(fc)
    return prob_of_structure(fc, structure)
end

# free energy change of folding

"""
    energy(fc, structure; [verbose])
    energy(sequence, structure; [verbose])

Calculate the free energy of folding for a given `sequence` and
`structure` in dot-bracket notation.
"""
function energy(fc::FoldCompound,
                structure::AbstractString;
                verbose::Bool=false,
                verbosity_level::Integer=LibRNA.VRNA_VERBOSITY_DEFAULT)
    length(fc) == length(structure) ||
        throw(ArgumentError("sequence and structure must have equal length"))
    if verbose
        en = LibRNA.vrna_eval_structure_v(fc.ptr,
                                          structure,
                                          verbosity_level,
                                          C_NULL)
    else
        en = LibRNA.vrna_eval_structure(fc.ptr, structure)
    end
    return en * en_unit
end

function energy(sequence::AbstractString,
                structure::AbstractString;
                verbose::Bool=false,
                verbosity_level::Integer=LibRNA.VRNA_VERBOSITY_DEFAULT)
    fc = FoldCompound(sequence; options=LibRNA.VRNA_OPTION_EVAL_ONLY)
    return energy(fc, structure; verbose, verbosity_level)
end

# minimum free energy (mfe) structure

"""
    mfe(fc)
    mfe(sequence)

Calculate the minimum free energy structure and energy.
"""
function mfe(fc::FoldCompound)
    seqlen = length(fc)
    cstr_structure = Ptr{Cchar}(LibRNA.vrna_alloc(seqlen + 1))
    en_mfe = LibRNA.vrna_mfe(fc.ptr, cstr_structure)
    structure = unsafe_string(cstr_structure)
    Libc.free(cstr_structure)
    return structure, en_mfe * en_unit
end

mfe(sequence::AbstractString) =
    mfe(FoldCompound(sequence; options=LibRNA.VRNA_OPTION_MFE))

# partition function

"""
    partfn(fc)
    partfn(sequence)

Calculate the partition function, returning a string representation of
secondary structure pairing probabilities and the ensemble free energy
`-RT log Z`.  The string consists of the letters `. , | { } ( )` which
denote bases that are essentially unpaired, weakly paired, strongly
paired without upstream/downstream preference, weakly
upstream/downstream paired, or strongly upstream/downstream paired.
"""
function partfn(fc::FoldCompound)
    seqlen = length(fc)
    cstr_structure = Ptr{Cchar}(LibRNA.vrna_alloc(seqlen + 1))
    RTlogZ = LibRNA.vrna_pf(fc.ptr, cstr_structure)
    structure = unsafe_string(cstr_structure)
    Libc.free(cstr_structure)
    return structure, RTlogZ * en_unit
end

partfn(sequence::AbstractString) = partfn(FoldCompound(sequence))

# basepair probabilities

"""
    bpp(fc)
    bpp(sequence)

Calculate the matrix of base-pair probabilities.
"""
function bpp(fc::FoldCompound)
    # TODO: return upper triangular matrix type
    fc.has_exp_matrices ||
        throw(ArgumentError("must call ViennaRNA.partfn(::FoldCompound) first"))
    n = length(fc)
    p = zeros(n,n)
    f = unsafe_load(fc.ptr)
    index = f.iindx
    probs = unsafe_load(f.exp_matrices).probs
    for i = 1:n-1
        for j = i+1:n
            p[i,j] = unsafe_load(probs, unsafe_load(index, i + 1) - j + 1)
            p[j,i] = p[i,j]
        end
    end
    return p
end

function bpp(sequence::AbstractString)
    fc = FoldCompound(sequence)
    partfn(fc)
    return bpp(fc)
end

# stochastic backtrack

"""
    pbacktrack(fc; [num_samples=1])
    pbacktrack(sequence; [num_samples=1])

Sample `num_samples` secondary structures according to their Boltzmann
probabilities.
"""
function pbacktrack(fc::FoldCompound;
                    num_samples::Integer=1,
                    options::Integer=LibRNA.VRNA_PBACKTRACK_DEFAULT)
    fc.has_exp_matrices ||
        throw(ArgumentError("must call ViennaRNA.partfn(::FoldCompound) first"))
    fc.uniq_ML ||
        throw(ArgumentError("must have fc.uniq_ML == true"))
    s = LibRNA.vrna_pbacktrack_num(fc.ptr, num_samples, options)
    samples = String[]
    s == C_NULL && return samples
    i = 1
    while unsafe_load(s, i) != C_NULL
        push!(samples, unsafe_string(unsafe_load(s, i)))
        i += 1
    end
    LibRNA_Helper.free_structure_list(s, num_samples)
    return samples
end

function pbacktrack(sequence::AbstractString; num_samples::Integer=1,
                    options::Integer=LibRNA.VRNA_PBACKTRACK_DEFAULT)
    fc = FoldCompound(sequence; uniq_ML=true)
    partfn(fc)
    return pbacktrack(fc; num_samples, options)
end

# maximum expected accuracy (MEA) structure

"""
    mea(fc); [gamma=1.0])
    mea(sequence); [gamma=1.0])

Maximum expected accuracy (MEA) structure.
"""
function mea(fc::FoldCompound; gamma=1.0)
    fc.has_exp_matrices ||
        throw(ArgumentError("must call ViennaRNA.partfn(::FoldCompound) first"))
    ptr_mea   = Ptr{Cfloat}(LibRNA.vrna_alloc(sizeof(Cfloat)))
    ptr_str   = LibRNA.vrna_MEA(fc.ptr, gamma, ptr_mea)
    mea_val   = unsafe_load(ptr_mea)
    mea_struc = unsafe_string(ptr_str)
    Libc.free(ptr_mea)
    Libc.free(ptr_str)
    return mea_struc, mea_val
end

function mea(sequence::AbstractString; gamma=1.0)
    fc = FoldCompound(sequence)
    partfn(fc)
    return mea(fc; gamma)
end

# centroid structure of ensemble

"""
    centroid(fc)
    centroid(sequence)

Centroid structure of the ensemble, the secondary structure with the
smallest Boltzmann probability weighted base-pair distance to all
other structures.
"""
function centroid(fc::FoldCompound)
    fc.has_exp_matrices ||
        throw(ArgumentError("must call ViennaRNA.partfn(::FoldCompound) first"))
    ptr_dist = Ptr{Cdouble}(LibRNA.vrna_alloc(sizeof(Cdouble)))
    ptr_str = LibRNA.vrna_centroid(fc.ptr, ptr_dist)
    cen_dist = unsafe_load(ptr_dist)
    cen_struc = unsafe_string(ptr_str)
    Libc.free(ptr_dist)
    Libc.free(ptr_str)
    return cen_struc, cen_dist
end

function centroid(sequence::AbstractString)
    fc = FoldCompound(sequence)
    partfn(fc)
    return centroid(fc)
end

# suboptimal structures

"""
    subopt(fc; [delta], [sorted])
    subopt(sequence; [delta], [sorted])

All suboptimal structures with an energy not more than `delta` above
the minimum free energy structure.
"""
function subopt(fc::FoldCompound; delta::Quantity, sorted::Bool=true)
#    fc.ptr.params.model_details.uniq_ML[] == 1 ||
#        throw(ArgumentError("must have fc.params.model_details.uniq_ML == 1"))
    delta_nounit = uconvert(u"kcal/mol", delta) / en_int_unit
    delta_int = round(Int, delta_nounit)
    # TODO: warn if this float -> int conversion leads to loss of
    # precision, e.g. in the case that delta == 0.0001 kcal/mol where
    # we would have delta_int == 0
    subs = Tuple{String,Quantity}[]
    s = LibRNA.vrna_subopt(fc.ptr, delta_int, sorted, C_NULL)
    s == C_NULL && return subs
    i = 1
    while (si = unsafe_load(s, i)).structure != C_NULL
        push!(subs, (unsafe_string(si.structure), si.energy * en_unit))
        i += 1
    end
    LibRNA_Helper.free_subopt_solutions(s)
    return subs
end

function subopt(sequence::AbstractString; delta::Quantity, sorted::Bool=true)
    fc = FoldCompound(sequence; uniq_ML=true, options=LibRNA.VRNA_OPTION_MFE)
    return subopt(fc; delta, sorted)
end

"""
    subopt_zuker(fc)
    subopt_zuker(sequence)

All suboptimal structures according to the Zuker algorithm.
"""
function subopt_zuker(fc::FoldCompound)
    subs = Tuple{String,Quantity}[]
    s = LibRNA.vrna_subopt_zuker(fc.ptr)
    s == C_NULL && return subs
    i = 1
    while (si = unsafe_load(s, i)).structure != C_NULL
        push!(subs, (unsafe_string(si.structure), si.energy * en_unit))
        i += 1
    end
    LibRNA_Helper.free_subopt_solutions(s)
    return subs
end

function subopt_zuker(sequence::AbstractString)
    fc = FoldCompound(sequence; options=LibRNA.VRNA_OPTION_MFE)
    return subopt_zuker(fc)
end

# inverse folding / sequence design

"""
    inverse_fold(start, target)

Inverse folding (aka sequence design) for a given `target` structure,
starting from the sequence `start`.

The stopping criterion for the search is that the sequence has the
target structure as its minimum free energy structure.
"""
function inverse_fold(start::AbstractString, target::AbstractString)
    length(start) == length(target) ||
        throw(ArgumentError("start and target must have same size"))
    # TODO: there is no copy(::String), so this is a substitute
    #       needed because inverse_fold mutates its first argument
    design = join(collect(start))
    dist = LibRNA.inverse_fold(design, target)
    return design, dist
end

"""
    inverse_pf_fold(start, target)

Inverse folding (aka sequence design) for a given `target` structure,
starting from the sequence `start`.

The search in sequence space tries to maximise the probability of the
target structure.
"""
function inverse_pf_fold(start::AbstractString, target::AbstractString)
    length(start) == length(target) ||
        throw(ArgumentError("start and target must have same size"))
    # TODO: there is no copy(::String), so this is a substitute
    #       needed because inverse_pf_fold mutates its first argument
    design = join(collect(start))
    RT_log_p = LibRNA.inverse_pf_fold(design, target)
    return design, RT_log_p * en_unit
end

# kinetics and move sets

"""
    neighbors(fc, pt)

All neighbors of a secondary structure `pt` given as a `Pairtable`.
"""
function neighbors(fc::FoldCompound, pt::Pairtable;
                   options::Int=LibRNA.VRNA_MOVESET_DEFAULT)
    mptr = LibRNA.vrna_neighbors(fc.ptr, pt.ptr, options)
    moves = Vector{Tuple{Int,Int}}[]
    mptr == C_NULL && return moves
    i = 1
    while (mptr_i = unsafe_load(mptr, i)).pos_5 != 0 && mptr_i.pos_3 != 0
        mlist = [(mptr_i.pos_5, mptr_i.pos_3)]
        if mptr_i.next != C_NULL
            # this move is not an atomic move, so follow the next
            # pointer and collect all moves
            next = mptr_i.next
            j = 1
            while (next_j = unsafe_load(next, j)).pos_5 != 0 && next_j.pos_3 != 0
                if next_j.next != C_NULL
                    @warn "unexpected extra branching in neighbors()"
                end
                push!(mlist, (next_j.pos_5, next_j.pos_3))
                j += 1
            end
        end
        push!(moves, mlist)
        i += 1
    end
    LibRNA.vrna_move_list_free(mptr)
    return moves
end

# plotting secondary structures

"""
    plot_coords(structure; [plot_type])

Plots a secondary structure, returning the coordinates. The optional
`plot_type` parameter can be `:simple`, `:naview`, `:circular`,
`:turtle`, or `:puzzler`, with the default being `:naview`.
"""
function plot_coords(structure; plot_type::Symbol=:simple)
    _vrna_plot_coords(structure::AbstractString, cx, cy, type) =
        LibRNA.vrna_plot_coords(structure, cx, cy, type)
    _vrna_plot_coords(structure::Pairtable, cx, cy, type) =
        LibRNA.vrna_plot_coords_pt(structure.ptr, cx, cy, type)

    if length(structure) == 0
        return Float32[], Float32[]
    end
    if plot_type === :simple
        type = LibRNA.VRNA_PLOT_TYPE_SIMPLE
    elseif plot_type === :naview
        type = LibRNA.VRNA_PLOT_TYPE_NAVIEW
    elseif plot_type === :circular
        type = LibRNA.VRNA_PLOT_TYPE_CIRCULAR
    elseif plot_type === :turtle
        type = LibRNA.VRNA_PLOT_TYPE_TURTLE
    elseif plot_type === :puzzler
        type = LibRNA.VRNA_PLOT_TYPE_PUZZLER
    else
        throw(ArgumentError("unknown plot_type $plot_type, options are:" *
                            " :simple, :naview, :circular, :turtle, :puzzler"))
    end
    ptr_cx = Ptr{Ptr{Cfloat}}(LibRNA.vrna_alloc(sizeof(Ptr{Cfloat})))
    ptr_cy = Ptr{Ptr{Cfloat}}(LibRNA.vrna_alloc(sizeof(Ptr{Cfloat})))
    n = _vrna_plot_coords(structure, ptr_cx, ptr_cy, type)
    if n == 0
        throw(ErrorException("error returned by vrna_plot_coords"))
    end
    cx = unsafe_load(ptr_cx)
    cy = unsafe_load(ptr_cy)
    x = [unsafe_load(cx, i) for i = 1:n]
    y = [unsafe_load(cy, i) for i = 1:n]
    Libc.free(cx)
    Libc.free(cy)
    Libc.free(ptr_cx)
    Libc.free(ptr_cy)
    return x, y
end

# heat capacity

"""
    heat_capacity(fc, Tmin, Tmax, [Tincrement=1u"°C"]; [mpoints=2])
    heat_capacity(sequence, Tmin, Tmax, [Tincrement=1u"°C"]; [mpoints=2])

Computes the specific heat of an RNA in a temperature range, given by
`Tmin`, `Tmax`, and `Tincrement`, from the partition function by
numeric differentiation.  The parameter `mpoints` determines how
smooth the curve should be: to calculate second derivatives, a
parabola is fit to `mpoints` + 1 data points, and increasing the
`mpoints` parameter produces a smoother curve.
"""
function heat_capacity(fc::FoldCompound, Tmin::Quantity, Tmax::Quantity,
                       Tincrement::Quantity=1u"°C"; mpoints::Integer=2)
    Tmin = ustrip(uconvert(u"°C", Tmin))
    Tmax = ustrip(uconvert(u"°C", Tmax))
    Tincrement = ustrip(uconvert(u"°C", Tincrement))
    type_T = typeof(1.0f0u"°C")
    type_c = typeof(1.0f0u"kcal/mol/K")
    ptr = LibRNA.vrna_heat_capacity(fc.ptr, Tmin, Tmax, Tincrement, mpoints)
    if ptr == C_NULL
        error("failure running vrna_heat_capacity")
    end
    hcs = Tuple{type_T, type_c}[]
    i = 1
    while (h = unsafe_load(ptr, i)).temperature >= Tmin
        T = h.temperature * u"°C"
        c = h.heat_capacity * u"kcal/mol/K"
        push!(hcs, (T, c))
        i += 1
    end
    Libc.free(ptr)
    return hcs
end

function heat_capacity(sequence::AbstractString, Tmin::Quantity, Tmax::Quantity,
                       Tincrement::Quantity=1.0u"°C"; mpoints::Integer=2)
    heat_capacity(FoldCompound(sequence), Tmin, Tmax, Tincrement; mpoints)
end

include("utils.jl")
include("plot_structure.jl")

end # module ViennaRNA


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

using Markdown
using InteractiveUtils

# ╔═╡ 4bba7820-f52f-11ea-2ffd-879f7444eb9c
using Pkg, DrWatson

# ╔═╡ 4bbabd56-f52f-11ea-0da9-55fcb93b0375
begin
	@quickactivate "StatisticalRethinkingTuring"
	using StatisticalRethinking
end

# ╔═╡ 0622be1c-f52f-11ea-026d-39e1d2e6b03c
md"## Clip-03-25-26t.jl"

# ╔═╡ 4bbb41ba-f52f-11ea-18bd-039df6e7d6cb
md"### snippet 3.25"

# ╔═╡ 4bc8abfc-f52f-11ea-1a7e-515b267bb1e1
begin
	N = 10000
	w = rand(Binomial(9, 0.6), N);
	h1 = histogram(w; normalize=:probability, leg=false)
end

# ╔═╡ 4bc93f7c-f52f-11ea-10c7-1595439f2b94
md"##### Generate the samples."

# ╔═╡ 4bd64e92-f52f-11ea-3326-fb4b8eb664a4
begin
	p_grid = range(0, step=0.001, stop=1)
	prior = ones(length(p_grid))
	likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
	samples = sample(p_grid, Weights(posterior), N)
end;

# ╔═╡ 4bd71c82-f52f-11ea-1706-876cbc5db9d5
md"### snippet 3.26"

# ╔═╡ 4bdfebf0-f52f-11ea-22df-43be326fcd25
begin
	d = rand.(Binomial.(9, samples));
	h2 = histogram(d; normalize=:probability,
	  bins=-0.5:1:9.5, leg=false, xticks=0:9, bar_width=0.2)
end

# ╔═╡ 4be073c2-f52f-11ea-2ad1-5b555bbb69be
plot(h1, h2, layout=(1, 2))

# ╔═╡ 4bec4ae4-f52f-11ea-34c1-cdbc9b668632
md"## End of clip-03-25-26t.jl"

# ╔═╡ Cell order:
# ╟─0622be1c-f52f-11ea-026d-39e1d2e6b03c
# ╠═4bba7820-f52f-11ea-2ffd-879f7444eb9c
# ╠═4bbabd56-f52f-11ea-0da9-55fcb93b0375
# ╟─4bbb41ba-f52f-11ea-18bd-039df6e7d6cb
# ╠═4bc8abfc-f52f-11ea-1a7e-515b267bb1e1
# ╟─4bc93f7c-f52f-11ea-10c7-1595439f2b94
# ╠═4bd64e92-f52f-11ea-3326-fb4b8eb664a4
# ╟─4bd71c82-f52f-11ea-1706-876cbc5db9d5
# ╠═4bdfebf0-f52f-11ea-22df-43be326fcd25
# ╠═4be073c2-f52f-11ea-2ad1-5b555bbb69be
# ╟─4bec4ae4-f52f-11ea-34c1-cdbc9b668632


text:
#=
test_trees:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-29
=#

module TestTrees
    # packages
    using Test

    # external modules
    include("binary_tree.jl")
    include("max_heap.jl")
    include("min_heap.jl")
    include("trie.jl")

    # methods
    function test()
        @testset "Trees" begin
            @testset "BinaryTree" begin
                @test_nowarn MyBinaryTree.main()
            end
            @testset "MyMaxHeap" begin
                @test_nowarn MyMaxHeap.main()
            end
            @testset "MyMinHeap" begin
                @test_nowarn MyMinHeap.main()
            end
            @testset "MyTrie" begin
                @test_nowarn MyTrie.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestTrees
    TestTrees.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Stheno: CrossKernel, ZeroKernel, ConstantKernel, pairwise, EmpiricalKernel,
    _pairwise_fallback

struct FooKernel <: CrossKernel end
(::FooKernel)(x, x′) = sum(abs2, x - x′)

@testset "kernel" begin

    # Test pairwise.
    let
        rng, P, Q, D = MersenneTwister(123456), 3, 2, 4
        X, X′ = ColsAreObs(randn(rng, D, P)), ColsAreObs(randn(rng, D, Q))
        x, x′ = randn(rng, P), randn(rng, Q)
        foo = FooKernel()

        @test _pairwise_fallback(foo, X, X′) ==
            [foo(X[1], X′[1]) foo(X[1], X′[2]);
             foo(X[2], X′[1]) foo(X[2], X′[2]);
             foo(X[3], X′[1]) foo(X[3], X′[2])]
        @test _pairwise_fallback(foo, x, x′) ==
            [foo(x[1], x′[1]) foo(x[1], x′[2]);
             foo(x[2], x′[1]) foo(x[2], x′[2]);
             foo(x[3], x′[1]) foo(x[3], x′[2])]

        @test pairwise(foo, X, X′) == _pairwise_fallback(foo, X, X′)
        @test pairwise(foo, x, x′) == _pairwise_fallback(foo, x, x′)
        @test pairwise(foo, X) == pairwise(foo, X, X)
        @test pairwise(foo, x) == pairwise(foo, x, x)

        XB = BlockData([X, X′])
        @test pairwise(foo, XB, XB) ==
            vcat(hcat(pairwise(foo, X, X), pairwise(foo, X, X′)),
                 hcat(pairwise(foo, X′, X), pairwise(foo, X′, X′)))
    end

    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        XB0, XB1, XB2 = BlockData([x0, X0]), BlockData([x1, X1]), BlockData([x2, X2])

        # Tests for ZeroKernel.
        k_zero = ZeroKernel{Float64}()
        @test isstationary(k_zero)
        @test k_zero(0, 0) === zero(Float64)
        @test size(k_zero, 1) == Inf && size(k_zero, 2) == Inf
        @test size(k_zero) == (Inf, Inf)
        @test ZeroKernel{Float64}() == ZeroKernel{Float32}()
        kernel_tests(k_zero, x0, x1, x2)
        kernel_tests(k_zero, X0, X1, X2)
        kernel_tests(k_zero, XB0, XB1, XB2)

        @test ZeroKernel{Float64}() + ZeroKernel{Float64}() == ZeroKernel{Float64}()
        @test ZeroKernel{Float64}() * ZeroKernel{Float64}() === ZeroKernel{Float64}()
        @test ZeroKernel{Float64}() * EQ() == ZeroKernel{Float64}()

        # Tests for ConstantKernel.
        k_const = ConstantKernel(randn(rng))
        @test isstationary(k_const)
        @test k_const == k_const
        @test size(k_const, 1) == Inf && size(k_const, 2) == Inf
        @test size(k_const) == (Inf, Inf)
        kernel_tests(k_const, x0, x1, x2)
        kernel_tests(k_const, X0, X1, X2)
        kernel_tests(k_const, XB0, XB1, XB2)

        zro = ZeroKernel{Float64}()
        @test k_const + k_const isa ConstantKernel
        @test k_const * k_const == ConstantKernel(k_const.c^2)
        @test zro + k_const === k_const
        @test k_const + zro === k_const
        @test zro * k_const === zro
        @test k_const * zro === zro

        # Tests for EQ.
        @test isstationary(EQ())
        @test size(EQ(), 1) == Inf && size(EQ(), 2) == Inf
        @test size(EQ()) == (Inf, Inf)
        kernel_tests(EQ(), x0, x1, x2)
        kernel_tests(EQ(), X0, X1, X2)
        kernel_tests(EQ(), BlockData([x0, x0]), BlockData([x1, x1]), BlockData([x2, x2]))
        kernel_tests(EQ(), BlockData([X0, X0]), BlockData([X1, X1]), BlockData([X2, X2]))

        # Tests for Exponential.
        @test isstationary(EQ())
        @test size(Exponential(), 1) == Inf && size(Exponential(), 2) == Inf
        @test size(Exponential()) == (Inf, Inf)
        kernel_tests(Exponential(), x0, x1, x2)
        kernel_tests(EQ(), BlockData([x0, x0]), BlockData([x1, x1]), BlockData([x2, x2]))

        # Tests for Linear.
        @test !isstationary(Linear)
        a, b = Linear(randn(rng)), Linear(randn(rng))
        @test a == a && a ≠ b
        kernel_tests(a, x0, x1, x2)
        kernel_tests(a, X0, X1, X2)
        kernel_tests(a, BlockData([x0, x0]), BlockData([x1, x1]), BlockData([x2, x2]))
        kernel_tests(a, BlockData([X0, X0]), BlockData([X1, X1]), BlockData([X2, X2]))

        # Tests for Noise
        @test isstationary(Noise(randn(rng)))
        @test Noise(5.0) == Noise(5)
        kernel_tests(Noise(5.0), x0, x1, x2)
        kernel_tests(Noise(5.0), X0, X1, X2)
        # kernel_tests(Noise(5.0), XB0, XB1, XB2)

        # Tests for EmpiricalKernel.
        A_ = randn(rng, N, N)
        A = LazyPDMat(A_ * A_' + 1e-6I)
        k = EmpiricalKernel(A)
        Ds0, Ds1, Ds2 = 1:N, 1:N, 1:N-1
        kernel_tests(k, Ds0, Ds1, Ds2)
        kernel_tests(k, BlockData([Ds0, Ds0]), BlockData([Ds1, Ds1]), BlockData([Ds2, Ds2]))
        @test size(k, 1) == N && size(k, 2) == N

        @test map(k, :) == diag(A)
        @test map(k, :, :) == diag(A)
        @test pairwise(k, :) == A
        @test pairwise(k, :, :) == A

        @test_throws AssertionError AbstractMatrix(EQ())
    end

    # # Tests for Rational Quadratic (RQ) kernel.
    # @test isstationary(RQ)
    # @test RQ(1.0)(1.0, 1.0) == 1
    # @test RQ(100.0)(1.0, 1000.0) ≈ 0
    # @test RQ(1.0) == RQ(1.0)
    # @test RQ(1.0) == RQ(1)
    # @test RQ(1.0) != RQ(5.0)
    # @test RQ(1000.0) != EQ()

    # # Tests for Linear kernel.
    # @test !isstationary(Linear)
    # @test Linear(0.0)(1.0, 1.0) == 1
    # @test Linear(1.0)(1.0, 1.0) == 0
    # @test Linear(0.0)(0.0, 0.0) == 0
    # @test Linear(0.0)(5.0, 4.0) ≈ 20
    # @test Linear(2.0)(5.0, 4.0) ≈ 6
    # @test Linear(2.0) == Linear(2.0)
    # @test Linear(1.0) != Linear(2.0)

    # # Tests for Polynomial kernel.
    # @test !isstationary(Poly)
    # @test Poly(2, -1.0)(1.0, 1.0) == 0.0
    # @test Poly(5, -1.0)(1.0, 1.0) == 0.0
    # @test Poly(5, 0.0)(1.0, 1.0) == 1.0
    # @test Poly(5, 0.0) == Poly(5, 0.0)
    # @test Poly(2, 1.0) != Poly(5, 1.0)

    # # Tests for Noise kernel.
    # @test isstationary(Noise)
    # @test Noise()(1.0, 1.0) == 1.0
    # @test Noise()(0.0, 1e-9) == 0.0
    # @test Noise() == Noise()
    # @test Noise() != RQ(1.0)

    # # Tests for Wiener kernel.
    # @test !isstationary(Wiener)
    # @test Wiener()(1.0, 1.0) == 1.0
    # @test Wiener()(1.0, 1.5) == 1.0
    # @test Wiener()(1.5, 1.0) == 1.0
    # @test Wiener() == Wiener()
    # @test Wiener() != Noise()

    # # Tests for WienerVelocity.
    # @test !isstationary(WienerVelocity)
    # @test WienerVelocity()(1.0, 1.0) == 1 / 3
    # @test WienerVelocity() == WienerVelocity()
    # @test WienerVelocity() != Wiener()
    # @test WienerVelocity() != Noise()

    # # Tests for Exponential.
    # @test isstationary(Exponential)
    # @test Exponential() == Exponential()
    # @test Exponential()(5.0, 5.0) == 1.0
    # @test Exponential() != EQ()
    # @test Exponential() != Noise()

    # if check_mem

    #     # Performance checks: Constant.
    #     @test memory(@benchmark Constant(1.0) seconds=0.1) == 0
    #     @test memory(@benchmark $(Constant(1.0))(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: EQ.
    #     @test memory(@benchmark EQ() seconds=0.1) == 0
    #     @test memory(@benchmark $(EQ())(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: RQ.
    #     @test memory(@benchmark RQ(1.0) seconds=0.1) == 0
    #     @test memory(@benchmark $(RQ(1.0))(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: Linear.
    #     @test memory(@benchmark Linear(1.0) seconds=0.1) == 0
    #     @test memory(@benchmark $(Linear(1.0))(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: Linear.
    #     @test memory(@benchmark Poly(2, 1.0) seconds=0.1) == 0
    #     @test memory(@benchmark $(Poly(2, 1.0))(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: White noise.
    #     @test memory(@benchmark Noise() seconds=0.1) == 0
    #     @test memory(@benchmark $(Noise())(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: Wiener process.
    #     @test memory(@benchmark Wiener() seconds=0.1) == 0
    #     @test memory(@benchmark $(Wiener())(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: Wiener process.
    #     @test memory(@benchmark WienerVelocity() seconds=0.1) == 0
    #     @test memory(@benchmark $(WienerVelocity())(1.0, 0.0) seconds=0.1) == 0

    #     # Performance checks: Exponential
    #     @test memory(@benchmark Exponential() seconds=0.1) == 0
    #     @test memory(@benchmark $(Exponential())(1.0, 0.0) seconds=0.1) == 0
    # end
end


text:
using BinaryBuilder

name = "Qemu"
version = v"6.2.0"

# Collection of sources required to build libffi
sources = [
    ArchiveSource("https://download.qemu.org/qemu-6.2.0.tar.xz",
                  "68e15d8e45ac56326e0b9a4afa8b49a3dfe8aba3488221d098c84698bca65b45"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/qemu-*

# Patch out usage of MADV_NOHUGEPAGE which does not exist in glibc 2.12.X
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_madv_nohugepage.patch"

# Patch to include `falloc` header in `strace.c`
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_falloc.patch"

if [[ "${target}" == *-*-musl ]]; then
    # Patch to fix messy header situation on musl
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_syscall.patch"
fi

# Patch to disable tests
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_disable_tests.patch"

# Patch to properly link to rt
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_link_rt.patch"

## Patch in adapter for `clock_gettime()` on macOS 10.12-
#atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_clock_gettime.patch"
#
## Patch to fix pointer mismatch between `size_t` and `uint64_t`
#atomic_patch -p1 "${WORKSPACE}/srcdir/patches/qemu_size_uint64.patch"

# Configure, ignoring some warnings that we don't need, etc...
./configure --host-cc="${HOSTCC}" --extra-cflags="-I${prefix}/include -Wno-unused-result" --disable-cocoa --prefix=$prefix

echo '#!/bin/true ' > /usr/bin/Rez
echo '#!/bin/true ' > /usr/bin/SetFile
chmod +x /usr/bin/Rez
chmod +x /usr/bin/SetFile
make -j${nproc} || true
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("qemu-system-xtensaeb"    , :qemu_system_xtensaeb    ),
    ExecutableProduct("qemu-system-xtensa"      , :qemu_system_xtensa      ),
    ExecutableProduct("qemu-system-tricore"     , :qemu_system_tricore     ),
    ExecutableProduct("qemu-system-sparc64"     , :qemu_system_sparc64     ),
    ExecutableProduct("qemu-system-sparc"       , :qemu_system_sparc       ),
    ExecutableProduct("qemu-system-x86_64"      , :qemu_system_x86_64      ),
    ExecutableProduct("qemu-system-sh4eb"       , :qemu_system_sh4eb       ),
    ExecutableProduct("qemu-system-s390x"       , :qemu_system_s390x       ),
    ExecutableProduct("qemu-system-sh4"         , :qemu_system_sh4         ),
    ExecutableProduct("qemu-system-rx"          , :qemu_system_rx          ),
    ExecutableProduct("qemu-system-riscv32"     , :qemu_system_riscv32     ),
    ExecutableProduct("qemu-system-riscv64"     , :qemu_system_riscv64     ),
    ExecutableProduct("qemu-system-ppc64"       , :qemu_system_ppc64       ),
    ExecutableProduct("qemu-system-ppc"         , :qemu_system_ppc         ),
    ExecutableProduct("qemu-system-or1k"        , :qemu_system_or1k        ),
    ExecutableProduct("qemu-system-nios2"       , :qemu_system_nios2       ),
    ExecutableProduct("qemu-system-microblazeel", :qemu_system_microblazeel),
    ExecutableProduct("qemu-system-microblaze"  , :qemu_system_microblaze  ),
    ExecutableProduct("qemu-system-m68k"        , :qemu_system_m68k        ),
    ExecutableProduct("qemu-system-i386"        , :qemu_system_i386        ),
    ExecutableProduct("qemu-system-cris"        , :qemu_system_cris        ),
    ExecutableProduct("qemu-system-hppa"        , :qemu_system_hppa        ),
    ExecutableProduct("qemu-system-avr"         , :qemu_system_avr         ),
    ExecutableProduct("qemu-system-arm"         , :qemu_system_arm         ),
    ExecutableProduct("qemu-system-alpha"       , :qemu_system_alpha       ),
    ExecutableProduct("qemu-system-aarch64"     , :qemu_system_aarch64     ),
    ExecutableProduct("qemu-system-mips"        , :qemu_system_mips        ),
    ExecutableProduct("qemu-system-mipsel"      , :qemu_system_mipsel      ),
    ExecutableProduct("qemu-system-mips64el"    , :qemu_system_mips64el    ),
    ExecutableProduct("qemu-system-mips64"      , :qemu_system_mips64      ),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Pixman_jll"),
    Dependency("Glib_jll"; compat="2.68.1"),
    Dependency("PCRE_jll"),
    BuildDependency("Gettext_jll"),
    Dependency("libcap_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8", julia_compat="1.6")


text:
using .Gadfly
using Colors
using StatsBase

function Gadfly.plot(
    zep::Zeppelin,
    rows;
    autoklms = true,
    klms = [],
    xmin = 0.0,
    xmax = missing,
    norm = :None,
    yscale = 1.05,
    ytransform = identity,
    style = NeXLSpectrum.NeXLSpectrumStyle,
    palette = NeXLCore.NeXLPalette,
)
    plot(
        Spectrum[filter(s->!ismissing(s), map(r -> spectrum(zep, r, false),rows))...],
        autoklms = autoklms,
        klms = klms,
        xmin = xmin,
        xmax = xmax,
        norm = norm,
        yscale = yscale,
        ytransform = ytransform,
        style = style,
        palette = palette,
    )
end

function Gadfly.plot(zeps::AbstractArray{Zeppelin}; offset=(0.0,0.0), point_size=1.5pt, coords=missing)
    names = NeXLSpectrum.name.(zeps)
    colors = map(i->NeXLPalette[1 + (i-1) % length(NeXLPalette)], eachindex(zeps))
    clsCx=collect(merge( [ StatsBase.countmap(repr.(z[:,:CLASS])) for z in zeps]...))
    sort!(clsCx, lt=(a,b)->a[2]<b[2], rev=true)
    shlen = length(Theme().point_shapes)
    if length(clsCx) >= shlen
        clsCx = clsCx[1:(shlen-1)]
        push!(clsCx,"Other"=>0)
    end
    shIdx=Dict(clsCx[i][1]=>i for i in eachindex(clsCx))
    xmin = minimum(z->minimum(z.data[:,:XABS]), zeps)
    xmax = maximum(z->maximum(z.data[:,:XABS]), zeps)
    ymin = minimum(z->minimum(z.data[:,:YABS]), zeps)
    ymax = maximum(z->maximum(z.data[:,:YABS]), zeps)
    dmax = ceil(max(xmax-xmin, ymax-ymin)/2)
    xc, yc = round((xmax+xmin)/2), round((ymax+ymin)/2)
    coords = Coord.cartesian(xmin=xc-dmax,xmax=xc+dmax,ymin=yc-dmax,ymax=yc+dmax)
    plot(
        Theme(point_size=1.5pt),
        map(enumerate(zeps)) do (i,z)
            layer(
                x=map(r->r.XABS+(i-1)*offset[1], eachrow(z.data)),
                y=map(r->r.YABS+(i-1)*offset[2], eachrow(z.data)),
                shape=map(r->get(shIdx, repr(r.CLASS), length(clsCx)), eachrow(z.data)),
                Geom.point,
                Theme(default_color=colors[i], point_size=point_size)
            )
        end...,
        Guide.xlabel("X (mm)"), Guide.ylabel("Y (mm)"),
        Guide.manual_color_key("Particle Data Set", names, colors),
        Guide.shapekey(title="Particle Class", labels=map(cc->cc[1],clsCx)),
        coords
    )
end

function Gadfly.plot(::Type{Histogram}, zeps::AbstractArray, sym::Symbol; bincount=40, limits=missing)
    names = NeXLSpectrum.name.(zeps)
    df = vcat(
        map(enumerate(zeps)) do (i, z)
            DataFrame(Dataset=fill(names[i],nrow(z.data)), v=z.data[:,sym])
        end...
    )
    if ismissing(limits)
        limits = ( min=floor(minimum(z->minimum(z.data[:,sym]), zeps)), max=ceil(maximum(z->maximum(z.data[:,sym]), zeps)) )
    end
    h = Geom.histogram(bincount=bincount, position=:identity, limits=limits)
    plot(df, x=:v, color=:Dataset, alpha=[0.5], h,
        Guide.xlabel(repr(sym)[2:end]), Guide.ylabel("Count")
    )
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
@doc raw"""
```
Pentland(slant::Real = Inf, tilt::Real = Inf)
```
Creates a `Pentland` fucntor which attempts to produce a heightmap from a grayscale
image using Pentland's algorithm. Under the assumptions of Lambertian surface,
orthographic projections, the surface being illuminated by distant light sources,
the surface is not self-shadowing and the surface has constant albedo, hence it
can be ignored. The algorithm employs Tayler series expansion and Fourier transforms
to compute a non-iterative solution.
# Output
Returns a `Pentland` functor which can be used to run the Pentland algorithm to
reconstruct a surface form a grayscale image.
# Details
In Pentlands algorithm the image irradiance is defined as:
```math
E(x,y)=R(p,q)=\dfrac{\rho(i_xp+i_yq-i_z)}{\sqrt{1+p^2+q^2}}=\dfrac{p\sin\sigma\cos
\tau+q\sin\sigma\sin\tau+\cos\sigma}{\sqrt{1+p^2+q^2}}
```
This can be reduced using the Taylor series expansion about ``p=p_0`` and ``p=p_0``
and ignoring the higher order terms becomes:
```math
E(x,y)=R(p,q)\approx R(p_0,q_0)+(p-p_0)\dfrac{\delta R}{\delta p}(p_0,q_0)+(q-q_0)
\dfrac{\delta R}{\delta q}(p_0,q_0)
```
which for ``p_0=q_0=0`` further reduces to:
```math
E(x,y)\approx\cos\sigma+p\cos\tau\sin\sigma+q\sin\tau\sin\sigma
```
This gives the following transform identities:
```math
\begin{gathered}
p=\dfrac{\delta}{\delta x}Z(x,y)\xleftrightarrow{\Im}(-j\omega_x)F_z(
\omega_x,\omega_y)\\q=\dfrac{\delta}{\delta y}Z(x,y)\xleftrightarrow{\Im}(-j
\omega_y)F_z(\omega_x,\omega_y)
\end{gathered}
```
By taking the Fourier transform of both sides if ``E(x,y)`` yields the following:
```math
F_E=(-j\omega_x)F_z(\omega_x,\omega_y)\cos\tau\sin\sigma+(-j\omega_y)F_z(\omega_x
,\omega_y)\sin\tau\sin\sigma
```
where ``F_z`` is the Fourier transform of ``Z(x,y)``.

These can be rearranged, and the Inverse Fourier transform used to recover the
surface ``Z(x,y)`` as per the following:
```math
\begin{gathered}
F_E=F_z(\omega_x,\omega_y)[-j\omega_x\cos\tau\sin\sigma-j\omega_y\sin\tau\sin
\sigma]\\\Rightarrow F_z(\omega_x,\omega_y)=\dfrac{F_E}{-j\omega_x\cos\tau\sin
\sigma-j\omega_y\sin\tau\sin\sigma}\\Z(x,y)=\Im^{-1}\{F_z(\omega_x,\omega_y)\}
\end{gathered}
```
# Arguments
The function arguments are described in more detail below.
## `slant`
A `Real` that specifies the slant value to be used by the algorithm. The `slant`
should be a value in the range [0,π/2]. If `slant` is specified to must the `tilt`.
Will defult to `Inf` which will trigger `estimate_img_properties` to be run.
## `tilt`
A `Real` that specifies the tilt value to be used by the algorithm. The `tilt`
should be a value in the range [0,2π]. If `tilt` is specified to must
the `slant`. Will defult to `Inf` which will trigger `estimate_img_properties` to
be run.
!!! note
    If `slant` and `tilt` are not defined (i.e. have defulted to `Inf`) they
    will be calculated at runtime using `estimate_img_properties`.
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.1,0,0.9])

#calculate the heightmap
pentland = Pentland(slant = 0.1106572211739, tilt = 0.0)
Z = pentland(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. A. Pentland, "Shape Information From Shading: A Theory About Human Perception," [1988 Proceedings] Second International Conference on Computer Vision, Tampa, FL, USA, 1988, pp. 404-413. [doi: 10.1109/CCV.1988.590017](https://doi.org/10.1109/ccv.1988.590017 )
"""
function (algorithm::Pentland)(img::AbstractArray)
    #find illumination and albedo
    σ = algorithm.slant
    τ = algorithm.tilt
    if σ == Inf || τ == Inf
        ρ, I, σ, τ = estimate_img_properties(img)
    end

    E = Complex{Float64}.(reinterpret(Float64, img))

    #take Fourier transform
    fft!(E)
    M, N = size(E)

    #setup wx and wy
    wx, wy = setup_transform_values_pentland(M, N)

    #using the illumination direction calculate the transformed Z
    Zₜ = zeros(Complex{Float64}, size(E))
    for i in CartesianIndices(Zₜ)
        Zₜ[i] = E[i] / (-1im * wx[i] * cos(τ) * sin(σ) - 1im * wy[i] * sin(τ)
            * sin(σ))
    end

    #recover Z
    ifft!(Zₜ)
    Z = zeros(Float64, size(E))
    for i in CartesianIndices(Z)
        Z[i] = abs(Zₜ[i])
    end
    return Z
end


text:
using DataFrames
# using Plots
using PyCall
import PyPlot

PyCall.pygui(:tk)

PyPlot.plot([12,3,4,5,6,7], [2,4,6,7,9,0])
PyPlot.savefig("plots.png")
PyPlot.show()

println("hello world")

readline()

max_stars_repo_name: Yobmod/dmlmung

max_issues_repo_name: Yobmod/dmlmung

max_forks_repo_name: Yobmod/dmlmung


text:
function Base.rand(method::S, sdm::M) where {S<:SpatialSampler,M<:AbstractMatrix}
    _generate!(method, sdm)
end


text:
# # Heat Equation
#
# ## Introduction
#
# The heat equation is the "Hello, world!" equation of finite elements.
# Here we solve the equation on a unit square, with a uniform internal source.
# The strong form of the (linear) heat equation is given by
#
# ```math
#  -\nabla \cdot (k \nabla u) = f  \quad x \in \Omega,
# ```
#
# where $u$ is the unknown temperature field, $k$ the heat conductivity,
# $f$ the heat source and $\Omega$ the domain. For simplicity we set
# $k = 1$. We will consider homogeneous Dirichlet boundary conditions such that
# ```math
# u(x) = 0 \quad x \in \partial \Omega,
# ```
# where $\partial \Omega$ denotes the boundary of $\Omega$.
#
# The resulting weak form is given by
# ```math
# \int_{\Omega} \nabla v \cdot \nabla u \ d\Omega = \int_{\Omega} f\cdot v \ d\Omega,
# ```
# where $v$ is a suitable test function.
#-
# ## Commented Program
using HDiscontinuousGalerkin
using Tensors
using SparseArrays

# We start  generating a simple grid with 20x20 quadrilateral elements
# using `generate_grid`. The generator defaults to the unit square,
# so we don't need to specify the corners of the domain.
mesh = rectangle_mesh(TriangleCell, (10,10), Vec{2}((0.0,0.0)), Vec{2}((1.0,1.0)));

# ### Initiate function Spaces
dim = 2
Wh = ScalarFunctionSpace(mesh, ContinuousLagrange{dim,RefTetrahedron,1}(); update_face_data = false)

# Declare variables
u_h = TrialFunction(Wh)

# ### Degrees of freedom
# Next we need to define a `DofHandler`, which will take care of numbering
# and distribution of degrees of freedom for our approximated fields.
# We create the `DofHandler` and then add a single field called `u`.
# Lastly we `close!` the `DofHandler`, it is now that the dofs are distributed
# for all the elements.
dh = DofHandler([u_h], mesh);

# Now that we have distributed all our dofs we can create our tangent matrix,
# using `create_sparsity_pattern`. This function returns a sparse matrix
# with the correct elements stored.
K = create_sparsity_pattern(dh);

# We can inspect the pattern using the `spy` function from `UnicodePlots.jl`.
# By default the stored values are set to $0$, so we first need to
# fill the stored values, e.g. `K.nzval` with something meaningful.
using UnicodePlots
fill!(K.nzval, 1.0);
spy(K; height = 15)

# ### Boundary conditions
dbc = Dirichlet(u_h, dh, "boundary", [0.0])

# ### RHS function
f(x::Vec{dim}) = 2*π^2*sin(π*x[1])*sin(π*x[2])

# ### Assembling the linear system
# Now we have all the pieces needed to assemble the linear system, $K u = f$.
function doassemble(Wh, K::SparseMatrixCSC, dh::DofHandler)
    # We allocate the element stiffness matrix and element force vector
    # just once before looping over all the cells instead of allocating
    # them every time in the loop.
    n_basefuncs = getnbasefunctions(Wh)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    b = zeros(ndofs(dh))
    cell_dofs = Vector{Int}(undef, ndofs_per_cell(dh))

    # Next we define the global force vector `f` and
    # create an assembler. The assembler
    # is just a thin wrapper around `f` and `K` and some extra storage
    # to make the assembling faster.
    assembler = start_assemble(K, b)

    # It is now time to loop over all the cells in our grid.
    @inbounds for (cellcount, cell) in enumerate(CellIterator(mesh))
        # Always remember to reset the element stiffness matrix and
        # force vector since we reuse them for all elements.
        fill!(Ke, 0)
        fill!(fe, 0)
        reinit!(Wh, cell)

        # It is now time to loop over all the quadrature points in the cell and
        # assemble the contribution to `Ke` and `fe`. The integration weight
        # can be queried from `cellvalues` by `getdetJdV`.
        for q_point in 1:getnquadpoints(Wh)
            dΩ = getdetJdV(Wh, q_point)
            fh = function_value(f, Wh, cell, q_point)
            # For each quadrature point we loop over all the (local) shape functions.
            # We need the value and gradient of the testfunction `v` and also the gradient
            # of the trial function `u`. We get all of these from `cellvalues`.
            for i in 1:n_basefuncs
                v  = shape_value(Wh, q_point, i)
                ∇v = shape_gradient(Wh, q_point, i)
                fe[i] += fh*v * dΩ
                for j in 1:n_basefuncs
                    ∇u = shape_gradient(Wh, q_point, j)
                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                end
            end
        end
        # The last step in the element loop is to assemble `Ke` and `fe`
        # into the global `K` and `f` with `assemble!`.
        celldofs!(cell_dofs, dh, cell)
        assemble!(assembler, cell_dofs, fe, Ke)
    end
    return K, b
end
#md nothing # hide

# ### Solution of the system
# The last step is to solve the system. First we call `doassemble`
# to obtain the global stiffness matrix `K` and force vector `f`.
K, b = doassemble(Wh, K, dh);

# To account for the boundary conditions we use the `apply!` function.
# This modifies elements in `K` and `f` respectively, such that
# we can get the correct solution vector `u` by using `\`.
apply!(K,b,dbc);
u = K \ b;

reconstruct!(u_h, u, dh)
#
# # ### Compute errors
u_ex(x::Vec{dim}) = sin(π*x[1])*sin(π*x[2])
Etu_h = errornorm(u_h, u_ex)
Etu_h <= 0.0002
#
### Plot Solution
# TODO: Turn into a recipe
using Makie, CairoMakie
#AbstractPlotting.__init__()
using ColorSchemes
connectivity = getcells_matrix(mesh)
coordinates = get_vertices_matrix(mesh)
color = u
scene = Scene(resolution = (400, 200), colorrange = (0.0, 1.0), colormap = ColorSchemes.viridis.colors)
#poly!(scene, coordinates, connectivity, color = color, strokecolor = (:black, 0.6), strokewidth = 1)
for row in 1:size(connectivity)[1]
	#read coordinates
	points = node(:poly, Point2f0[coordinates[i,:] for i in connectivity[row,:]])
	colors = get(ColorSchemes.viridis, [color[i] for i in connectivity[row,:]])
	poly!(scene, points, strokewidth = 1, strokecolor = (:black,0.6), color = colors, show_axis = true, scale_plot = false)
end
display(scene)


text:
function merge(A,B)
    C = promote_type(eltype(A), eltype(B))[]
    i, j = 1, 1
    n, m = 1, 1
    while i <= length(A)  || j <= length(B)
        while i <= length(A) && ( A[i] <= B[m] || j > length(B))
            push!(C, A[i])
            i += 1
            n += 1
        end
        if n > length(A)
            n -= 1
        end
        while j <= length(B) && (B[j] <= A[n] || i > length(A))
            push!(C,B[j])
            j += 1
            m += 1
        end
        if m > length(B)
            m -= 1
        end
    end

    return C
end

function merge1(A,B)
    C = Vector{promote_type(eltype(A), eltype(B))}(undef, length(A) + length(B))
    k = 1
    i, j = 1, 1
    n, m = 1, 1
    while i <= length(A)  || j <= length(B)
        while i <= length(A) && ( A[i] <= B[m] || j > length(B))
            C[k] = A[i]
            i += 1
            n += 1
            k += 1
        end
        if n > length(A)
            n -= 1
        end
        while j <= length(B) && (B[j] <= A[n] || i > length(A))
            C[k] = B[j]
            j += 1
            m += 1
            k += 1
        end
        if m > length(B)
            m -= 1
        end
    end

    return C
end

function merge!(A, B, C)
    i, j, k = 1, 1, 1
    D = promote_type(eltype(A), eltype(B),eltype(C))[]
    while i <= length(A) || j <= length(B) || k <= length(C)
        min_El = typemax(promote_type(eltype(A), eltype(B),eltype(C)))
        if i <= length(A) && A[i] <= min_El
            min_El = A[i]
            determinator = 1
        end
        if j <= length(B) && B[j] <= min_El
            min_El = B[j]
            determinator = 2
        end
        if k <= length(C) && C[k] <= min_El
            min_El = C[k]
            determinator = 3
        end
        push!(D, min_El)
        if determinator == 1
            i += 1
        elseif  determinator == 2
            j += 1
        else
            k += 1
        end
    end
    C = D
    return C
end


text:
mutable struct Food
    x::Real
    y::Real
    color::String
    Food() = new(0, 0)
end

function spawn!(food::Food)
    food.x = rand(0:w - 1)
    food.y = rand(0:h - 1)
    food.color = rand(["orange", "yellow", "yellowgreen", "indigo"])
end

function draw(food::Food)
    sethue(food.color)
    rect(food.x, food.y, 1, 1, :fill)
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
# returns range-reduced x in [0,pi] for odd n
function range_reduce_odd(x::Float64)
    if x < zero(x)
        x = -x
    end

    if x >= 2.0*pi
        x = mod2pi(x)
    end

    if x > pi
        p0 = 6.28125
        p1 = 0.0019353071795864769253
        x = (p0 - x) + p1
    end

    x
end

# returns (x, sign) with range-reduced x in [0,pi] for even n
function range_reduce_even(x::Float64)
    sgn = one(x)

    if x < zero(x)
        x = -x
        sgn = -one(x)
    end

    if x >= 2.0*pi
        x = mod2pi(x)
    end

    if x > pi
        p0 = 6.28125
        p1 = 0.0019353071795864769253
        x = (p0 - x) + p1
        sgn = -sgn
    end

    (x, sgn)
end

# returns (x, sign) with range-reduced x in [0,pi]
function range_reduce(n::Integer, x::Float64)
    if iseven(n)
        range_reduce_even(x)
    else
        (range_reduce_odd(x), one(x))
    end
end


text:
# Since we are not installing the package, we manually push the
# src dir to load path so that the supporting code can be used.
# Note we use @everywhere to ensure all the worker knows the correct
# load path
@everywhere append!(LOAD_PATH,
                    [joinpath(dirname(@__FILE__), x) for x in ("src", "algor")])

using DistLearn
using ASGD
using COCOA

if length(ARGS) < 1
  println("Usage:")
  println()
  println("  julia -p NWORKER run.jl (asgd|cocoa|...) --options...")
  println()
  println("e.g.  julia -p 2 run.jl asgd --help")
  exit(-1)
end

if nprocs() <= 1
  error("Please start julia with multiple workers.")
end

algor_name = ARGS[1]
cmdline = ARGS[2:end]
s = DistLearn.make_arg_parser()

if algor_name == "asgd"
  args = ASGD.parse_cmdline(cmdline, s)
  ASGD.asgd(args)
elseif algor_name == "cocoa"
  args = COCOA.parse_cmdline(cmdline, s)
  COCOA.cocoa(args)
else
  error("Unknown algorithm: ", algor_name)
end


text:
#' ## What is learning?
#'
#' Computers read data, as we saw in notebooks 1 and 2. We can then build
#' functions that model that data to make decisions, as we saw in notebooks
#' 3 and 5.
#'
#' But how do you make sure that the model actually fits the data well?
#' In the last notebook, we saw that we can fiddle with the parameters of our
#' function defining the model to reduce the loss function. However,
#' we don't want to have to pick the model parameters ourselves.
#' Choosing parameters ourselves works *well enough* when we have a
#' simple model and only a few data points, but can quickly become
#' extremely complex for more detailed models and larger data sets.
#'
#' Instead, we want our machine to *learn* the parameters that fit the model
#' to our data, without needing us to fiddle with the parameters ourselves.
#' In this notebook, we'll talk about the "learning" in machine learning.

#' ### Motivation: Fitting parameters by hand
#'
#' Let's go back to our example of fitting parameters from notebook 3.
#' Recall that we looked at whether the amount of green in the pictures
#' could distinguish between an apple and a banana, and used a sigmoid
#' function to model our choice of "apple or banana" using the amount of
#' green in an image.

using Plots
gr()
using Images, Statistics

sigma(x,w,b) = 1/(1 + exp(-w*x + b))

apple = load("data/10_100.jpg")
banana = load("data/104_100.jpg")
apple_g = mean(Float64.(green.(apple)))
banana_g = mean(Float64.(green.(banana)))

w = 10.0 # Try manipulating w between 0 and 30 to see how the plot changes
b = 15.0 # Try manipulating b bewteen 0 and 30

pl = plot(x->sigma(x,w,b), xlims=(0, 1))
scatter!([apple_g], [0.0], label="Apple")
scatter!([banana_g], [1.0], label="Banana")

#' Intuitively, how did you tweak the sliders so that way the model sends
#' apples to 0 and bananas to 1? Most likely, you did the following:

#' #### Move the sliders a bit, see whether the curve moves in the right
#' direction, and if it did, keep doing it.

#' For a machine, "learning" is that same process, translated into math!

#' ## "Learning by nudging": The process of descent
#'
#' Let's start to formalize this idea. In order to push the curve in the
#' "right direction", we need some measurement of "how right" and "how wrong"
#' the model is. When we translate the idea of a "right direction" into math,
#' we end up with a **loss function**, `L(w, b)`, as we saw in notebook 5.
#' We say that the loss function is lowest when the model `σ(x, w, b)`
#' performs the best.
#'
#' Now we want to create a loss function that is the lowest when the apple is
#' at `0` and the banana is at `1`. If the data (the amount of green) for
#' our apple is $x_1$, then our model will output $σ(x_1,w, b)$ for our apple.
#' So, we want the difference $0 - σ(x_1, w, b)$ to be small. Similarly, if our
#' data for our banana (the banana's amount of green) is $x_2$, we want the
#' difference $1 - σ(x_2, w, b)$ to be small.
#'
#' To create our loss function, let's add together the squares of the
#' difference of the model's output from the desired output for the apple
#' and the banana. We get
#'
#' $$ L(w,b) = (0 - σ(x_1, w, b))^2 + (1 - σ(x_2, w, b))^2. $$
#'
#' $L(w, b)$ is lowest when it outputs `0` for the apple and `1` for the
#' banana, and thus the cost is lowest when the model "is correct".
#'
#' We can visualize this function by plotting it in 3D with the
# `surface` function or in 2D with contour lines

# plotly() # The plotly backend is nice for 3d surface plots
# gr() # The GR backend is good for faster interactive plots

L(w,b) = (0 - sigma(apple_g,w,b))^2 + (1 - sigma(banana_g,w,b))^2

w_range = 10:0.1:13
b_range = 0:1:20
L_values = [L(w,b) for b in b_range, w in w_range]

w = 11.5 # Try manipulating w with values from w_range (between 10 and 13)
b = 7 # Try manipulating b with values from b_range (between 0 and 20)
##     p1 = surface(w_range, b_range, L_values, xlabel="w", ylabel="b", cam=(70,40), cbar=false, leg=false)
##     scatter!(p1, [w], [b], [L(w,b)+1e-2], markersize=5, color = :blue)
p1 = contour(w_range, b_range, L_values, levels = 0.05:0.1:1,
    xlabel="w", ylabel="b", cam=(70,40), cbar=false, legend=false)
scatter!(p1, [w], [b], markersize=5, color=RGB(0.05, 0.05, 0.90))

p2 = plot(x->sigma(x,w,b), 0, 1, label="Model", legend=(0.1,0.9), lw=3)
scatter!(p2, [apple_g], [0.0], label="Apple", markersize=8)
scatter!(p2, [banana_g], [1.0], label="Banana", markersize=8, xlim=(0,1),
    ylim=(-0.1,1.1))
plot(p1, p2, layout=(2,1))

#' The blue ball on the 3D plot shows the current parameter choices,
#' plotted as `(w,b)`. Shown below the 3D plot is a 2D plot of the
#' corresponding model with those parameters. Notice that as the blue ball
#' rolls down the hill, the model becomes a better fit. Our loss function
#' gives us a mathematical notion of a "hill", and the process of
#' "learning by nudging" is simply rolling the ball down that hill.

#' To do this mathematically, we need to know which direction is "downhill".
#' Recall from calculus that the derivative of `L` with respect to `b` tells
#' you how `L` changes when `b` changes. Thus to roll downhill, we should go
#' in the direction where the derivative is negative (the function goes down)
#' for each parameter. This direction is the negative of what's called the
#' **gradient**, $\nabla L$. This means that the "learn by nudging method"
#' can be rephrased in mathematical terms as:
#'
#' 1. Calculate the gradient
#' 2. Move a little bit in the direction of the negative gradient
#' 3. Repeat

#' This process of rolling the ball in the direction of the negative gradient
#' is called **gradient descent**; written mathematically, it is
#'
#' $$p_{n+1} = p_n - \eta \nabla L(p_n).$$
#'
#' Here, $p_n$ represents the vector of current parameters
#' $(w, b)$; $\nabla L(p_n)$ is the gradient of the loss function, given
#' those parameters. We start from $p_n$ and change it by $\eta \nabla L(p_n)$,
#' where $\eta$ is a small step size that determines how far we move the
#' parameters in the direction of the negative gradient; notice that if you
#' step too far, you'll overshoot the minimum!. The result is
#' $p_{n+1}$, the new vector of parameters.
#'
#' [Picture of Gradient Descent Vectors]
#'
#' If we repeat this process, then we will end up at parameters where the
#' model correctly labels apples as `0` and bananas as `1`. When this happens,
#' the model has learned from the data and can then read pictures and tell
#' you whether they are apples or bananas!

#' #### Exercise 1
#'
#' Use the following terms to fill in the sentences below.
#' Terms may be used more than once or not at all:
#' > gradient, loss function, derivative, gradient descent, learning.
#'
#' * We can think of a _(A)_ as a 1D version of a _(B)_.
#' * We can visualize a _(C)_ as a hill.
#' * In the explanation above, rolling downhill is called _(D)_ and means traveling along the _(E)_.
#' * To quantify the correctness of a model we use a _(F)_.
#' * When our program can minimize a _(G)_ on its own, we say it is _(H)_.
#'
#'

#'
#' A)

#' B)

#' C)

#' D)

#' E)

#' F)

#' G)

#' H)