text:
# enable parallel computing
using Distributed
addprocs(...) # set number of cores you want to make available for the computations

# needed packages
@everywhere using Distributions, Bootstrap, Statistics, LinearAlgebra, SharedArrays
using DataFrames, HypothesisTests, DelimitedFiles, StatsBase, Colors, Gadfly, Cairo, Compose

# default graphics
Gadfly.push_theme(:default)
set_default_plot_size(9inch, 9inch/MathConstants.golden)

function gen_brew_colors(n) # to create your own colors, here based on one of the brewer series
    cs = distinguishable_colors(n,
        [colorant"#66c2a5", colorant"#fc8d62", colorant"#8da0cb", colorant"#e78ac3",
            colorant"#a6d854", colorant"#ffd92f", colorant"#e5c494", colorant"#b3b3b3"],
        lchoices=Float64[58, 45, 72.5, 90],
        transform=c->deuteranopic(c, 0.1),
        cchoices=Float64[20,40],
        hchoices=[75,51,35,120,180,210,270,310]
    )
    convert(Vector{Color}, cs)
end

# set parameters, define priors, etc.
@everywhere begin
    const numb_hyp = 11
    const numb_toss = 500
    const numb_sim = 1000
    const prior = fill(Float32(1/numb_hyp), numb_hyp)
    const likelihood_heads = range(0f0, stop=1, length=numb_hyp)
    const likelihood_tails = range(1f0, stop=0, length=numb_hyp)
end

@everywhere datFunc(bias) = rand(Bernoulli(bias), numb_toss)

# Bayes' rule
@everywhere function b_upd(probs::Array{Float32,1}, dat::Array{Bool,1}, toss_num::Int64)
    if dat[toss_num] == true
        @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)
    else
        @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)
    end
end

# EXPL
@everywhere function expl_upd(probs::Array{Float32,1}, dat::Array{Bool, 1}, toss_num::Int64, bonus::Float32=0.1)
    val::Float32 = mean(dat[1:toss_num]) * 10 + 1
    vec::Array{Float32,1} = if dat[toss_num] == true
            @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)
        else
            @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)
        end

    if val % 1 == .5
        vec[floor(Int64, val)] += .5*bonus
        vec[ceil(Int64, val)] += .5*bonus
    else
        vec[round(Int64, val, RoundNearestTiesAway)] += bonus
    end

    return vec / (1 + bonus)
end

# Good's rule
@everywhere function good_bonus(probs::Array{Float32,1}, res::Bool, λ=2) # with λ=2, we obtain the rule L2 from Douven and Schupbach, Frontiers ...

    pE::Float32 = res == true ? dot(probs, likelihood_heads) : dot(probs, likelihood_tails)
    gb::Array{Float32,1} = res == true ? log.(likelihood_heads ./ pE) : log.(likelihood_tails ./ pE)

    function rsc(i)
        if i >= 0
            1 - exp(2λ^2 * -i^2)
        else
            -1 + exp(2λ^2 * -i^2)
        end
    end

    return map(rsc, gb)

end

@everywhere function good_upd(probs::Array{Float32,1}, dat::Array{Bool,1}, toss_num::Int64, γ::Float32=0.1)

    res::Bool = dat[toss_num]

    probvec::Array{Float32,1} = if res == true
        @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)
    else
        @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)
    end

    goodvec::Array{Float32,1} = probvec + γ .* (probvec .* good_bonus(probs, res))

    return goodvec / sum(goodvec)

end

# Popper's rule
@everywhere function pop_bonus(probs::Array{Float32,1}, res::Bool)

    pE::Float32 = res == true ? dot(probs, likelihood_heads) : dot(probs, likelihood_tails)
    pb::Array{Float32, 1} = res == true ? (likelihood_heads .- pE) ./ (likelihood_heads .+ pE) : (likelihood_tails .- pE) ./ (likelihood_tails .+ pE)

 end

@everywhere function pop_upd(probs::Array{Float32,1}, dat::Array{Bool, 1}, toss_num::Int64, γ::Float32=0.1)

    res::Bool = dat[toss_num]

    probvec::Array{Float32,1} = if res == true
        @. (probs * likelihood_heads) / $dot(probs, likelihood_heads)
    else
        @. (probs * likelihood_tails) / $dot(probs, likelihood_tails)
    end

    popvec::Array{Float32,1} = probvec + γ .* (probvec .* pop_bonus(probs, res))

    return popvec / sum(popvec)

end

# simulations
@everywhere const numb_agents = 200
@everywhere const numb_generations = 250
const numb_simulations = 50

@everywhere function survWei(upds::Array{Float32,2}, # modeling probability of death, based on Weibull distribution
                             hyp::Int64,
                             a::Float64,
                             b::Float64,
                             thresh::Float32,
                             shape::Float64=rand(Uniform(.5, 5)),
                             scale::Float64=rand(Uniform(50, 250))
                             )

    t = something(findfirst(upds .> thresh), (numb_toss, 0)) # where in the matrix with probability updates do we find the first value above thresh?
    c = t[2]
    p = t[1]

    # cdf(Weibull(shape, scale), p) below gives the probability of death at the relevant time

    if c == hyp
        1 - (cdf(Weibull(shape, scale), p) / a) # probability goes down if right intervention is made (which is made when the truth is assigned a probability above thresh)
    elseif c == 0
        1 - cdf(Weibull(shape, scale), numb_toss + 1) # if no intervention is made, output survival probability at last time step
    else
        (1 - cdf(Weibull(shape, scale), p)) / b # probability goes down if wrong intervention is made (which happens if a false hypothesis is assigned a probabilty above thresh)
    end
end

@everywhere function survGam(upds::Array{Float32,2}, # modeling probability of death, based on Gamma distribution
                             hyp::Int64,
                             a::Float64,
                             b::Float64,
                             thresh::Float32,
                             shape::Float64=rand(Uniform(10, 16)),
                             scale::Float64=rand(Uniform(10, 16))
                             )

    t = something(findfirst(upds .> thresh), (numb_toss, 0)) # where in the matrix with probability updates do we find the first value above thresh?
    c = t[2]
    p = t[1]

    if c == hyp
        1 - (cdf(Gamma(shape, scale), p) / a) # the probability goes down if the right intervention is made (and the right intervention is made if the truth is assigned a probability above thresh)
    elseif c == 0
        1 - cdf(Gamma(shape, scale), numb_toss + 1) # if no intervention is made, output survival probability at last time step
    else
        (1 - cdf(Gamma(shape, scale), p)) / b # the probability goes down if the wrong intervention is made (which happens if a false hypothesis is assigned a probabilty above thresh)
    end
end

@everywhere function patient(rule_index::Float32, c_value::Float32, thresh::Float32, dist::Function)

    rand_hyp::Int64 = rand(1:11) # pick α hypothesis ("what's wrong with the patient")
    right = rand(Uniform(1, 10)) # effect of right intervention
    wrong = rand(Uniform(1, 10)) # effect of wrong intervention

    data::Array{Bool, 1} = datFunc((rand_hyp - 1) / (numb_hyp - 1)) # generate synthetic data for this pick (the test results for the patient)

    updates = Array{Float32,2}(undef, numb_toss + 1, numb_hyp) # initialize array for probabilities

    updates[1, :] = prior # set prior

    if rule_index == 1.0f0
        @fastmath @inbounds for t in 1:numb_toss # generate updates
            updates[t + 1, :] = b_upd(updates[t, :], data, t)
        end
    elseif rule_index == 2.0f0
        @fastmath @inbounds for t in 1:numb_toss # generate updates
            updates[t + 1, :] = expl_upd(updates[t, :], data, t, c_value)
        end
    elseif rule_index == 3.0f0
        @fastmath @inbounds for t in 1:numb_toss # generate updates
            updates[t + 1, :] = good_upd(updates[t, :], data, t, c_value)
        end
    else
        @fastmath @inbounds for t in 1:numb_toss # generate updates
            updates[t + 1, :] = pop_upd(updates[t, :], data, t, c_value)
        end
    end

    return dist(updates, rand_hyp, right, wrong, thresh)
end

#= tests doctor on 100 patients and calculates average survival score obtained by doctor, so
average probability that patients would survive =#
@everywhere function avScore(rule_index::Float32, c_value::Float32, thresh::Float32, dist::Function)
    tot = @distributed (+) for i in 1:100
        patient(rule_index, c_value, thresh, dist)
    end
    return tot / 100
end

#= tests all doctors in a population and selects the 50 percent fittest; outputs those
as well as a copy of each of them =#
function population_upd_rep(pop::Array{Float32,2}, dist::Function)
    agent_scores = SharedArray{Float32,1}(numb_agents)
    @sync @distributed for i in 1:numb_agents
        @inbounds agent_scores[i] = avScore(pop[i, :]..., dist)
    end
    best_index = findall(agent_scores .>= Statistics.median(agent_scores))
    best = pop[best_index[1:Int(numb_agents/2)], :]
    return vcat(best, best), agent_scores
end

run(`mkdir data`)

function sim_run(dist::Function)
    k = 1
    while k < numb_simulations + 1
        groupID = repeat(1.0:4.0, inner=div(numb_agents, 4))
        population_start = vcat(fill(0, div(numb_agents, 4)), rand(Uniform(0, .25), 3*div(numb_agents, 4)))
        agent_thresh = rand(Uniform(.5, 1), numb_agents)
        pop_start = Array{Float32,2}(hcat(groupID, population_start, agent_thresh))
        pop_upd_c_a = Array{Float32,3}(undef, numb_agents, 3, numb_generations + 1)
        pop_upd_f = Array{Float32,2}(undef, numb_agents, numb_generations + 1)
        pop_upd_c_a[:, :, 1] = pop_start

        @inbounds for i in 1:numb_generations
            pop_upd_c_a[:, :, i + 1], pop_upd_f[:, i] = population_upd_rep(pop_upd_c_a[:, :, i], dist)
        end

        @inbounds for i in 1:numb_agents
            pop_upd_f[i, numb_generations + 1] = avScore(pop_upd_c_a[i, 1:3, numb_generations + 1]..., dist)
        end

        res_a = Array{Int32,2}(undef, numb_agents, numb_generations + 1)
        res_c = Array{Float32,2}(undef, numb_agents, numb_generations + 1)
        res_t = Array{Float32,2}(undef, numb_agents, numb_generations + 1)

        @inbounds for i in 1:(numb_generations + 1)
            res_a[:, i], res_c[:, i], res_t[:, i] = pop_upd_c_a[:, 1, i], pop_upd_c_a[:, 2, i], pop_upd_c_a[:, 3, i]
        end

        writedlm("data/type$dist$k.txt", res_a)
        writedlm("data/c_value$dist$k.txt", res_c)
        writedlm("data/thresh$dist$k.txt", res_t)
        writedlm("data/fit$dist$k.txt", pop_upd_f)

        population_start = nothing
        pop_start = nothing
        pop_upd_c_a = nothing
        pop_upd_f = nothing
        res_a = nothing
        res_c = nothing
        GC.gc()

        k += 1
    end
end

sim_run(survWei)

fullType = readdlm("data/typesurvWei1.txt")

ks = [ keys(sort(countmap(fullType[:,i]))) for i in 1:numb_generations + 1 ]
vls = [ values(sort(countmap(fullType[:,i]))) for i in 1:numb_generations + 1 ]

group = []
freq = []
gen = []

for i in 1:numb_generations + 1
    append!(group, collect(ks[i]))
    append!(freq, collect(vls[i]))
    append!(gen, fill(i, length(collect(ks[i]))))
end

bar_df = DataFrame(group=group, freq=freq, gen=gen)

bar_df[!, :name] = map(bar_df[!, :group]) do x
    if x == 1
        return "Bayes"
    elseif x == 2
        return "EXPL"
    elseif x == 3
        return "Good"
    else
        return "Popper"
    end
end

rename!(bar_df, [Symbol("$i") for i in ["Group", "Count", "Generation", "Rule"]])

p1 = plot(bar_df, x=:Generation, y=:Count, color=:Rule, Geom.bar(position=:stack),
    Coord.cartesian(xmin=1, xmax=numb_generations + 1),
    Guide.colorkey(title="Rule"),
    Scale.color_discrete_manual(gen_brew_colors(4)...),
    style(minor_label_font_size=10pt, major_label_font_size=14pt,
            key_label_font_size=11pt, key_title_font_size=13pt,
            grid_color=colorant"#222831",
            colorkey_swatch_shape=:square))

p1 |> PDF("exa_with_var_thresh.pdf")

full_results_type = Array{Int32, 3}(undef, numb_agents, numb_generations + 1, numb_simulations)

for i in 1:numb_simulations
    full_results_type[:, :, i] = readdlm("nthresh_data/typesurvWei$i.txt")
end

percentages = Array{Float64, 3}(undef, numb_simulations, numb_generations + 1, 4)

for k in 1:4
    for j in 1:numb_generations + 1
        percentages[:, j, k] = [ length(findall(full_results_type[:, j, i] .== k)) / numb_agents for i in 1:numb_simulations ]
    end
end

@everywhere function bs(x, n=1000)
    bootstrap(mean, x, BasicSampling(n))
end

bayesBoot = SharedArray{Float64}(numb_generations + 1, 3);

@distributed for i in 1:numb_generations + 1
    bayesBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 1]), BasicConfInt(.95))[1]...]
end

explBoot = SharedArray{Float64}(numb_generations + 1, 3);

@distributed for i in 1:numb_generations + 1
    explBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 2]), BasicConfInt(.95))[1]...]
end

goodBoot = SharedArray{Float64}(numb_generations + 1, 3);

@distributed for i in 1:numb_generations + 1
    goodBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 3]), BasicConfInt(.95))[1]...]
end

popperBoot = SharedArray{Float64}(numb_generations + 1, 3);

@distributed for i in 1:numb_generations + 1
    popperBoot[i, :] = [Bootstrap.confint(bs(percentages[:, i, 4]), BasicConfInt(.95))[1]...]
end

Generation = repeat(collect(1:numb_generations + 1), outer=4)
Rule = repeat(["Bayes", "EXPL", "Good", "Popper"], inner=numb_generations + 1)
res = vcat(bayesBoot, explBoot, goodBoot, popperBoot)
type_df = DataFrame(res, :auto)
type_df.Generation=Generation
type_df.Rule=Rule
rename!(type_df, [Symbol("$i") for i in ["y", "ymin", "ymax", "Generation", "Rule"]]);

p2 = plot(type_df, x=:Generation, y=:y, ymin=:ymin, ymax=:ymax, color=:Rule, Geom.line, Geom.ribbon,
    Guide.ylabel("Average percentage"),
    Coord.cartesian(xmin=1, xmax=numb_generations + 1, ymin=-.001),
    Scale.color_discrete_manual(gen_brew_colors(4)...),
    style(line_width=2pt, lowlight_color=c->RGB{Float32}(c.r, c.g, c.b), alphas=[.2,.2,.2,.2],
          minor_label_font_size=10pt, major_label_font_size=14pt,
          key_label_font_size=11pt, key_title_font_size=13pt))

p2 |> PDF("perc_with_var_thresh.pdf")

p = plot(type_df, x=:Generation, y=:y, color=:Rule, Geom.line,
    Guide.ylabel("Average percentage"),
    Coord.cartesian(xmin=1, xmax=numb_generations + 1, ymin=-.001),
    Scale.color_discrete_manual(gen_brew_colors(4)...),
    style(line_width=2pt,
            minor_label_font_size=10pt, major_label_font_size=14pt,
            key_label_font_size=11pt, key_title_font_size=13pt,
            colorkey_swatch_shape=:square))

p |> PDF("percentages_thresh.pdf")

full_results_t = Array{Float64, 3}(undef, numb_agents, numb_generations + 1, numb_simulations)

for i in 1:numb_simulations
    full_results_t[:, :, i] = readdlm("nthresh_data/threshsurvWei$i.txt")
end

thresh_res = Array{Float64, 2}(undef, numb_generations + 1, 3)

for k in 1:numb_generations + 1
    cr = Float64[]
    for j in 1:numb_simulations
        push!(cr, mean([ full_results_t[:, k, i] for i in 1:numb_simulations ][j]))
    end
    cr = cr[.!isnan.(cr)]
    thresh_res[k, :] = vcat(mean(cr), [confint(OneSampleTTest(cr))...])
end

mns = mean_and_std.([ full_results_t[:, i, :][:] for i in 1:numb_generations + 1 ])

thresh_df = DataFrame(mean=first.(mns), std=last.(mns))
thresh_df.gen = 1:size(thresh_df, 1)
thresh_df.min = thresh_df.mean .- thresh_df.std
thresh_df.max= thresh_df.mean .+ thresh_df.std

q = plot(thresh_df, x=:gen, y=:mean, ymin=:min, ymax=:max, Geom.line, Geom.ribbon,
         #Coord.cartesian(xmin=1, xmax=25),
         Scale.x_log2,
         Guide.xlabel("Generation"), Guide.ylabel("Average threshold value"),
         Theme(line_width=2pt, default_color = RGB(102/255,194/255,165/255),
               minor_label_font_size=10pt, major_label_font_size=14pt, alphas=[.2],
               key_label_font_size=11pt, key_title_font_size=13pt))

q |> PDF("thresholds.pdf")

# thresholds per agent type
bayes_thresh = [ mean_and_std(vcat([ full_results_t[:, j, i][full_results_type[:, j, i] .== 1] for i in 1:numb_simulations ]...)) for j in 1:numb_generations + 1 ]
expl_thresh = [ mean_and_std(vcat([ full_results_t[:, j, i][full_results_type[:, j, i] .== 2] for i in 1:numb_simulations ]...)) for j in 1:numb_generations + 1 ]
good_thresh = [ mean_and_std(vcat([ full_results_t[:, j, i][full_results_type[:, j, i] .== 3] for i in 1:numb_simulations ]...)) for j in 1:numb_generations + 1 ]
popper_thresh = [ mean_and_std(vcat([ full_results_t[:, j, i][full_results_type[:, j, i] .== 4] for i in 1:numb_simulations ]...)) for j in 1:numb_generations + 1 ]

Generation = repeat(collect(1:numb_generations + 1), outer=4)
Rule = repeat(["Bayes", "EXPL", "Good", "Popper"], inner=numb_generations + 1)
t_res = vcat(bayes_thresh, expl_thresh, good_thresh, popper_thresh)
t_df = DataFrame(y = first.(t_res), ymin = first.(t_res) .- last.(t_res), ymax = first.(t_res) .+ last.(t_res), Generation = Generation, Rule = Rule)

p_thresh = plot(t_df, x=:Generation, y=:y, ymin=:ymin, ymax=:ymax, color=:Rule, Geom.line, Geom.ribbon,
    Guide.ylabel("Average threshold value"),
    #Coord.cartesian(xmin=1, xmax=numb_generations + 1, ymin=-.001),
    Scale.x_log2,
    Scale.color_discrete_manual(gen_brew_colors(4)...),
    style(line_width=2pt, lowlight_color=c->RGB{Float32}(c.r, c.g, c.b), alphas=[.2,.2,.2,.2],
          minor_label_font_size=10pt, major_label_font_size=14pt,
          key_label_font_size=11pt, key_title_font_size=13pt))

p_thresh |> PDF("thresh_per_type.pdf")


text:
@otestset "SurfaceDefs" begin
    @testset "Knots" begin
        # find span
        knots = KnotVector{Int64}([0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5])
        function apply(knots::KnotVector, curveorder, u, correctindex)
            knotnum = findspan(knots, curveorder, u)
            return knotnum == correctindex
        end
        @test apply(knots, 2, 2.5, 5) && apply(knots, 2, 0, 3) && apply(knots, 2, 4, 6) && apply(knots, 2, 5, 8)

        # insert knots
        knots = [0, 0, 0, 0, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7, 7, 7]
        insertions = knotstoinsert(knots, 3)
        @test insertions == [(5, 2), (6, 2), (7, 1), (9, 1), (11, 1), (13, 2)]
    end

    @testset "BSpline" begin
        # bspline conversion
        correctverts = [[0.0, 0.0, 0.0], [0.0, 0.49625, 0.0], [0.6625, 0.49625, 0.65625], [0.0, 0.0, 0.0], [0.6625, 0.49625, 0.65625], [0.6625, 0.0, 0.0], [0.6625, 0.0, 0.0], [0.6625, 0.49625, 0.65625], [1.33, 0.49625, 0.0], [0.6625, 0.0, 0.0], [1.33, 0.49625, 0.0], [1.33, 0.0, 0.0], [0.0, 0.49625, 0.0], [0.0, 1.0, 0.0], [0.6625, 1.0, 0.0], [0.0, 0.49625, 0.0], [0.6625, 1.0, 0.0], [0.6625, 0.49625, 0.65625], [0.6625, 0.49625, 0.65625], [0.6625, 1.0, 0.0], [1.33, 1.0, 0.0], [0.6625, 0.49625, 0.65625], [1.33, 1.0, 0.0], [1.33, 0.49625, 0.0]]
        correcttris = [
            0x00000001 0x00000002 0x00000003
            0x00000004 0x00000005 0x00000006
            0x00000007 0x00000008 0x00000009
            0x0000000a 0x0000000b 0x0000000c
            0x0000000d 0x0000000e 0x0000000f
            0x00000010 0x00000011 0x00000012
            0x00000013 0x00000014 0x00000015
            0x00000016 0x00000017 0x00000018
        ]
        surf = TestData.bsplinesurface()
        points, triangles = makiemesh(makemesh(surf, 2))
        segments = tobeziersegments(surf) # TODO just testing that this evaluates for now
        @test triangles == correcttris
        @test all(isapprox.(points, correctverts, rtol = RTOLERANCE, atol = ATOLERANCE))
    end # testset BSpline

    @testset "PowerBasis" begin
        # polynomial is 3x^2 + 2x + 1
        coeff = [1.0 2.0 3.0]
        curve = PowerBasisCurve{OpticSim.Euclidean,Float64,1,2}(coeff)
        @test !all(isapprox.(coeff, coefficients(curve, 1), rtol = RTOLERANCE, atol = ATOLERANCE)) # TODO not sure if this is correct/what this is testing?
        correct = true
        for x in -10.0:0.1:10
            exact = 3 * x^2 + 2 * x + 1
            calculated = point(curve, x)[1]
            @test isapprox(exact, calculated, rtol = RTOLERANCE, atol = ATOLERANCE)
        end
    end # testset PowerBasis

    @testset "BezierSurface" begin
        surf = TestData.beziersurface()
        fdm = central_fdm(10, 1)
        for u in 0:0.1:1, v in 0:0.1:1
            (du, dv) = partials(surf, u, v)
            fu(u) = point(surf, u, v)
            fv(v) = point(surf, u, v)
            # compute accurate finite difference approximation to derivative
            fdu, fdv = (fdm(fu, u), fdm(fv, v))
            @test isapprox(du, fdu, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dv, fdv, rtol = 1e-12, atol = 2 * eps(Float64))
        end
    end # testset BezierSurface

    @testset "ZernikeSurface" begin
        surf = TestData.zernikesurface1a() # with normradius
        fdm = central_fdm(10, 1)
        for ρ in 0.05:0.1:0.95, ϕ in 0:(π / 10):(2π)
            (dρ, dϕ) = partials(surf, ρ, ϕ)
            fρ(ρ) = point(surf, ρ, ϕ)
            fϕ(ϕ) = point(surf, ρ, ϕ)
            # compute accurate finite difference approximation to derivative
            fdρ, fdϕ = (fdm(fρ, ρ), fdm(fϕ, ϕ))
            @test isapprox(dρ, fdρ, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dϕ, fdϕ, rtol = 1e-12, atol = 2 * eps(Float64))
        end

        @test !any(isnan.(normal(TestData.conicsurface(), 0.0, 0.0)))
    end # testset ZernikeSurface

    @testset "QTypeSurface" begin
        # test predefined values against papers
        # Forbes, G. W. "Robust, efficient computational methods for axially symmetric optical aspheres." OpticSim express 18.19 (2010): 19700-19712.
        # Forbes, G. W. "Characterizing the shape of freeform optics." OpticSim express 20.3 (2012): 2483-2499.

        f0_true = [2, sqrt(19 / 4), 4 * sqrt(10 / 19), 1 / 2 * sqrt(509 / 10), 6 * sqrt(259 / 509), 1 / 2 * sqrt(25607 / 259)]
        for i in 1:length(f0_true)
            @test isapprox(OpticSim.QType.f0(i - 1), f0_true[i], rtol = 2 * eps(Float64))
        end
        g0_true = [-1 / 2, -5 / (2 * sqrt(19)), -17 / (2 * sqrt(190)), -91 / (2 * sqrt(5090)), -473 / (2 * sqrt(131831))]
        for i in 1:length(g0_true)
            @test isapprox(OpticSim.QType.g0(i - 1), g0_true[i], rtol = 2 * eps(Float64))
        end
        h0_true = [-1 / 2, -6 / sqrt(19), -3 / 2 * sqrt(19 / 10), -20 * sqrt(10 / 509)]
        for i in 1:length(h0_true)
            @test isapprox(OpticSim.QType.h0(i - 1), h0_true[i], rtol = 2 * eps(Float64))
        end

        F_true = [1/4 1/2 27/32 5/4; 15/32 7/8 35/16 511/128; 17/72 35/36 35/16 23/6; 29/40 67/40 243/80 12287/2560]
        N, M = size(F_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.F(m, n), F_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        G_true = [1/4 3/8 15/32 35/64; -1/24 -5/48 -7/32 -21/64; -7/40 -7/16 -117/160 -33/32]
        N, M = size(G_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.G(m, n), G_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        f_true = [1/2 1/sqrt(2) 3 / 4*sqrt(3 / 2) sqrt(5)/2; 1 / 4*sqrt(7 / 2) 1 / 4*sqrt(19 / 2) 1 / 4*sqrt(185 / 6) 3 / 32*sqrt(427); 1 / 6*sqrt(115 / 14) 1 / 2*sqrt(145 / 38) 1 / 4*sqrt(12803 / 370) 1 / 2*sqrt(2785 / 183); 1 / 5*sqrt(3397 / 230) 1 / 4*sqrt(6841 / 290) 9 / 4*sqrt(14113 / 25606) 1 / 16*sqrt(1289057 / 1114)]
        N, M = size(f_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.f(m, n), f_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        g_true = [1/2 3/(4 * sqrt(2)) 5/(4 * sqrt(6)) 7 * sqrt(5)/32; -1/(3 * sqrt(14)) -5/(6 * sqrt(38)) -7 / 4*sqrt(3 / 370) -1 / 2*sqrt(7 / 61); -21 / 10*sqrt(7 / 230) -7 / 4*sqrt(19 / 290) -117 / 4*sqrt(37 / 128030) -33 / 16*sqrt(183 / 2785)]
        N, M = size(g_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.g(m, n), g_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        A_true = [2 3 5 7 9 11; -4/3 2/3 2 76/27 85/24 106/25; 9/5 26/15 2 117/50 203/75 108/35; 55/28 66/35 2 536/245 135/56 130/49; 161/81 122/63 2 515/243 143/63 161/66]
        N, M = size(A_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.A(m, n), A_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        B_true = [-1 -2 -4 -6 -8 -10; -8/3 -4 -4 -40/9 -5 -28/5; -24/5 -4 -4 -21/5 -112/25 -24/5; -30/7 -4 -4 -144/35 -30/7 -220/49; -112/27 -4 -4 -110/27 -88/21 -13/3]
        N, M = size(B_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.B(m, n), B_true[n + 1, m], rtol = 2 * eps(Float64))
            end
        end
        C_true = [NaN NaN NaN NaN NaN NaN; -11/3 -3 -5/3 -35/27 -9/8 -77/75; 0 5/9 7/15 21/50 88/225 13/35; 27/28 21/25 27/35 891/1225 39/56 33/49; 80/81 45/49 55/63 1430/1701 40/49 1105/1386]
        N, M = size(C_true)
        for m in 1:M
            for n in 0:(N - 1)
                @test isapprox(OpticSim.QType.C(m, n), C_true[n + 1, m], rtol = 2 * eps(Float64)) || (isnan(OpticSim.QType.C(m, n)) && isnan(C_true[n + 1, m]))
            end
        end

        # test deriv with finite differences
        surf = TestData.qtypesurface1()
        fdm = central_fdm(10, 1)
        for ρ in 0.05:0.1:0.95, ϕ in 0:(π / 10):(2π)
            (dρ, dϕ) = partials(surf, ρ, ϕ)
            fρ(ρ) = point(surf, ρ, ϕ)
            fϕ(ϕ) = point(surf, ρ, ϕ)
            # compute accurate finite difference approximation to derivative
            fdρ, fdϕ = (fdm(fρ, ρ), fdm(fϕ, ϕ))
            @test isapprox(dρ, fdρ, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dϕ, fdϕ, rtol = 1e-12, atol = 2 * eps(Float64))
        end
    end # testset QTypeSurface

    @testset "GridSag" begin
        surf = TestData.gridsagsurfacelinear()
        fdm = central_fdm(10, 1)
        # doesn't enforce C1 across patch boundaries meaning that finite differences won't match at all
        # just test within a patch to make sure partials() is working
        for ρ in 0.02:0.01:0.18, ϕ in 0:(π / 30):(2π)
            (dρ, dϕ) = partials(surf, ρ, ϕ)
            fρ(ρ) = point(surf, ρ, ϕ)
            fϕ(ϕ) = point(surf, ρ, ϕ)
            # compute accurate finite difference approximation to derivative
            fdρ, fdϕ = (fdm(fρ, ρ), fdm(fϕ, ϕ))
            @test isapprox(dρ, fdρ, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dϕ, fdϕ, rtol = 1e-12, atol = 2 * eps(Float64))
        end

        surf = TestData.gridsagsurfacebicubic()
        fdm = central_fdm(10, 1)
        # doesn't enforce C2 across patch boundaries meaning that finite differences won't match exactly
        # just test within a patch to make sure partials() is working
        for ρ in 0.02:0.01:0.18, ϕ in 0:(π / 30):(2π)
            (dρ, dϕ) = partials(surf, ρ, ϕ)
            fρ(ρ) = point(surf, ρ, ϕ)
            fϕ(ϕ) = point(surf, ρ, ϕ)
            # compute accurate finite difference approximation to derivative
            fdρ, fdϕ = (fdm(fρ, ρ), fdm(fϕ, ϕ))
            @test isapprox(dρ, fdρ, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dϕ, fdϕ, rtol = 1e-12, atol = 2 * eps(Float64))
        end

        surf = TestData.gridsagsurfacebicubiccheby()
        fdm = central_fdm(10, 1)
        # not valid at the very boundary of the surface as stuff gets weird outside |u| < 1 and |v| < 1
        for u in -0.3:0.01:0.3, v in -0.15:0.01:0.15
            (du, dv) = partials(surf, u, v)
            fu(u) = point(surf, u, v)
            fv(v) = point(surf, u, v)
            # compute accurate finite difference approximation to derivative
            fdu, fdv = (fdm(fu, u), fdm(fv, v))
            @test isapprox(du, fdu, rtol = 1e-12, atol = 2 * eps(Float64))
            @test isapprox(dv, fdv, rtol = 1e-12, atol = 2 * eps(Float64))
        end
    end # testset GridSag

    @testset "ChebyshevSurface" begin
        surf = TestData.chebyshevsurface1() # with normradius
        fdm = central_fdm(10, 1)
        # not valid at the very boundary of the surface as stuff gets weird outside |u| < 1 and |v| < 1
        for u in -0.99:0.1:0.99, v in -0.99:0.1:0.99
            (du, dv) = partials(surf, u, v)
            fu(u) = point(surf, u, v)
            fv(v) = point(surf, u, v)
            # compute accurate finite difference approximation to derivative
            fdu, fdv = (fdm(fu, u), fdm(fv, v))
            @test isapprox(du, fdu, rtol = 1e-11, atol = 2 * eps(Float64)) #changed rtol from 1e-12 to 1e-11. FiniteDifferences approximation to the derivative had larger than expected error.
            @test isapprox(dv, fdv, rtol = 1e-11, atol = 2 * eps(Float64)) #changed rtol from 1e-12 to 1e-11. FiniteDifferences approximation to the derivative had larger than expected error.
        end
    end # testset ChebyshevSurface

end # testset SurfaceDefs

ext: jl


text:
using Term
import Term.Renderables: Renderable
import Term: int
import Term: typestree, inspect
using Term.Layout

import MyterialColors: orange_light, blue_light, green_light

"""
Convert HSL values to RGB values, return a markup color string.
"""
function hsl2rgb(h, s, l)  # pragma: no cover
    C = (1 - abs(2 * l - 1)) * s
    X = C * (1 - abs(mod(h / 60, 2) - 1))
    M = l - C / 2

    if 0 ≤ h < 60
        r, g, b = C, X, 0
    elseif 60 ≤ h < 120
        r, g, b = X, C, 0
    elseif 120 ≤ h < 180
        r, g, b = 0, C, X
    elseif 180 ≤ h < 240
        r, g, b = 0, X, C
    elseif 240 ≤ h < 300
        r, g, b = X, 0, C
    elseif 300 ≤ h ≤ 360
        r, g, b = C, 0, X
    end

    r = (Int ∘ round)((r + M) * 255)
    g = (Int ∘ round)((g + M) * 255)
    b = (Int ∘ round)((b + M) * 255)

    return "($r, $g, $b)"
end

function make_julia_circles()
    circle = """
    oooo
 oooooooooo
oooooooooooo
oooooooooooo
 oooooooooo
    oooo    """

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

    hspacer = Spacer(green.measure.w / 2 + 1, green.measure.h)
    line = Spacer(green.measure.w * 2 + 6, 1)
    circles =
        line / (hspacer * green * hspacer) / (red * Spacer(2, purple.measure.h) * purple)
    return circles
end

function make_rgb_colors(; max_width = 88)
    colors = ""

    for y in 0:5
        for x in 0:max_width
            h = x / max_width
            l = 0.1 + (((5 - y) / 5) * 0.7)
            color = hsl2rgb(h * 360, 0.9, l)
            bg = hsl2rgb(h * 360, 0.9, l + 0.7 / 10)

            colors *= "{$color on_$bg}▄{/$color on_$bg}"
            # colors *= "{$color on_$bg}▬▄x{/$color on_$bg}"
        end
        colors *= "\n"
    end

    return colors
end

function rainbow_maker(N)
    R = range(30, 200; length = N)
    G = range(200, 60; length = N)
    B = range(50, 200; length = N)
    colors = []
    for n in 1:N
        r, g, b = int(R[n]), int(G[n]), int(B[n])
        push!(colors, "($r,$g,$b)")
    end
    return colors
end

circles = make_julia_circles() # 42 x 17

_code_style = "yellow italic bold"

bfc = rainbow_maker(9)
basic_features = Panel(
    """
    {bright_red bold underline}Features{/bright_red bold underline}

{bold $(bfc[1])}✔{/bold $(bfc[1])}{white} {blue}Colored text{/blue}{/white}
{bold $(bfc[2])}✔{/bold $(bfc[2])}{white} {italic}italic{/italic}, {bold}bold{/bold}, {underline}underline{/underline}, {striked}striked{/striked}, {inverse}inverse{/inverse}{/white}
{bold $(bfc[3])}✔{/bold $(bfc[3])}{white} styling {$_code_style}@macros{/$_code_style}{/white}
{bold $(bfc[4])}✔{/bold $(bfc[4])}{white} {italic white}markup{/italic white} style syntax{/white}
{bold $(bfc[5])}✔{/bold $(bfc[5])}{white} progress bars{/white}
{bold $(bfc[6])}✔{/bold $(bfc[6])}{white} {$_code_style}`Expr`{/$_code_style} and {$_code_style}`Type`{/$_code_style} introspection{/white}
{bold $(bfc[7])}✔{/bold $(bfc[7])}{white} logging{/white}
{bold $(bfc[8])}✔{/bold $(bfc[8])}{white} stacktraces{/white}
{bold $(bfc[9])}✔{/bold $(bfc[9])}{white} syntax highlighting{/white}
""";
    width = 70,
    padding = (2, 2, 1, 2),
    justify = :center,
    style = "default blue dim",
    title = "Term.jl",
    title_style = "default bright_blue bold",
    subtitle = "https://github.com/FedeClaudi/Term.jl",
    subtitle_style = "default dim",
    subtitle_justify = :right,
)

colors_info = TextBox(
    """{bold bright_green}Colors{/bold bright_green}
      {$(green_light) bold}✔{/$(green_light) bold} 8-bit colors
      {$(green_light) bold}✔{/$(green_light) bold} 16-bit colors
      {$(green_light) bold}✔{/$(green_light) bold} hex colors
      {$(green_light) bold}✔{/$(green_light) bold} rgb colors
      {$(green_light) bold}✔{/$(green_light) bold} colors conversion
    """;
    width = 25,
    padding = (0, 0, 0, 0),
)
colors = make_rgb_colors(; max_width = 103)

_lorem = """朗眠裕安無際集正聞進士健音社野件草売規作独特認権価官家複入豚末告設悟自職遠氷育教載最週場仕踪持白炎組特曲強真雅立覧自価宰身訴側善論住理案者券真犯著避銀楽験館稿告
"""
lorem_description = TextBox(
    """{bold bright_yellow}Text reshaping{/bold bright_yellow}
      {$(orange_light) bold}✔{/$(orange_light) bold} reshaping.
      {$(orange_light) bold}✔{/$(orange_light) bold} justification
      {$(orange_light) bold}✔{/$(orange_light) bold} Asian languages support
    """;
    width = 30,
    padding = (0, 0, 0, 0),
)
lorem1 = TextBox(_lorem; width = 62, fit = false, padding = (0, 0, 0, 0))
lorem2 = TextBox(_lorem; width = 42, fit = false, padding = (0, 0, 0, 0))

expr = :(2x + 2π / θ)
tree = Renderable(sprint(typestree, Float64))
dendo = Renderable(sprint(inspect, expr))
renderables_info = TextBox(
    """{bold bright_blue}Renderables types{/bold bright_blue}
      {bold $(blue_light)}✔{/bold $(blue_light)} Panel
      {bold $(blue_light)}✔{/bold $(blue_light)} TextBox
      {bold $(blue_light)}✔{/bold $(blue_light)} hLine
      {bold $(blue_light)}✔{/bold $(blue_light)} vLine
      {bold $(blue_light)}✔{/bold $(blue_light)} Tree
      {bold $(blue_light)}✔{/bold $(blue_light)} Dendogram

    """;
    width = 30,
    padding = (0, 0, 0, 0),
)

line = hLine(140; style = "bold dim grey35", box = :HEAVY)

layout_text = TextBox(
    """{bold bright_cyan}Renderables layout{/bold bright_cyan}
      {bold bright_cyan}✔{/bold bright_cyan} Horizontal stacking
      {bold bright_cyan}✔{/bold bright_cyan} Vertical stacking
      {bold bright_cyan}✔{/bold bright_cyan} Left justify
      {bold bright_cyan}✔{/bold bright_cyan} Center justify
      {bold bright_cyan}✔{/bold bright_cyan} Right justify
      {bold bright_cyan}✔{/bold bright_cyan} Shorthands {bold bright_blue}*,/{/bold bright_blue}
    """;
    width = 30,
    padding = (0, 0, 0, 0),
)

p1 = Panel(; width = 20, style = "#80bbe8")
p1b = Panel(; width = 34, style = "#80bbe8")
p2 = Panel(; width = 28, style = "#5692bf")
p3 = Panel(; width = 34, style = "#316c99 bold")
p3b = Panel(; width = 20, style = "#316c99 bold")
_space = Spacer(1, layout_text.measure.h)

layout_example =
    hstack(lvstack(p1, p2, p3), cvstack(p1b, p2, p3b), rvstack(p1, p2, p3); pad = 1)

layout_example = cvstack(
    "" / RenderableText(
        "{bold #81bae6} Left/center/right justify and stack renderables to create layouts",
    ),
    layout_example,
)

# ----------------------------------- print ---------------------------------- #
# width is 140
print("\n"^3)

readme =
    (
        Spacer(10, circles.measure.h) *
        circles *
        Spacer(8, circles.measure.h) *
        basic_features
    ) / (Spacer(140, 2) / line) / (
        Spacer(3, colors_info.measure.h) *
        colors_info *
        Spacer(3, colors_info.measure.h) *
        (Spacer(100, 1) / colors)
    ) / line /
    (Spacer(3, lorem_description.measure.h) * lorem_description * lorem1 * lorem2) / line /
    (
        Spacer(3, tree.measure.h) *
        renderables_info *
        tree *
        Spacer(5, tree.measure.h) *
        dendo
    ) / line / (Spacer(3, layout_text.measure.h) * layout_text * layout_example)

print(Spacer(10, readme.measure.h) * readme)


text:
include("hrb_utils.jl")

type HBMeta
  # Metadata attached to a Harwell-Boeing matrix.
  title :: AbstractString
  key   :: AbstractString

  totcrd :: Int
  ptrcrd :: Int
  indcrd :: Int
  valcrd :: Int
  rhscrd :: Int

  mxtype :: AbstractString
  nrow :: Int
  ncol :: Int
  nnzero :: Int
  neltvl :: Int

  hermitian :: Bool
  assembled :: Bool

  ptrfmt :: AbstractString
  indfmt :: AbstractString
  valfmt :: AbstractString
  rhsfmt :: AbstractString

  nrhs   :: Int
  rhstyp :: AbstractString
  nrhsix :: Int
end

@compat RHSType = Union{Array{Float64,2}, Array{Complex64,2}, SparseMatrixCSC}

type HarwellBoeingMatrix

  meta   :: HBMeta
  matrix :: SparseMatrixCSC
  rhs    :: RHSType  # Right-hand side, if any.
  guess  :: RHSType  # Initial guesses, if any.
  sol    :: RHSType  # Solutions, if any.

  function HarwellBoeingMatrix(file_name :: AbstractString)
    hb = open(file_name)

    # Read header.
    seekstart(hb)
    line  = readline(hb)
    title = strip(line[1:72])      # A72
    key   = strip(line[73:end-1])  # A8
    totcrd, ptrcrd, indcrd, valcrd, rhscrd = map(s -> parse(Int, s),
                                                 split(chomp(readline(hb))))

    line = readline(hb)
    mxtype = line[1:3]   # A3
    # Skip 11 blanks.
    nrow, ncol, nnzero, neltvl = map(s -> parse(Int, s),
                                     [line[14 + 14*(i-1)+1 : 14 + 14*i] for i = 1 : 4])

    pattern_only = (mxtype[1] == 'P')
    is_complex = (mxtype[1] == 'C')
    hermitian = (mxtype[2] == 'S')
    assembled = (mxtype[3] == 'A')
    data_type = is_complex ? Complex64 : Float64

    line = readline(hb)
    ptrfmt = strip(line[1:16])  # 2A16
    indfmt = strip(line[17:32])

    valfmt = pattern_only ? "" : strip(line[33:52])  # 2A20
    rhsfmt = (pattern_only | rhscrd == 0) ? "" : strip(line[53:72])

    # Read right-hand side info, if present.
    if rhscrd > 0
      line = readline(hb)
      rhstyp = line[1:3]  # A3
      nrhs, nrhsix = map(s -> parse(Int, s), split(chomp(line[15:end])))

    else
      rhstyp = "   "
      nrhs = 0
      nrhsix = 0
    end

    # Read integer data.
    ip = read_array(hb, ncol+1, ptrfmt)
    ind = read_array(hb, nnzero, indfmt)

    # Read values, if supplied.
    if pattern_only
      vals = ones(Float64, nnzero)  # Just so we can build a SparseMatrixCSC.
    else
      vallen = assembled ? nnzero : neltvl
      vals = read_array(hb, vallen, valfmt, is_complex=is_complex)
    end
    # Ensure row indices are sorted in each column.
    sortsparse!(ip, ind, vals)
    matrix = SparseMatrixCSC(nrow, ncol, ip, ind, vals)

    # Read right-hand sides, if any.
    if (nrhs > 0) & (rhstyp[1] != 'F') & assembled  # Sparse rhs.
        rhsptr  = read_array(hb, nrhs+1, ptrfmt)
        rhsind  = read_array(hb, nrhsix, indfmt)
        rhsvals = read_array(hb, nrhsix, rhsfmt, is_complex=is_complex)
        # Ensure row indices are sorted in each column.
        sortsparse!(rhsptr, rhsind, rhsvals)
        rhs = SparseMatrixCSC(nrow, nrhs, rhsptr, rhsind, rhsvals)

    else
      lenrhs = rhstyp[1] == 'F' ? nrow : nnzero
      rhs    = zeros(data_type, (nrow, nrhs))

      for i = 1 : nrhs
        rhs[:,i] = read_array(hb, lenrhs, rhsfmt, is_complex=is_complex)
      end
    end

    guess  = zeros(data_type, (nrow, rhstyp[2] == 'G' ? nrhs : 0))
    sol    = zeros(data_type, (nrow, rhstyp[3] == 'X' ? nrhs : 0))

    # Read initial guesses, if any.
    if rhstyp[2] == 'G'
      for i = 1 : nrhs
        guess[:,i] = read_array(hb, nrow, rhsfmt, is_complex=is_complex)
      end
    end

    # Read solution vectors, if any.
    if rhstyp[3] == 'X'
      for i = 1 : nrhs
        sol[:,i] = read_array(hb, nrow, rhsfmt, is_complex=is_complex)
      end
    end

    close(hb)

    meta = HBMeta(title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd,
                  mxtype, nrow, ncol, nnzero, neltvl, hermitian, assembled,
                  ptrfmt, indfmt, valfmt, rhsfmt, nrhs, rhstyp, nrhsix)

    new(meta, matrix, rhs, guess, sol)
  end
end

# Displaying Harwell-Boeing matrix instances.

import Base.show, Base.print
function show(io :: IO, hb :: HarwellBoeingMatrix)
  s  = @sprintf("Harwell-Boeing matrix %s of type %s\n", hb.meta.key, hb.meta.mxtype)
  s *= @sprintf("%d rows, %d cols, %d nonzeros\n", hb.meta.nrow, hb.meta.ncol, hb.meta.nnzero)
  s *= @sprintf("%d right-hand sides, %d guesses, %d solutions\n",
                size(hb.rhs,2), size(hb.guess,2), size(hb.sol,2))
  print(io, s)
end

function print(io :: IO, hb :: HarwellBoeingMatrix)
  @printf("Harwell-Boeing matrix %s\n", hb.meta.key)
  @printf("%s\n", hb.meta.title)
  @printf("%d rows, %d cols, %d nonzeros\n", hb.meta.nrow, hb.meta.ncol, hb.meta.nnzero)
  if hb.meta.mxtype[1] == 'P'
    dtype = "pattern only"
  elseif hb.meta.mxtype[1] == 'R'
    dtype = "real"
  else
    dtype = "complex"
  end
  herm = hb.meta.hermitian ? "hermitian" : "non-hermitian"
  assm = hb.meta.assembled ? "assembled" : "elemental"
  @printf("(%s, %s, %s)\n", dtype, herm, assm)
  @printf("%d right-hand sides, %d guesses, %d solutions\n",
           size(hb.rhs,2), size(hb.guess,2), size(hb.sol,2))
  if size(hb.rhs,2) > 0
    @printf("Right-hand side(s):\n");
    display(hb.rhs'); @printf("\n")
  end
  if size(hb.guess,2) > 0
    @printf("Initial guess(es):\n");
    display(hb.guess'); @printf("\n")
  end
  if size(hb.sol,2) > 0
    @printf("Solution(s):\n");
    display(hb.sol'); @printf("\n")
  end
end


text:
module Io
  # this module contains input/output functions

  # Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
  #
  #  MIT license.
  #
  # Permission is hereby granted, free of charge, to any person
  # obtaining a copy of this software and associated documentation files (the
  # "Software"), to deal in the Software without restriction, including without
  # limitation the rights to use, copy, modify, merge, publish, distribute,
  # sublicense, and/or sell copies of the Software, and to permit persons to whom
  # the Software is furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included in all
  # copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  # SOFTWARE.
  export load

  import MAT
  import FixedPointNumbers

  #wrapper MATLAB/Octave load function based on MAT.jl library,
  #for mat files returns the closes analogue to a struct array - Dictionary type
  function load(matFile, treatAs = "-mat", variable = "all")
    if treatAs == "-mat"
      #load all variables as a Dictionary type - the closes analogue to struct array
      allData = MAT.matread(matFile)

      if variable == "all"
        matFileDict = allData
      elseif variable != "all"
        if variable in keys(allData)
          matFileDict = allData[variable]
        else
          error("specified variable doesn't exist in the mat file")
        end
      end

    elseif treatAs == "-ascii"
      # load the ascii file and return an array of Float64 type
      fileContent = open(readlines, matFile)
      matFileDict = map(parse, fileContent)
    else
      error("unknown type specified in the second argument")
    end

    return matFileDict
  end
end #End of Io


text:
using ProbNumDiffEq
using Random
using Statistics
using Plots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax);

u0 = [990.0,10.0,0.0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

sol_ode_ek0 = solve(prob_ode,
                EK0(prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true),
                dt=δt,
                abstol=1e-1,
                reltol=1e-2);

sol_ode_ek1 = solve(prob_ode,
                EK1(prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true),
                dt=δt,
                abstol=1e-1,
                reltol=1e-2);

s20_ek0 = sol_ode_ek0(20.0)
[mean(s20_ek0) std(s20_ek0)]

s20_ek1 = sol_ode_ek1(20.0)
[mean(s20_ek1) std(s20_ek1)]

num_samples = 100
samples_ode_ek0 = ProbNumDiffEq.sample(sol_ode_ek0, num_samples);
samples_ode_ek1 = ProbNumDiffEq.sample(sol_ode_ek1, num_samples);

p_ek0 = plot(sol_ode_ek0.t,
         samples_ode_ek0[:, :, 1],
         label=["S" "I" "R"],
         color=[:blue :red :green],
         xlabel="Time",
         ylabel="Number",
         title="EK0")
for i in 2:num_samples
    plot!(p_ek0,
          sol_ode_ek0.t,
          samples_ode_ek0[:, :, i],
          label="",
          color=[:blue :red :green])
end;

p_ek1 = plot(sol_ode_ek1.t,
         samples_ode_ek1[:, :, 1],
         label=["S" "I" "R"],
         color=[:blue :green],
         xlabel="Time",
         ylabel="Number",
         title="EK1")
for i in 2:num_samples
    plot!(p_ek1,
          sol_ode_ek1.t,
          samples_ode_ek1[:, :, i],
          label="",
          color=[:blue :red :green],)
end;

plot(p_ek0, p_ek1, layout = (1,2), xlim=(15,20),ylim=(100,1000),yaxis=:log10)

plot(p_ek0, p_ek1, layout = (1,2), xlim=(35,40),ylim=(10,1000),yaxis=:log10)

@benchmark solve(prob_ode,
                 EK0(prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true),
                 abstol=1e-1,
                 reltol=1e-2)

@benchmark solve(prob_ode,
                 EK1(prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true),
                 abstol=1e-1,
                 reltol=1e-2)


text:
# # CO2 and temperature

using AbstractGPs
using AbstractGPsMakie
using CairoMakie
using CSV
using DataDeps
using DataFrames
using Dates
using Distributions
using KernelFunctions
using LinearAlgebra
using Optim
using ParameterHandling
using Random
using Statistics
using Stheno
using Zygote

using CairoMakie: RGB

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

# ## Obtain the data

# Use DataDeps.jl to obtain the data.
register(DataDep(
    "mauna_loa",
    "The classic Mauna Loa CO2 data.",
    "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/monthly/" *
        "monthly_in_situ_co2_mlo.csv",
))
register(DataDep(
    "hadcrut5_global",
    "Globally-averaged surface temperatures",
    "https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/analysis/diagnostics/" *
        "HadCRUT.5.0.1.0.analysis.summary_series.global.monthly.csv",
))

# Load and pre-process the data.
ml_df, Ttr_df, Tte_df = let

    mauna_loa_co2 = let
        mauna_loa_data_raw = CSV.read(
            joinpath(datadep"mauna_loa", "monthly_in_situ_co2_mlo.csv"), DataFrame;
            skipto=58, header=false,
        )

        data = DataFrame(
            date=map(Date, mauna_loa_data_raw.Column1, mauna_loa_data_raw.Column2),
            CO2=map(x -> x == -99.99 ? missing : x, mauna_loa_data_raw.Column5),
        )
        data = filter(:CO2 => !ismissing, data)

        ## standardise CO2 data to have unit mean.
        m = mean(data.CO2)
        σ = std(data.CO2)
        data.CO2 = Float64.(data.CO2 .- m) ./ σ

        ## Convert dates into floats.
        data.float_date = map(x -> x.value, data.date .- Date(1)) ./ 365

        data
    end

    temperature = let
        hadcrut5_data_raw = CSV.read(
            joinpath(
                datadep"hadcrut5_global",
                "HadCRUT.5.0.1.0.analysis.summary_series.global.monthly.csv",
            ),
            DataFrame;
            header=1,
        )

        tmp = DataFrame(
            date=hadcrut5_data_raw[!, "Time"],
            temperature=hadcrut5_data_raw[!, "Anomaly (deg C)"],
        )

        ## Convert dates into floats.
        tmp.float_date = map(x -> x.value, tmp.date .- Date(1)) ./ 365

        ## Remove any data from before the first Mauna Loa data.
        time_min = minimum(mauna_loa_co2.date)
        filter(:date => date -> date >= time_min, tmp)
    end

    ## Split temperature data into train and test.
    Ntr = 425
    Ttr_raw = temperature[1:Ntr, :]
    Tte_raw = temperature[(Ntr + 1):end, :]

    ## Standard temperature data using training data.
    mtr = mean(Ttr_raw.temperature)
    σtr = std(Ttr_raw.temperature)
    Ttr = DataFrame(
        date=Ttr_raw.date,
        float_date=Ttr_raw.float_date,
        temperature=(Ttr_raw.temperature .- mtr) ./ σtr,
    )
    Tte = DataFrame(
        date=Tte_raw.date,
        float_date=Tte_raw.float_date,
        temperature=(Tte_raw.temperature .- mtr) ./ σtr,
    )

    mauna_loa_co2, Ttr, Tte
end;

# ## Build and train a model for the data

# Model for CO2 and temperature as a function of time.
function build_gp(θ)
    return @gppp let

        ## Shared trend process.
        f_trend = stretch(GP(SEKernel()), θ.λ_trend)

        ## Specify model for CO2.
        f_co2_latent = θ.CO2.σ_latent * f_trend
        f_co2_wiggle = θ.CO2.σ_wiggle * stretch(GP(SEKernel()), θ.CO2.λ_wiggle)
        f_co2_period = θ.CO2.σ_period * GP(SEKernel() ∘ PeriodicTransform(θ.CO2.freq))
        f_co2 = f_co2_latent + f_co2_wiggle + f_co2_period + θ.CO2.σm * GP(ConstantKernel())

        ## Specify model for temperature.
        f_T_trend = θ.T.σ_trend * f_trend
        f_T_wiggle = θ.T.σ_wiggle * stretch(GP(SEKernel()), θ.T.λ_wiggle)
        f_T = f_T_trend + f_T_wiggle + θ.T.σm * GP(ConstantKernel())
    end
end

x_co2 = ml_df.float_date;
x_T_tr = Ttr_df.float_date;
x_T_te = Tte_df.float_date;

y_co2 = ml_df.CO2;
y_T_tr = Ttr_df.temperature;
y_T_te = Tte_df.temperature;

x = BlockData(GPPPInput(:f_co2, x_co2), GPPPInput(:f_T, x_T_tr))
y = vcat(y_co2, y_T_tr)

function build_obs_cov(θ)
    σ²_CO2 = θ.CO2.σ_noise^2
    σ²_T = θ.T.σ_noise^2
    return Diagonal(vcat(fill(σ²_CO2, length(y_co2)), fill(σ²_T, length(y_T_tr))))
end

function nlml(θ)
    f = build_gp(θ)
    return -logpdf(f(x, build_obs_cov(θ)), y)
end

init_params = (
    λ_trend = positive(1e-1),
    CO2 = (
        σ_latent = positive(0.3),
        σ_wiggle = positive(0.3),
        λ_wiggle = positive(5.0),
        σ_period = positive(0.3),
        freq = fixed(1.0),
        σ_noise = positive(0.1),
        σm = positive(1.0),
    ),
    T = (
        σ_trend = positive(1.0),
        σ_wiggle = positive(1.0),
        λ_wiggle = positive(0.5),
        σ_noise = positive(0.1),
        σm = positive(1.0),
    ),
)

default_optimizer = LBFGS(;
    alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
    linesearch=Optim.LineSearches.BackTracking(),
)

function optimize_loss(loss, θ_init; optimizer=default_optimizer, maxiter=1_000)
    options = Optim.Options(; iterations=maxiter, show_trace=true)

    θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init)
    loss_packed = loss ∘ unflatten

    ## https://julianlsolvers.github.io/Optim.jl/stable/#user/tipsandtricks/#avoid-repeating-computations
    function fg!(F, G, x)
        if F !== nothing && G !== nothing
            val, grad = Zygote.withgradient(loss_packed, x)
            G .= only(grad)
            return val
        elseif G !== nothing
            grad = Zygote.gradient(loss_packed, x)
            G .= only(grad)
            return nothing
        elseif F !== nothing
            return loss_packed(x)
        end
    end

    result = optimize(Optim.only_fg!(fg!), θ_flat_init, optimizer, options; inplace=false)

    return unflatten(result.minimizer), result
end

θ_opt, result = optimize_loss(nlml, init_params)

# ## Plot the resulting model fit.

function colours()
    return Dict(
        :blue => RGB(0/255, 107/255, 164/255),
        :cyan => RGB(75/255, 166/255, 251/255),
        :red => RGB(200/255, 82 / 255, 0 / 255),
        :pink => RGB(169/255, 90/255, 161/255),
        :black => RGB(0.0, 0.0, 0.0),
        :orange => RGB(245/255, 121/255, 58/255),
    )
end

band_alpha() = 0.3
sample_alpha() = 0.2
point_alpha() = 1.0

x_plot = range(
    minimum(vcat(ml_df.float_date, Ttr_df.float_date)),
    maximum(vcat(ml_df.float_date, Tte_df.float_date));
    length=1_500,
);

f_post = posterior(build_gp(θ_opt)(x, build_obs_cov(θ_opt)), y);

function plot_band!(ax, x_plot, fx, colour, label)
    ms = marginals(fx)
    symband!(
        ax, x_plot, mean.(ms), std.(ms);
        bandscale=3,
        color=(colours()[colour], 0.5 * band_alpha()),
        label=label,
    )
end

function plot_sample!(ax, x_plot, fx, colour)
    gpsample!(
        ax, x_plot, fx;
        samples=4, color=(colours()[colour], sample_alpha()),
    )
end

function plot_gp!(ax, x_plot, fx, colour, label)
    plot_band!(ax, x_plot, fx, colour, label)
    plot_sample!(ax, x_plot, fx, colour)
end

let
    fig = Figure()
    ga = fig[1, 1] = GridLayout()

    jitter = 1e-6

    ## Plot CO2 posterior.
    let
        ax = Axis(ga[1, 1])
        f_post_x = f_post(GPPPInput(:f_co2, x_plot), jitter)

        plot_band!(ax, x_plot, f_post(GPPPInput(:f_co2_latent, x_plot), jitter), :orange, "CO2 (latent)")
        plot_band!(ax, x_plot, f_post(GPPPInput(:f_co2_wiggle, x_plot), jitter), :black, "CO2 (wiggle)")
        plot_band!(ax, x_plot, f_post(GPPPInput(:f_co2_period, x_plot), jitter), :pink, "CO2 (period)")

        plot_gp!(ax, x_plot, f_post_x, :blue, "CO2")
        scatter!(ax, x_co2, y_co2; color=(colours()[:blue], point_alpha()), markersize=3)

        axislegend(ax; position=:lt)

        hidexdecorations!(ax; grid=false)
    end

    ## Plot temperature posterior.
    let
        ax = Axis(ga[2, 1])
        f_post_x = f_post(GPPPInput(:f_T, x_plot), jitter)

        plot_band!(ax, x_plot, f_post(GPPPInput(:f_T_trend, x_plot), jitter), :orange, "T (trend)")
        plot_band!(ax, x_plot, f_post(GPPPInput(:f_T_wiggle, x_plot), jitter), :black, "T (wiggle)")

        plot_gp!(ax, x_plot, f_post_x, :blue, "T")
        scatter!(ax, x_T_tr, y_T_tr; color=(colours()[:blue], point_alpha()), markersize=3)
        scatter!(
            ax, x_T_te, y_T_te;
            color=(colours()[:pink], point_alpha()),
            markersize=5,
            marker=:xcross,
            label="T (test)"
        )

        axislegend(ax; position=:lt)
    end

    rowgap!(ga, 5)

    fig
end


text:
using Distributions
using Base.Test

# Test for gradlogpdf on univariate distributions

@test_approx_eq_eps gradlogpdf(Beta(1.5, 3.0), 0.7) -5.9523809523809526 1.0e-8
@test_approx_eq_eps gradlogpdf(Chi(5.0), 5.5) -4.7727272727272725 1.0e-8
@test_approx_eq_eps gradlogpdf(Chisq(7.0), 12.0) -0.29166666666666663 1.0e-8
@test_approx_eq_eps gradlogpdf(Exponential(2.0), 7.0) -0.5 1.0e-8
@test_approx_eq_eps gradlogpdf(Gamma(9.0, 0.5), 11.0) -1.2727272727272727 1.0e-8
@test_approx_eq_eps gradlogpdf(Gumbel(3.5, 1.0), 4.0) -1.6065306597126334 1.0e-8
@test_approx_eq_eps gradlogpdf(Laplace(7.0), 34.0) -1.0 1.0e-8
@test_approx_eq_eps gradlogpdf(Logistic(-6.0), 1.0) -0.9981778976111987 1.0e-8
@test_approx_eq_eps gradlogpdf(LogNormal(5.5), 2.0) 1.9034264097200273 1.0e-8
@test_approx_eq_eps gradlogpdf(Normal(-4.5, 2.0), 1.6) -1.525 1.0e-8
@test_approx_eq_eps gradlogpdf(TDist(8.0), 9.1) -0.9018830525272548 1.0e-8
@test_approx_eq_eps gradlogpdf(Weibull(2.0), 3.5) -6.714285714285714 1.0e-8

# Test for gradlogpdf on multivariate distributions

@test_approx_eq_eps gradlogpdf(MvNormal([1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9]) [0.191919191919192, 1.080808080808081] 1.0e-8
@test_approx_eq_eps gradlogpdf(MvTDist(5., [1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9]) [0.2150711513583442, 1.2111901681759383] 1.0e-8


text:
function sort_by_bits(nums)
    #result = repeat([Array{Integer,1}()],16)
    result = [Array{Integer,1}() for _=1:16]
    for num in nums
        let temp = result[count_ones(num) + 1]
            result[count_ones(num) + 1] = push!(temp, num)
        end
    end
    collect(Iterators.flatten(map(x -> sort(x), filter(!isempty, result))))
end

using Test
@test sort_by_bits([0,1,2,3,4,5,6,7,8]) == [0,1,2,4,8,3,5,6,7]
@test sort_by_bits([2,3,5,7,11,13,17,19]) == [2,3,5,17,7,11,13,19]
@test sort_by_bits([10,100,1000,10000]) == [10,100,10000,1000]
@test sort_by_bits([1024,512,256,128,64,32,16,8,4,2,1]) == [1,2,4,8,16,32,64,128,256,512,1024]

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
export Rocchio, RocchioBagging, fit, predict, transform
import SimilaritySearch: KnnResult
using StatsBase: countmap

mutable struct Rocchio
    protos::Vector{BOW}  # prototypes
    pops::Vector{Int} # population per class
end

const RocchioBagging = Vector{Rocchio}

function fit(::Type{Rocchio}, X::AbstractVector{BOW}, y::AbstractVector{Int}; nclasses=0)
    if nclasses == 0
        nclasses = unique(y) |> length
    end

    prototypes = [BOW() for i in 1:nclasses]
    populations = zeros(Int, nclasses)
    println(stderr, "fitting Rocchio classifier with $(length(X)) items; and $nclasses classes")
    for i in 1:length(X)
        i % 1000 == 0 && print(stderr, "*")
        i % 100000 == 0 && println(stderr, " adv: $(i / length(X))")
        c = y[i]
        add!(prototypes[c], X[i])
        populations[c] += 1
    end

    for bow in prototypes
        normalize!(bow)
    end

    Rocchio(prototypes, populations)
end

function predict(rocchio::Rocchio, x::BOW)
    res = KnnResult(1)
    for i in 1:length(rocchio.protos)
        d = cosine_distance(rocchio.protos[i], x)
        push!(res, i, d)
    end

    first(res).objID
end

function transform(rocchio::Rocchio, x::BOW)
    [dot(rocchio.protos[i], x) for i in 1:length(rocchio.protos)]
end

function broadcastable(rocchio::Rocchio)
    (rocchio,)
end

function fit(::Type{RocchioBagging}, X::AbstractVector{BOW}, y::AbstractVector{Int}, nrocchios=5; nclasses=0)::RocchioBagging
    if nclasses == 0
        nclasses = unique(y) |> length
    end

    s = ceil(Int, length(X) / 2)
    L = RocchioBagging()
    for i in 1:nrocchios
        I = rand(1:length(X), s)
        rocchio = fit(Rocchio, X[I], y[I], nclasses=nclasses)
        push!(L, rocchio)
    end

    L
end

function predict(rocchio::RocchioBagging, x::BOW)
    P = [predict(r, x) for r in rocchio] |> countmap |> collect
    sort!(P, by=x->x[end])
    P[end][1]
end

function transform(rocchio::RocchioBagging, x::BOW)
    sum(transform(r, x) for r in rocchio) |> normalize!
end

function broadcastable(rocchio::RocchioBagging)
    (rocchio,)
end


text:
@testset "generic.jl" begin
    @testset "Simple example" begin
        # sample predictions and outcomes
        nsamples = 1_000
        predictions = rand(nsamples)
        outcomes = rand(Bool, nsamples)

        # create bin with all predictions and outcomes
        bin = CalibrationErrors.Bin(predictions, outcomes)

        # check statistics
        @test bin.nsamples == nsamples
        @test bin.mean_predictions ≈ mean(predictions)
        @test bin.proportions_targets == mean(outcomes)

        # compare with adding data
        bin2 = CalibrationErrors.Bin(predictions[1], outcomes[1])
        for i in 2:nsamples
            CalibrationErrors.adddata!(bin2, predictions[i], outcomes[i])
        end
        @test bin2.nsamples == bin.nsamples
        @test bin2.mean_predictions ≈ bin.mean_predictions
        @test bin2.proportions_targets ≈ bin.proportions_targets
    end

    @testset "Simple example ($nclasses classes)" for nclasses in (2, 10, 100)
        # sample predictions and targets
        nsamples = 1_000
        dist = Dirichlet(nclasses, 1.0)
        predictions = [rand(dist) for _ in 1:nsamples]
        targets = rand(1:nclasses, nsamples)

        # create bin with all predictions and targets
        bin = CalibrationErrors.Bin(predictions, targets)

        # check statistics
        @test bin.nsamples == nsamples
        @test bin.mean_predictions ≈ mean(predictions)
        @test bin.proportions_targets == proportions(targets, nclasses)
        @test sum(bin.mean_predictions) ≈ 1
        @test sum(bin.proportions_targets) ≈ 1

        # compare with adding data
        bin2 = CalibrationErrors.Bin(predictions[1], targets[1])
        for i in 2:nsamples
            CalibrationErrors.adddata!(bin2, predictions[i], targets[i])
        end
        @test bin2.nsamples == bin.nsamples
        @test bin2.mean_predictions ≈ bin.mean_predictions
        @test bin2.proportions_targets ≈ bin.proportions_targets
    end
end


text:
module OCReract

export
    run_tesseract

include("tesseract.jl")

"""
OCReract is a simple Julia wrapper of the well-known OCR engine called Tesseract.

Here, a simple example of usage:

# Example
```julia-repl
julia> using Images
julia> using OCReract
julia> img_path = "/path/to/img.png";
# In disk
julia> run_tesseract(img_path, "/tmp/res.txt", psm=3, oem=1)
# In memory
julia> img = Images.load(img_path);
julia> res_text = run_tesseract(img, psm=3, oem=1);
julia> println(strip(res_text));
```

For more information, check the homepage in https://github.com/leferrad/OCReract.jl.
"""
OCReract

end # module


text:
module EMISAgentSimulation

#################################################################################
# Exports

# Export Structs
export CaseDefinition

export AgentSimulation

export Investor

export KalmanFilter
export InvestorBelief
export Scenario
export Forecast
export Perfect
export Imperfect

export RiskPreference
export RiskNeutral
export RiskAverse

export MarketPrices

export Finance

export Product
export OperatingProduct
export ReserveDirection
export ReserveUpEMIS
export ReserveDownEMIS
export ReserveUpEMIS
export OperatingReserve
export Energy
export Capacity
export REC

export BuildPhase
export Existing
export Planned
export Queue
export Option
export Retired

export DeviceEMIS
export Project
export GeneratorEMIS
export StorageEMIS

export BatteryTech
export BatteryEMIS

export HydroTech
export HydroGenEMIS

export RenewableTech
export RenewableGenEMIS

export ThermalTech
export ThermalGenEMIS

export ZonalLine

export MarketProject
export EnergyMarket
export ReserveUpMarket
export ReserveDownMarket
export CapacityMarket
export RECMarket
export MarketCollection
export MarketClearingProblem

# Export Functions

# Export Core functions
export calculate_annualized_costs
export calculate_capcost_and_wacc
export calculate_derating_data
export calculate_npv
export capacity_market_clearing
export cem
export cluster
export create_agent_simulation
export create_capacity_demand_curve
export create_expected_marketdata
export create_PSY_generator
export create_realized_marketdata
export economicdispatch
export energy_mkt_clearing
export finish_construction!
export make_investments!
export rec_market_clearing
export retire_old!
export retire_unprofitable!
export run_agent_simulation
export start_construction!

# Export Utility functions
export calculate_operating_profit
export calculate_required_processes
export create_parallel_workers
export dir_exists
export leaftypes
export make_case_data_dir
export make_results_dir
export read_data
export remove_leap_day!
export size_in_MW
export update_derating_factor!
export update_installed_cap!
export update_PSY_timeseries!
export update_simulation_derating_data!
export write_data

# Export Getter Functions
export get_accepted_perc
export get_active_power_limit
export get_active_power_limits
export get_activeprojects
export get_allprojects
export get_all_techs
export get_annual_cashflow
export get_annual_growth
export get_base_dir
export get_belief_update
export get_bus
export get_capacity_forward_years
export get_capacity_price
export get_capacity_bid
export get_capacity_factors
export get_cap_cost_multiplier
export get_capex_years
export get_case
export get_constant
export get_construction_year
export get_data_dir
export get_decidedprojects
export get_decision_year
export get_derating
export get_derating_data
export get_device_size
export get_discount_rate
export get_effective_investment_cost
export get_efficiency
export get_end_life_year
export get_energy_price
export get_error_covariance_estimate
export get_existing
export get_expected_npv
export get_expected_utility
export get_finance_data
export get_fixed_OM_cost
export get_forecast
export get_forecast_type
export get_from_zone
export get_fuel
export get_heterogeneity
export get_hour_weight
export get_info_symmetry
export get_input_active_power_limits
export get_investment_cost
export get_investor_data
export get_investors
export get_iteration_year
export get_kalman_filter
export get_lag_time
export get_life_time
export get_line_rating
export get_lines
export get_marginal_cost
export get_market_prices
export get_markets
export get_max_annual_projects
export get_maxcap
export get_max_limit
export get_maxperc
export get_measurement_covariance
export get_mincap
export get_multiplier
export get_name
export get_num_rep_days
export get_operation_cost
export get_options
export get_output_active_power_limits
export get_ownedby
export get_parallel_investors
export get_parallel_scenarios
export get_parameter_multipliers
export get_parameter_values
export get_peak_load
export get_planned
export get_prices
export get_probability
export get_process_covariance
export get_products
export get_projects
export get_queue
export get_queue_cost
export get_ramp_limits
export get_realized_profit
export get_rec_certificates
export get_rec_price
export get_rec_bid
export get_rep_hour_weight
export get_reservedown_price
export get_reserveup_price
export get_results_dir
export get_retired
export get_retirement_year
export get_risk_aversion
export get_risk_coefficient
export get_risk_preference
export get_rolling_horizon
export get_scenario_npv
export get_scenario_profit
export get_scenario_utility
export get_scenario_data
export get_siip_market_clearing
export get_simulation_years
export get_soc
export get_solver
export get_start_year
export get_state_estimate
export get_state_measurement_matrix
export get_state_transition_matrix
export get_storage_capacity
export get_system_ED
export get_system_UC
export get_system_services
export get_tech
export get_sys_dir
export get_time_limits
export get_to_zone
export get_total_horizon
export get_type
export get_uncertainty
export get_zone
export get_zones
export update_belief!

# Export Find functions
export find_active_invested_projects
export find_energy_product
export find_operating_products
export find_option_projects

# Export Setter Functions
export set_investors!
export set_iteration_year!
export set_peak_load!

#################################################################################
# Imports
import AxisArrays
import Clustering
import ClusterManagers
import CSV
import DataFrames
import Dates
import Distributed
import FileIO
import InteractiveUtils
import JuMP
export optimizer_with_attributes
import JLD2
import LinearAlgebra
import PowerSystems
import PowerSimulations
import InfrastructureSystems

const PSY = PowerSystems
const PSI = PowerSimulations
const IS = InfrastructureSystems

import Random
import UUIDs
import StatsBase
import Statistics
import TimeSeries

const TS = TimeSeries

import Base.convert
import JuMP.value
import MathOptInterface
import MathOptInterface: AbstractOptimizer

const MOI = MathOptInterface

################################################################################
# Includes

# Include all structs:

include("structs/products/Product.jl")
include("structs/products/Energy.jl")
include("structs/products/OperatingReserve.jl")
include("structs/products/Capacity.jl")
include("structs/products/REC.jl")

include("structs/Finance.jl")

include("structs/devices/BuildPhase.jl")
include("structs/devices/Project.jl")
include("structs/devices/BatteryEMIS.jl")
include("structs/devices/HydroGenEMIS.jl")
include("structs/devices/RenewableGenEMIS.jl")
include("structs/devices/ThermalGenEMIS.jl")

include("structs/CaseDefinition.jl")
include("structs/MarketPrices.jl")
include("structs/Scenario.jl")

include("structs/InvestorBelief.jl")

include("structs/Forecast.jl")
include("structs/RiskPreference.jl")
include("structs/ZonalLine.jl")

include("structs/Investor.jl")

include("structs/AgentSimulation.jl")

include("structs/market_structs/MarketProject.jl")
include("structs/market_structs/EnergyMarket.jl")
include("structs/market_structs/ReserveUpMarket.jl")
include("structs/market_structs/ReserveDownMarket.jl")
include("structs/market_structs/CapacityMarket.jl")
include("structs/market_structs/RECMarket.jl")
include("structs/market_structs/MarketCollection.jl")
include("structs/market_structs/MarketClearingProblem.jl")

# Include utilities:
include("utils/general.jl")                 # General utilities.
include("utils/market_utils.jl")           # Utils for actual and expected market clearing modules.
include("utils/product_utils.jl")          # Utils for products provided by the projects
include("utils/siip_psy_utils.jl")         # Utils for SIIP PowerSystems.
include("utils/conversion_utils.jl")       # Define new convert functions for changing buildphase of projects.
include("utils/read_and_write_utils.jl")   # Read and write utils.
include("utils/parallel_utils.jl")         # Utils for parallelizing price prediction runs.
include("utils/finance_utils.jl")          # Functions for calculating adjusted CAPEX and WACC

#Include files containing functions for creating the simulation structs from the given data.
include("struct_creators/simulation_structs/product_creator.jl")
include("struct_creators/simulation_structs/project_creator.jl")
include("struct_creators/simulation_structs/project_vector_creator.jl")
include("struct_creators/simulation_structs/investor_creator.jl")
include("struct_creators/simulation_structs/rts_psy_creator.jl")
include("struct_creators/simulation_structs/siip_psy_device_creator.jl")
include("struct_creators/simulation_structs/agent_simulation_creator.jl")

#Include files containing functions for creating the structs used for market clearing modules.
include("struct_creators/market_structs/market_project_creator.jl")
include("struct_creators/market_structs/cem_creator.jl")
include("struct_creators/market_structs/economic_dispatch_creator.jl")

#Include Test System Parsers
include("test_system_parsers/test_system_reader.jl")
include("test_system_parsers/rts_reader.jl")

#Include Expected and Actual Market Simultion functions.
include("markets_simulation/cem.jl")
include("markets_simulation/expected_market_simulation.jl")
include("markets_simulation/economicdispatch.jl")
include("markets_simulation/siip_simulation_definition.jl")
include("markets_simulation/actual_energy_mkt_clearing.jl")
include("markets_simulation/actual_capacity_mkt_clearing.jl")
include("markets_simulation/actual_rec_mkt_clearing.jl")
include("markets_simulation/actual_market_simulation.jl")

#Include Investor functions
    include("investor_functions/investor_iteration.jl")            # Runs investors annual iteration.

    #### Predictions #####################
    include("investor_functions/prediction/prediction_methodology.jl")  # Functions for running investors' price prediction methodology.
    include("investor_functions/prediction/operating_profit.jl")        # Functions for calculating expected operating market profits.
    include("investor_functions/prediction/capacity_profit.jl")         # Functions for calculating expected capacity market profits.
    include("investor_functions/prediction/REC_profit.jl")              # Functions for calculating expected REC market profits.
    include("investor_functions/prediction/total_profit.jl")            # Functions for updating the expected profits from all markets.

    #### Decisions ######################
    include("investor_functions/decisions/buildphase_conversion.jl")         # Functions for evaluating when to convert the buildphase of projects.
    include("investor_functions/decisions/investment_decision.jl")           # Functions for making investment decisions
    include("investor_functions/decisions/retirement_decision.jl")           # Functions for making retirement decisions

    #### Decision Metrics ######
    include("investor_functions/decision_metrics/npv_functions.jl")          # Functions for calculating NPV.
    include("investor_functions/decision_metrics/utility_functions.jl")      # Functions for calculating expected utility.

    ### Realized Profits and Updates ######
    include("investor_functions/realized_profits_calculator.jl")   # Functions for calculating realized profits from different markets
    include("investor_functions/annual_updates.jl")                # Functions for updating investor revenues and forecasts each year.

#Include derating factor updating methodology
include("derating_factor_updates/derating_factor_calculator.jl")  # Derating factor calculation for renewables

#Include k-medoids clustering for representative days selection.
include("representative_days/clustering.jl")

#Include main investment simulation function.
include("investment_simulation.jl")

end


text:
using Random, Test
rng = Random.MersenneTwister(12345)

function random_psd(n)

    A = sprandn(rng,n,n,0.2)
    A = A+A';
    A = A + Diagonal(((sum(abs.(A),dims=1)[:]))) #make diagonally dominant

end

@testset "All Unit Tests" begin

  include("./UnitTests/basic.jl")
  include("./UnitTests/refactoring.jl")
  include("./UnitTests/non-quasidef.jl")
  include("./UnitTests/inertia.jl")
  include("./UnitTests/regularization.jl")
  include("./UnitTests/update_values.jl")

end
nothing


text:
using HDF5, DataFrames, CSV, JLD
using StatsBase, Polynomials, CurveFit

kep_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv")
gaia_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "gaiadr2_keplerdr25_crossref.csv")
mast_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "KeplerMAST_TargetProperties.csv")
stellar_catalog_file_out = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1q17_dr25_gaia_fgk.jld")

kep_df = CSV.read(kep_filename)
gaia_df = CSV.read(gaia_filename, rows_for_type_detect=30000)
mast_df = CSV.read(mast_filename)

dup_gaiaid = find(nonunique(DataFrame(x = gaia_df[:source_id])))
deleterows!(gaia_df, dup_gaiaid)

println("Total crossref target stars = ", length(gaia_df[:kepid]))

mag_diff = gaia_df[:phot_g_mean_mag].-gaia_df[:kepmag]
quant_arr = quantile(mag_diff, [0.067,0.933])   # 1.5-sigma cut
mag_match = find(x->quant_arr[1]<=x<=quant_arr[2], mag_diff)
gaia_df = gaia_df[mag_match,:]

gaia_col = [:kepid,:source_id,:parallax,:parallax_error,:astrometric_gof_al,:astrometric_excess_noise_sig,:bp_rp,:priam_flags,:teff_val,:teff_percentile_lower,:teff_percentile_upper,:radius_val,:radius_percentile_lower,:radius_percentile_upper,:lum_val,:lum_percentile_lower,:lum_percentile_upper]
df = join(kep_df, gaia_df[:,gaia_col], on=:kepid)
df = join(df, mast_df, on=:kepid)
kep_df = nothing
gaia_df = nothing

df = join(df, mast_df, on=:kepid)
kep_df = nothing
gaia_df = nothing

println("Total target stars (KOIs) matching magnitude = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

df[:teff] = df[:teff_val]
df[:teff_err1] = df[:teff_percentile_upper].-df[:teff_val]
df[:teff_err2] = df[:teff_percentile_lower].-df[:teff_val]
delete!(df, :teff_val)
delete!(df, :teff_percentile_upper)
delete!(df, :teff_percentile_lower)
df[:radius_err1] = df[:radius_percentile_upper].-df[:radius_val]
df[:radius_err2] = df[:radius_percentile_lower].-df[:radius_val]
df[:radius] = df[:radius_val]
delete!(df, :radius_val)
delete!(df, :radius_percentile_upper)
delete!(df, :radius_percentile_lower)

not_binary = (df[:astrometric_gof_al] .<= 20) .& (df[:astrometric_excess_noise_sig] .<= 5)
astro_good = []
for x in 1:length(df[:kepid])
    if !(ismissing(df[x,:priam_flags]))
        pflag = string(df[x,:priam_flags])
         if (pflag[2] == '0') & (pflag[3] == '0')
             push!(astro_good, true)
         else
             push!(astro_good, false)
         end
     else
         push!(astro_good, false)
     end
end
#astro_good = astro_good .& (df[:parallax_error] .< 0.05*df[:parallax])
astro_good = astro_good .& (df[:parallax_error] .< 0.1*df[:parallax])
planet_search = df[:kepmag] .<= 16.

has_mass = .! (ismissing.(df[:mass]) .| ismissing.(df[:mass_err1]) .| ismissing.(df[:mass_err2]))
has_radius = .! (ismissing.(df[:radius]) .| ismissing.(df[:radius_err1]) .| ismissing.(df[:radius_err2]))
has_dens = .! (ismissing.(df[:dens]) .| ismissing.(df[:dens_err1]) .| ismissing.(df[:dens_err2]))

has_cdpp = .! (ismissing.(df[:rrmscdpp01p5]) .| ismissing.(df[:rrmscdpp02p0]) .| ismissing.(df[:rrmscdpp02p5]) .| ismissing.(df[:rrmscdpp03p0]) .| ismissing.(df[:rrmscdpp03p5]) .| ismissing.(df[:rrmscdpp04p5]) .| ismissing.(df[:rrmscdpp05p0]) .| ismissing.(df[:rrmscdpp06p0]) .| ismissing.(df[:rrmscdpp07p5]) .| ismissing.(df[:rrmscdpp09p0]) .| ismissing.(df[:rrmscdpp10p5]) .| ismissing.(df[:rrmscdpp12p0]) .| ismissing.(df[:rrmscdpp12p5]) .| ismissing.(df[:rrmscdpp15p0]))
has_rest = .! (ismissing.(df[:dataspan]) .| ismissing.(df[:dutycycle]))
#has_limbd = .! (ismissing.(df[:limbdark_coeff1]) .| ismissing.(df[:limbdark_coeff2]) .| ismissing.(df[:limbdark_coeff3]) .| ismissing.(df[:limbdark_coeff4]))

mast_cut =.&(df[:numLCEXqtrs].>0,df[:numLCqtrs].>4)

is_usable = .&(has_radius, has_mass, has_rest, has_cdpp, has_dens, astro_good, not_binary, planet_search, mast_cut)

df = df[find(is_usable),:]
println("Total stars (KOIs) with valid parameters = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

fgk_color = (0.5 .<= df[:bp_rp] .<= 1.7)
df = df[find(fgk_color),:]
is_FGK = []
coeff = [2.5,-3.6,0.9]
for i in 1:6
    is_FGK = (log10.(df[:lum_val]).< map(x->polyval(Poly(coeff),x),df[:bp_rp]) + log10(1.75))
    coeff = poly_fit(df[is_FGK,:bp_rp],log10.(df[is_FGK,:lum_val]),2)
end
FGK = find(is_FGK)

println("Total FGK stars (KOIs) with valid parameters = ", length(FGK), " (", sum(df[:nkoi]),")")

# See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
# TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
symbols_to_keep = [ :kepid, :source_id, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp01p5, :rrmscdpp02p0, :rrmscdpp02p5, :rrmscdpp03p0, :rrmscdpp03p5, :rrmscdpp04p5, :rrmscdpp05p0, :rrmscdpp06p0, :rrmscdpp07p5, :rrmscdpp09p0, :rrmscdpp10p5, :rrmscdpp12p0, :rrmscdpp12p5, :rrmscdpp15p0, :dataspan, :dutycycle ]
delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
df = df[FGK, symbols_to_keep]
tmp_df = DataFrame()
for col in names(df)
    tmp_df[col] = collect(skipmissing(df[col]))
end
df = tmp_df

save(stellar_catalog_file_out,"stellar_catalog", df)


text:
using MacroTools

function Base.convert(::Type{Expression}, ex::Expr)
    ex.head === :if && (ex = Expr(:call, ifelse, ex.args...))
    ex.head === :call || throw(ArgumentError("internal representation does not support non-call Expr"))

    op = eval(ex.args[1])  # HACK
    args = convert.(Expression, ex.args[2:end])

    return Operation(op, args)
end
Base.convert(::Type{Expression}, x::Expression) = x
Base.convert(::Type{Expression}, x::Number) = Constant(x)

function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end

function build_function(rhss, vs, ps = (), args = (), conv = simplified_expr, expression = Val{true}; constructor=nothing)
    _vs = map(x-> x isa Operation ? x.op : x, vs)
    _ps = map(x-> x isa Operation ? x.op : x, ps)
    var_pairs   = [(u.name, :(u[$i])) for (i, u) ∈ enumerate(_vs)]
    param_pairs = [(p.name, :(p[$i])) for (i, p) ∈ enumerate(_ps)]
    (ls, rs) = zip(var_pairs..., param_pairs...)

    var_eqs = Expr(:(=), build_expr(:tuple, ls), build_expr(:tuple, rs))

    fname = gensym(:ModelingToolkitFunction)

    X = gensym(:MTIIPVar)
    ip_sys_exprs = [:($X[$i] = $(conv(rhs))) for (i, rhs) ∈ enumerate(rhss)]
    ip_let_expr = Expr(:let, var_eqs, build_expr(:block, ip_sys_exprs))

    sys_expr = build_expr(:tuple, [conv(rhs) for rhs ∈ rhss])
    let_expr = Expr(:let, var_eqs, sys_expr)

    fargs = ps == () ? :(u,$(args...)) : :(u,p,$(args...))

    oop_ex = :(
        ($(fargs.args...),) -> begin
            X = $let_expr
            T = promote_type(map(typeof,X)...)
            convert.(T,X)
            construct = $(constructor === nothing ? :(u isa ModelingToolkit.StaticArrays.StaticArray ? ModelingToolkit.StaticArrays.similar_type(typeof(u), eltype(X)) : x->(du=similar(u, T, $(size(rhss)...)); vec(du) .= x; du)) : constructor)
            construct(X)
        end
    )

    iip_ex = :(
        ($X,$(fargs.args...)) -> begin
            $ip_let_expr
            nothing
        end
    )

    if expression == Val{true}
        return oop_ex, iip_ex
    else
        return GeneralizedGenerated.mk_function(@__MODULE__,oop_ex), GeneralizedGenerated.mk_function(@__MODULE__,iip_ex)
    end
end

is_constant(::Constant) = true
is_constant(::Any) = false

is_operation(::Operation) = true
is_operation(::Any) = false

is_derivative(O::Operation) = isa(O.op, Differential)
is_derivative(::Any) = false

Base.occursin(t::Expression) = Base.Fix1(occursin, t)
Base.occursin(t::Expression, x::Operation ) = isequal(x, t) || any(occursin(t), x.args)
Base.occursin(t::Expression, x::Expression) = isequal(x, t)

clean(x::Variable) = x
clean(O::Operation) = isa(O.op, Variable) ? O.op : throw(ArgumentError("invalid variable: $(O.op)"))

vars(exprs) = foldl(vars!, exprs; init = Set{Variable}())
function vars!(vars, O)
    isa(O, Operation) || return vars
    for arg ∈ O.args
        if isa(arg, Operation)
            isa(arg.op, Variable) && push!(vars, arg.op)
            vars!(vars, arg)
        end
    end

    return vars
end


text:
module DecFP
export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str

using Compat

const libbid = joinpath(dirname(@__FILE__), "..", "deps", "libbid$(Sys.WORD_SIZE)")

const _buffer = Array(UInt8, 1024)

import Base.promote_rule
import Core.Intrinsics: box, unbox, bswap_int

# global pointers and dicts must be initialized at runtime (via __init__)
function __init__()
    global const rounding = cglobal((:__bid_IDEC_glbround, libbid), Cuint) # rounding mode
    global const flags = cglobal((:__bid_IDEC_glbflags, libbid), Cuint) # exception status
    unsafe_store!(flags, 0)

    # rounding modes, from bid_functions.h
    global const rounding_c2j = [RoundNearest, RoundDown, RoundUp, RoundToZero, RoundFromZero]
    global const rounding_j2c = Dict{RoundingMode, UInt32}([(rounding_c2j[i], Cuint(i-1)) for i in 1:length(rounding_c2j)])
end

# status flags from bid_functions.h:
const INVALID    = 0x01
const UNNORMAL   = 0x02
const DIVBYZERO  = 0x04
const OVERFLOW   = 0x08
const UNDERFLOW  = 0x10
const INEXACT    = 0x20

bidsym(w,s...) = string("__bid", w, "_", s...)

abstract DecimalFloatingPoint <: AbstractFloat
Base.get_rounding{T<:DecimalFloatingPoint}(::Type{T}) = rounding_c2j[unsafe_load(rounding)+1]
Base.set_rounding{T<:DecimalFloatingPoint}(::Type{T}, r::RoundingMode) = unsafe_store!(rounding, rounding_j2c[r])

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    @eval bitstype $w $BID <: DecimalFloatingPoint
end

# quickly check whether s begins with "±nan"
function isnanstr(s::AbstractString)
    i = start(s)
    while !done(s, i)
        c, i = next(s, i)
        isspace(c) || break
    end
    done(s, i) && return false
    if (c == '+' || c == '-')
        c, i = next(s, i)
        done(s, i) && return false
    end
    lowercase(c) == 'n' || return false
    c, i = next(s, i)
    (!done(s, i) && lowercase(c) == 'a') || return false
    c, i = next(s, i)
    (done(s, i) && lowercase(c) == 'n') || return false
    return true
end

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    T = eval(BID)
    Ti = eval(Symbol(string("UInt",w)))

    # hack: we need an internal parsing function that doesn't check exceptions, since
    # flags isn't defined until __init__ runs.  Similarly for nextfloat/prevfloat
    @eval begin
        _parse(::Type{$BID}, s::AbstractString) =
            ccall(($(bidsym(w,"from_string")), libbid), $BID, (Ptr{UInt8},), s)
        _nextfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128), x, pinf128)
        _prevfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128), x, minf128)
        _sub(x::$BID, y::$BID) = ccall(($(bidsym(w,"sub")), libbid), $BID, ($BID,$BID), x, y)
    end

    @eval begin
        function Base.parse(::Type{$BID}, s::AbstractString)
            x = _parse($BID, s)
            if isnan(x) && !isnanstr(s)
                throw(ArgumentError("invalid number format $s"))
            end
            return xchk(x, InexactError)
        end

        function Base.show(io::IO, x::$BID)
            ccall(($(bidsym(w,"to_string")), libbid), Void, (Ptr{UInt8}, $BID), _buffer, x)
            unsafe_write(io, pointer(_buffer), ccall(:strlen, Csize_t, (Ptr{UInt8},), _buffer))
        end

        Base.fma(x::$BID, y::$BID, z::$BID) = nox(ccall(($(bidsym(w,"fma")), libbid), $BID, ($BID,$BID,$BID), x, y, z))
        Base.muladd(x::$BID, y::$BID, z::$BID) = fma(x,y,z) # faster than x+y*z

        Base.one(::Union{Type{$BID},$BID}) = $(_parse(T, "1"))
        Base.zero(::Union{Type{$BID},$BID}) = $(_parse(T, "0"))

        Base.signbit(x::$BID) = $(zero(Ti)) != $(Ti(1) << (Ti(w - 1))) & reinterpret($Ti, x)
        Base.sign(x::$BID) = ifelse(signbit(x), $(_parse(T, "-1")), $(_parse(T, "1")))

        Base.nextfloat(x::$BID) = nox(_nextfloat(x))
        Base.prevfloat(x::$BID) = nox(_prevfloat(x))
        Base.eps(x::$BID) = ifelse(isfinite(x), xchk(nextfloat(x) - x, OVERFLOW), $(_parse(T, "NaN")))

        # the meaning of the exponent is different than for binary FP: it is 10^n, not 2^n:
        # Base.exponent(x::$BID) = nox(ccall(($(bidsym(w,"ilogb")), libbid), Cint, ($BID,), x))
        # Base.ldexp(x::$BID, n::Integer) = nox(ccall(($(bidsym(w,"ldexp")), libbid), $BID, ($BID,Cint), x, n))
    end

    for (f,c) in ((:isnan,"isNaN"), (:isinf,"isInf"), (:isfinite,"isFinite"), (:issubnormal,"isSubnormal"))
        @eval Base.$f(x::$BID) = ccall(($(bidsym(w,c)), libbid), Cint, ($BID,), x) != 0
    end

    for (f,c) in ((:+,"add"), (:-,"sub"), (:*,"mul"), (:/, "div"), (:hypot,"hypot"), (:atan2,"atan2"), (:^,"pow"), (:copysign,"copySign"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,$BID), x, y))
    end

    for f in (:exp,:log,:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,:atanh,:log1p,:expm1,:log10,:log2,:exp2,:exp10,:erf,:erfc,:lgamma,:sqrt,:cbrt,:abs)
        @eval Base.$f(x::$BID) = xchk(ccall(($(bidsym(w,f)), libbid), $BID, ($BID,), x), INVALID)
    end
    for (f,c) in ((:gamma,"tgamma"), (:-,"negate"), (:round,"nearbyint"))
        @eval Base.$f(x::$BID) = xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,), x), INVALID)
    end

    for (f,c) in ((:(==),"quiet_equal"), (:>,"quiet_greater"), (:<,"quiet_less"), (:(>=), "quiet_greater_equal"), (:(<=), "quiet_less_equal"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), Cint, ($BID,$BID), x, y) != 0)
    end

    for Tf in (Float32,Float64)
        bT = string("binary",sizeof(Tf)*8)
        @eval begin
            Base.convert(::Type{$Tf}, x::$BID) = nox(ccall(($(bidsym(w,"to_",bT)), libbid), $Tf, ($BID,), x))
            Base.convert(::Type{$BID}, x::$Tf) = nox(ccall(($(string("__",bT,"_to_","bid",w)), libbid), $BID, ($Tf,), x))
        end
    end

    for c in (:π, :e, :γ, :catalan, :φ)
        @eval begin
            Base.convert(::Type{$BID}, ::Irrational{$(QuoteNode(c))}) = $(_parse(T, @compat setprecision(256) do
                                                                                string(BigFloat(eval(c)))
                                                                                end))
        end
    end

    @eval promote_rule(::Type{$BID}, ::Type{Irrational}) = $BID

    for w′ in (32,64,128)
        BID′ = Symbol(string("Dec",w′))
        if w > w′
            @eval promote_rule(::Type{$BID}, ::Type{$BID′}) = $BID
        end
        if w != w′
            @eval Base.convert(::Type{$BID}, x::$BID′) = xchk(ccall(($(string("__bid",w′,"_to_","bid",w)), libbid), $BID, ($BID′,), x), INEXACT)
        end

        # promote binary*decimal -> decimal, for consistency with other operations above
        # (there doesn't seem to be any clear standard for this)
        if w′ <= 64
            FP′ = Symbol(string("Float",w′))
            @eval promote_rule(::Type{$BID}, ::Type{$FP′}) = $(Symbol(string("Dec",max(w,w′))))
            for i′ in ("Int$w′", "UInt$w′")
                Ti′ = eval(Symbol(i′))
                @eval begin
                    Base.convert(::Type{$BID}, x::$Ti′) = nox(ccall(($(bidsym(w,"from_",lowercase(i′))), libbid), $BID, ($Ti′,), x))
                end
            end
        end
    end

    for w′ in (8,16,32,64)
        for i′ in ("Int$w′", "UInt$w′")
            Ti′ = eval(Symbol(i′))
            @eval begin
                Base.floor(::Type{$Ti′}, x::$BID) = xchk(ccall(($(bidsym(w,"to_",lowercase(i′),"_xfloor")), libbid), $Ti′, ($BID,), x), InexactError, INVALID | OVERFLOW)
                Base.ceil(::Type{$Ti′}, x::$BID) = xchk(ccall(($(bidsym(w,"to_",lowercase(i′),"_xceil")), libbid), $Ti′, ($BID,), x), InexactError, INVALID | OVERFLOW)
                Base.convert(::Type{$Ti′}, x::$BID) = xchk(ccall(($(bidsym(w,"to_",lowercase(i′),"_xfloor")), libbid), $Ti′, ($BID,), x), InexactError)
            end
        end
    end

    @eval Base.bswap(x::$BID) = reinterpret($BID, bswap(reinterpret($Ti, x)))
    @eval Base.convert(::Type{Float16}, x::$BID) = convert(Float16, convert(Float32, x))
end # widths w

# the complex-sqrt function in base doesn't work for use, because it requires base-2 ldexp
function Base.sqrt{T<:DecimalFloatingPoint}(z::Complex{T})
    x, y = reim(z)
    x==y==0 && return Complex(zero(x),y)
    ρ = sqrt((abs(x) + hypot(x,y)) * 0.5)
    ξ = ρ
    η = y
    if isfinite(η) η=(η/ρ)/2 end
    if x<0
        ξ = abs(η)
        η = copysign(ρ,y)
    end
    Complex(ξ,η)
end

# used for next/prevfloat:
const pinf128 = _parse(Dec128, "+Inf")
const minf128 = _parse(Dec128, "-Inf")

for T in (Dec32,Dec64,Dec128)
    @eval begin
        Base.eps(::Type{$T}) = $(_sub(_nextfloat(one(T)), one(T)))
        Base.typemax(::Type{$T}) = $(_parse(T, "+inf"))
        Base.typemin(::Type{$T}) = $(_parse(T, "-inf"))
        Base.realmax(::Type{$T}) = $(_prevfloat(_parse(T, "+inf")))
        Base.realmin(::Type{$T}) = $(_nextfloat(zero(T)))
    end
end

Base.convert{F<:DecimalFloatingPoint}(T::Type{F}, x::Union{Int8,UInt8,Int16,UInt16}) = F(Int32(x))
Base.convert{F<:DecimalFloatingPoint}(T::Type{F}, x::Float16) = F(Float32(x))
promote_rule{F<:DecimalFloatingPoint}(::Type{F}, ::Type{Float16}) = F
promote_rule{F<:DecimalFloatingPoint,T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64}}(::Type{F}, ::Type{T}) = F

# so that mathconsts get promoted to Dec32, not Dec64, like Float32
promote_rule{s,F<:DecimalFloatingPoint}(::Type{Irrational{s}}, ::Type{F}) = F
promote_rule{s,F<:DecimalFloatingPoint}(::Type{Irrational{s}}, T::Type{Complex{F}}) = T

macro d_str(s, flags...) parse(Dec64, s) end
macro d32_str(s, flags...) parse(Dec32, s) end
macro d64_str(s, flags...) parse(Dec64, s) end
macro d128_str(s, flags...) parse(Dec128, s) end

# clear exception flags and return x
function nox(x)
    unsafe_store!(flags, 0)
    return x
end

# check exception flags in mask & throw, otherwise returning x;
# always clearing exceptions
function xchk(x, mask::Integer=0x3f)
    f = unsafe_load(flags)
    unsafe_store!(flags, 0)
    if f & mask != 0
        f & INEXACT != 0 && throw(InexactError())
        f & OVERFLOW != 0 && throw(OverflowError())
        f & DIVBYZERO != 0 && throw(DivideError())
        f & INVALID != 0 && throw(DomainError())
        f & UNDERFLOW != 0 && error("underflow")
        f & UNNORMAL != 0 && error("unnormal")
    end
    return x
end

function xchk{E<:Exception}(x, exc::Type{E}, mask::Integer=0x3f)
    f = unsafe_load(flags)
    unsafe_store!(flags, 0)
    f & mask != 0 && throw(exc())
    return x
end

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

mutable struct Table
    rows::Vector{Vector{Any}}
    align::Vector{Symbol}
end

function parserow(stream::IO)
    withstream(stream) do
        line = readline(stream)
        row = split(line, r"(?<!\\)\|")
        length(row) == 1 && return
        isempty(row[1]) && popfirst!(row)
        map!(x -> strip(replace(x, "\\|" => "|")), row, row)
        isempty(row[end]) && pop!(row)
        return row
    end
end

function rowlength!(row, len)
    while length(row) < len push!(row, "") end
    while length(row) > len pop!(row) end
    return row
end

const default_align = :r

function parsealign(row)
    align = Symbol[]
    for s in row
        (length(s) ≥ 3 && s ⊆ Set("-:")) || return
        push!(align,
              s[1] == ':' ? (s[end] == ':' ? :c : :l) :
              s[end] == ':' ? :r :
              default_align)
    end
    return align
end

function github_table(stream::IO, md::MD)
    withstream(stream) do
        skipblank(stream)
        rows = []
        cols = 0
        align = nothing
        while (row = parserow(stream)) !== nothing
            if length(rows) == 0
                isempty(row[1]) && return false
                cols = length(row)
            end
            if align === nothing && length(rows) == 1 # Must have a --- row
                align = parsealign(row)
                (align === nothing || length(align) != cols) && return false
            else
                push!(rows, map(x -> parseinline(x, md), rowlength!(row, cols)))
            end
        end
        length(rows) <= 1 && return false
        push!(md, Table(rows, align))
        return true
    end
end

function html(io::IO, md::Table)
    withtag(io, :table) do
        for (i, row) in enumerate(md.rows)
            withtag(io, :tr) do
                for c in md.rows[i]
                    withtag(io, i == 1 ? :th : :td) do
                        htmlinline(io, c)
                    end
                end
            end
        end
    end
end

mapmap(f, xss) = map(xs->map(f, xs), xss)

colwidths(rows; len = length, min = 0) =
    reduce((x,y) -> max.(x,y), [min; convert(Vector{Vector{Int}}, mapmap(len, rows))])

padding(width, twidth, a) =
    a == :l ? (0, twidth - width) :
    a == :r ? (twidth - width, 0) :
    a == :c ? (floor(Int, (twidth-width)/2), ceil(Int, (twidth-width)/2)) :
    error("Invalid alignment $a")

function padcells!(rows, align; len = length, min = 0)
    widths = colwidths(rows, len = len, min = min)
    for i = 1:length(rows), j = axes(rows[1],1)
        cell = rows[i][j]
        lpad, rpad = padding(len(cell), widths[j], align[j])
        rows[i][j] = " "^lpad * cell * " "^rpad
    end
    return rows
end

_dash(width, align) =
    align == :l ? ":" * "-"^width * " " :
    align == :r ? " " * "-"^width * ":" :
    align == :c ? ":" * "-"^width * ":" :
    throw(ArgumentError("Invalid alignment $align"))

function plain(io::IO, md::Table)
    cells = mapmap(md.rows) do each
        replace(plaininline(each), "|" => "\\|")
    end
    padcells!(cells, md.align, len = length, min = 3)
    for i = axes(cells,1)
        print(io, "| ")
        join(io, cells[i], " | ")
        println(io, " |")
        if i == 1
            print(io, "|")
            join(io, [_dash(length(cells[i][j]), md.align[j]) for j = axes(cells[1],1)], "|")
            println(io, "|")
        end
    end
end

function rst(io::IO, md::Table)
    cells = mapmap(rstinline, md.rows)
    padcells!(cells, md.align, len = length, min = 3)
    single = ["-"^length(c) for c in cells[1]]
    double = ["="^length(c) for c in cells[1]]
    function print_row(row, row_sep, col_sep)
        print(io, col_sep, row_sep)
        join(io, row, string(row_sep, col_sep, row_sep))
        println(io, row_sep, col_sep)
    end
    print_row(single, '-', '+')
    for i = 1:length(cells)
        print_row(cells[i], ' ', '|')
        i ≡ 1 ? print_row(double, '=', '+') :
                print_row(single, '-', '+')
    end
end

function term(io::IO, md::Table, columns)
    cells = mapmap(x -> terminline_string(io, x), md.rows)
    padcells!(cells, md.align, len = ansi_length)
    for i = 1:length(cells)
        join(io, cells[i], " ")
        if i == 1
            println(io)
            join(io, ["–"^ansi_length(cells[i][j]) for j = 1:length(cells[1])], " ")
        end
        i < length(cells) && println(io)
    end
end

function latex(io::IO, md::Table)
    wrapblock(io, "tabular") do
        align = md.align
        println(io, "{$(join(align, " | "))}")
        for (i, row) in enumerate(md.rows)
            for (j, cell) in enumerate(row)
                j != 1 && print(io, " & ")
                latexinline(io, cell)
            end
            println(io, " \\\\")
            if i == 1
                println(io, "\\hline")
            end
        end
    end
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia