text:
module NK_baseline

    using NLsolve

    function steadystate(model) #,ys,exo)

        new_calibration = copy(model.calibration)

        for p in keys(model.calibration)
            v = model.calibration[p]
            expr = :($p=$v)
            eval(expr)
        end
        eta = model.calibration[:eta] # workaround

        PI=PIbar
        u=1
        q=1
        d=1
        phi=1
        m=0
        zeta=1
        mu_z=exp(LambdaYd)
        mu_I=exp(Lambdamu)
        mu_A=exp(LambdaA)

        # %set the parameter Lambdax
        Lambdax=mu_z

        # %set the parameter gammma1
        gammma1=mu_z*mu_I/betta-(1-delta)

        r=1*gammma1
        R=1+(PI*mu_z/betta-1)

        Rbar=R

        PIstar=((1-thetap*PI^(-(1-epsilon)*(1-chi)))/(1-thetap))^(1/(1-epsilon))
        PIstarw=((1-thetaw*PI^(-(1-chiw)*(1-eta))*mu_z^(-(1-eta)))/(1-thetaw))^(1/(1-eta))

        mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^((1-chi)*epsilon))/(1-betta*thetap*PI^(-(1-epsilon)*(1-chi)))*PIstar
        w=(1-alppha)*(mc*(alppha/r)^alppha)^(1/(1-alppha))
        wstar=w*PIstarw
        vp=(1-thetap)/(1-thetap*PI^((1-chi)*epsilon))*PIstar^(-epsilon)
        vw=(1-thetaw)/(1-thetaw*PI^((1-chiw)*eta)*mu_z^eta)*PIstarw^(-eta)
        tempvaromega=alppha/(1-alppha)*w/r*mu_z*mu_I

        fobj(ld::Float64) = ( (1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-betta*thetaw*mu_z^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))
        -(eta-1)/eta*wstar/(varpsi*PIstarw^(-eta*gammma)*ld^gammma)*((1-h*mu_z^(-1))^(-1)-betta*h*(mu_z-h)^(-1))*
        ((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1) )
        fobj1(ld::Vector) = [fobj(ld[1])]

        sol = nlsolve(not_in_place(fobj1), [0.25])
        ld = sol.zero[1]

        l=vw*ld
        k=tempvaromega*ld
        x=(1-(1-delta)*(mu_z*mu_I)^(-1))*k
        yd=(mu_A/mu_z*k^alppha*ld^(1-alppha)-Phi)/vp
        c=(mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi
        lambda=(1-h*betta*mu_z^(-1))*(1-h/mu_z)^(-1)/c
        F=yd-1/(1-alppha)*w*ld
        f=(eta-1)/eta*wstar*PIstarw^(-eta)*lambda*ld/(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))
        f2=varpsi*d*phi*PIstarw^(-eta*(1+gammma))*ld^(1+gammma)/(1-betta*thetaw*(PI^chiw/PI)^(-eta*(1+gammma))*(wstar/wstar*mu_z)^(eta*(1+gammma)))

        g1=lambda*mc*yd/(1-betta*thetap*PI^((1-chi)*epsilon))
        g2=epsilon/(epsilon-1)*g1

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

        new_calibration[:PI] = PI
        new_calibration[:u] = u
        new_calibration[:q] = q
        new_calibration[:d] = d
        new_calibration[:phi] = phi
        new_calibration[:m] = m
        new_calibration[:zeta] = zeta
        new_calibration[:mu_z] = mu_z
        new_calibration[:mu_I] = mu_I
        new_calibration[:mu_A] = mu_A
        new_calibration[:Lambdax] = Lambdax
        new_calibration[:gammma1] = gammma1
        new_calibration[:r] = r
        new_calibration[:R] = R
        new_calibration[:Rbar] = Rbar
        new_calibration[:PIstar] = PIstar
        new_calibration[:PIstarw] = PIstarw
        new_calibration[:mc] = mc
        new_calibration[:w] = w
        new_calibration[:wstar] = wstar
        new_calibration[:vp] = vp
        new_calibration[:vw] = vw
        new_calibration[:tempvaromega] = tempvaromega
        new_calibration[:ld] = ld
        new_calibration[:l] = l
        new_calibration[:k] = k
        new_calibration[:x] = x
        new_calibration[:yd] = yd
        new_calibration[:c] = c
        new_calibration[:lambda] = lambda
        new_calibration[:F] = F
        new_calibration[:f] = f
        new_calibration[:f2] = f2
        new_calibration[:g1] = g1
        new_calibration[:g2] = g2

        return new_calibration

    end

    # not very satisfying and not to be imitated
    function update_calibration(model, new_calibration)
        for k in keys(new_calibration)
            model.calibration[k] = new_calibration[k]
        end
        for i=1:length(model.symbols[:endogenous])
            s = model.symbols[:endogenous][i]
            v = model.calibration[s]
            model.calibration_grouped[:endogenous][i] = v
        end
        for i=1:length(model.symbols[:parameters])
            s = model.symbols[:parameters][i]
            v = model.calibration[s]
            model.calibration_grouped[:parameters][i] = v
        end
    end

end


text:
function __init__()
    @require CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" include("Initialization/init_CDDLib.jl")
    @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" include("Initialization/init_Distributions.jl")
    @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("Initialization/init_Expokit.jl")
    @require Javis = "78b212ba-a7f9-42d4-b726-60726080707e" include("Initialization/init_Javis.jl")
    @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("Initialization/init_Makie.jl")
    @require MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" include("Initialization/init_MiniQhull.jl")
    @require Optim = "429524aa-4258-5aef-a3af-852621145aeb" include("Initialization/init_Optim.jl")
    @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" include("Initialization/init_Polyhedra.jl")
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("Initialization/init_StaticArrays.jl")
    @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("Initialization/init_Symbolics.jl")
    @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("Initialization/init_WriteVTK.jl")
end

ext: jl


text:
@safetestset "Sessions id test" begin
    using Genie

    @test !isempty(Genie.secret_token())
    @test Genie.Sessions.id() != Genie.Sessions.id()
end;


text:
"""
    rankinequalities(m::Model; max_ranks::Int = 5, suppress_warnings::Bool = true, U::Int = 10000, diag::Bool = false, active_constr::Bool = false, real_recourse_dim::Int = -1, conic_δ::Bool = false)
Returns ranking data in the form `Vector{Dict}` where each dictionary corresponds to a particular rank level and
contains the optimal flexibility_index, active constraint indexes, and flexibility model. The function will iteratively
solve copies of the flexibility model via `solve` where the prior active constraints are turned off in order rank the most limiting
constraints. The user can specify the maximum number of rank levels and the flexibility index problem will be repeatedly
solved until that maximum is acheived or the problem becomes unbounded, which occurs first.

**Arguments**
- `m::Model` The flexibility model.

**Keyword Arguments**
- `max_ranks::Int = 5` The maximum number of rank levels. 2
- `suppress_warnings::Bool = false` Indicates if solver warnings should be suppressed.
- `U::Union{Int, Float64} = 10000` The slack variable upper bound.
- `diag::Bool = false` Indicates whether or not to diagnonalize ellipsoidal uncertainty set (this is only active when an ellipsoidal set is used).
- `active_constr::Bool = false` Indicates if the optional active constraint should be used which enforces how many inequalities are active at the solution, this must be set to `true` for systems without control variables and/or contain state variables.
- `real_recourse_dim::Int = -1` The actual number of recourse variables in case state variables are included as recourse variables. This is mandatory if `active_constr = true` and no state variables are provided.
- `conic_δ::Bool = false` This should be set to `true` if a conic solver is used such as Pajarito.jl.

```julia
julia> rankinequalities(m, active_constr = true)
2-element Array{Dict,1}:
 Dict{String,Any}(Pair{String,Any}("flexibility_index", 3.59938),Pair{String,Any}("model", Feasibility problem with:
 * 6 linear constraints
 * 6 variables
Solver is Pavito),Pair{String,Any}("active_constraints", [3, 6]))
 Dict{String,Any}(Pair{String,Any}("flexibility_index", 9.58983),Pair{String,Any}("model", Feasibility problem with:
 * 6 linear constraints
 * 6 variables
Solver is Pavito),Pair{String,Any}("active_constraints", [1, 5]))
```
"""
function rankinequalities(m::Model; max_ranks::Int = 5, suppress_warnings::Bool = true, U::Int = 10000, diag::Bool = false, active_constr::Bool = false, real_recourse_dim::Int = -1, conic_δ::Bool = false)
    # Check the max_ranks
    if max_ranks < 1
        error("max_ranks set to $max_ranks, but must be an integer greater than 0.")
    end

    # Initialize data
    m_copy = deepcopy(m)
    ranked_data = Vector{Dict}(undef, max_ranks)
    ranked_constrs = []
    inactives = []
    counter = 0

    # Get the inequality information
    constr_bounds = JuMP.prepConstrBounds(m_copy)
    inequal_inds = findall(constr_bounds[1] .!= constr_bounds[2])

    # Iteratively solve the flexibility index problem and extract data
    for i = 1:max_ranks
        status = solve(m_copy, suppress_warnings = suppress_warnings, U = U, diag = diag, active_constr = active_constr, real_recourse_dim = real_recourse_dim, conic_δ = conic_δ, inactives = inactives)
        if status == :Optimal
            ranked_constrs = [findall(inequal_inds .== getactiveconstraints(m_copy)[j])[1] for j = 1:length(getactiveconstraints(m_copy))]
            inactives = unique([inactives; ranked_constrs])
            ranked_data[i] = Dict("flexibility_index" => getflexibilityindex(m_copy), "active_constraints" => getactiveconstraints(m_copy), "model" => m_copy)
            counter += 1
        else
            break
        end
        if length(inactives) == length(inequal_inds)
            break
        end
    end

    # Return the results
    return ranked_data[1:counter]
end

"""
    getconfidencelevel(m::Model)
Return the confidence level provided that the flexibility model `m` has been solved with an
ellipsoidal uncertainty set. This is equivalent to the quantile associated with the optimized
uncertainty set. Note that this assumes a multivariate Gaussian distribution with mean
`FlexibilityData.RVmeans` and covariance `FlexibilityData.covariance`.

**Arguments**
- `m::Model` The flexibility model.

```julia
julia> getconfidencelevel(m)
0.5370703369769008
```
"""
function getconfidencelevel(m::Model)
    flex_data = getflexibilitydata(m)
    if flex_data.uncertainty_set.name == :Ellipsoid && flex_data.flexibility_index != nothing
        return cdf(Chisq(flex_data.numRVs), flex_data.flexibility_index)
    elseif flex_data.flexibility_index == nothing
        error("Cannot compute the confidence level without first solving the flexibility model.")
    else
        error("Invalid uncertainty set, must solve the flexibility model with an ellipsoidal set.")
    end
end

"""
    findstochasticflexibility(m::Model; num_pts::Int = 10000, toler::Number = 1e-5, solver = Clp.ClpSolver(), only_positive::Bool = false, use_vulnerability_model::Bool = false, use_flexibility_index::Bool = false, seed::Int = -1)
Returns the estimated stochastic flexibility index that is evaluated via Monte Carlo sampling. At default
this estimation is carried out by evaluating the feasibility of each Monte Carlo sample. The samples are
generating from a multivariate Gaussian distribution with mean `FlexibilityData.RVmeans` and covariance
`FlexibilityData.covariance`. The vulnerability model also tests the feasibility of the samples, but does
so in one large optimization problem instead of evaluating each sample indiviually. The optimized flexibility
index can also be used to reduce the number of samples that need to be evaluated.

**Arguments**
- `m::Model` The flexibility model.

**Keyword Arguments**
- `num_pts::Int = 10000` Number of Monte Carlo samples.
- `toler::Number = 1e-5` The feasibility check tolerance.
- `solver = Clp.ClpSolver()` The solver, any LP or NLP solver should work.
- `only_positive::Bool = false` Indicates if only positive samples should be used.
- `use_vulnerability_model::Bool = false` Indicates if the vulnerability model should be used.
- `use_flexibility_index::Bool = false` Indicates if the optimal flexibility index should be used.
- `seed::Int = -1` Random seed for sample collection, any negative value will turn off the random seed.

```julia
julia> findstochasticflexibility(m)
0.9687

julia> findstochasticflexibility(m, use_vulnerability_model = true)
0.9705

julia> findstochasticflexibility(m, num_pts = 5000, use_flexibility_index = true)
0.973
```
"""
function findstochasticflexibility(m::Model; num_pts::Int = 10000, toler::Number = 1e-5, solver = Clp.ClpSolver(), only_positive::Bool = false, use_vulnerability_model::Bool = false, use_flexibility_index::Bool = false, seed::Int = -1)
    # Parse data
    flex_data = getflexibilitydata(m)

    if length(flex_data.covariance) == 0
        error("The covariance matrix needs to be set.")
    end

    # Load in the neceassary data
    input_dict = MakeInputDict(m)
    θ_nom = input_dict["theta_nom"]

    # Parse the necessary data dimensions
    n_z = input_dict["dims"][1]
    n_θ = input_dict["dims"][2]
    n_x = input_dict["dims"][3]
    n_f = input_dict["dims"][4]
    n_h = input_dict["dims"][5]

    # Setup a multivariate normal distribution and Monte Carlo samples
    if seed >= 0
        Random.seed!(seed)
    end
    d = MvNormal(Vector{Float64}(flex_data.RVmeans), Matrix{Float64}(getcovariance(m)))
    if only_positive
        samples = Array{Float64}(undef, flex_data.numRVs, num_pts)
        not_done = true
        counter = 1
        while not_done
            sample = rand(d, 1)
            if all(sample .>= 0)
                samples[:, counter] = sample
                counter += 1
            end
            if counter > num_pts
                not_done = false
            end
        end
    else
        samples = rand(d, num_pts)
    end

    # Determine which points are outside of the set and check for F
    F = getflexibilityindex(m)
    if use_flexibility_index && F == nothing
        @warn "Flexibility index hasn't yet been computed. Setting use_flexibility_index = false."
        use_flexibility_index = false
    elseif use_flexibility_index && flex_data.uncertainty_set.name == :Ellipsoid
        inv_covar = inv(getcovariance(m))
        outside_set = [(samples[:, k] - θ_nom)' * inv_covar * (samples[:, k] - θ_nom) > F for k = 1:num_pts]
    elseif use_flexibility_index && flex_data.uncertainty_set.name == :PNorm
        outside_set = [norm(samples[:, k] - θ_nom, flex_data.uncertainty_set.p) > F for k = 1:num_pts]
    elseif use_flexibility_index && flex_data.uncertainty_set.name == :Hyperbox
        outside_set = Vector{Bool}(undef, num_pts)
        for k = 1:num_pts
            all_inside = true
            for i = 1:flex_data.numRVs
                if samples[i, k] < θ_nom[i] - flex_data.uncertainty_set.neg_dev[i] * F || samples[i, k] > θ_nom[i] + flex_data.uncertainty_set.pos_dev[i] * F
                    all_inside = false
                end
            end
            outside_set[k] = !all_inside
        end
    end

    # Determine the points to test
    if use_flexibility_index
        test_pts = samples[:, outside_set]
        num_test_pts = size(test_pts)[2]
    else
        num_test_pts = num_pts
        test_pts = samples
    end

    # Stop if all the test points are feasible.
    if num_test_pts == 0
        return 1.0
    end

    # Setup the model to be solved
    m_solve = Model(solver = solver)

    if use_vulnerability_model
        # Initialize variables
        @variable(m_solve, θ[1:n_θ, 1:num_test_pts])
        @variable(m_solve, u[1:num_test_pts] >= 0)
        if n_z != 0
            @variable(m_solve, z[1:n_z, 1:num_test_pts])
        end
        if n_x != 0
            @variable(m_solve, x[1:n_x, 1:num_test_pts])
        end

        # Set the objective
        @objective(m_solve, Min, 1 / num_test_pts * sum(u[k] for k = 1:num_test_pts))

        # Add the constraints
        exprs = AddSystemExpressions(m_solve, input_dict, num_test_pts)
        @constraint(m_solve, f[j = 1:n_f, k = 1:num_test_pts], exprs[1][j, k] <= u[k])
        if n_h != 0
            @constraint(m_solve, h[i = 1:n_h, k = 1:num_test_pts], exprs[2][i, k] == 0)
        end
    else
        # Initialize variables
        @variable(m_solve, θ[1:n_θ])
        @variable(m_solve, u)
        if n_z != 0
            @variable(m_solve, z[1:n_z])
        end
        if n_x != 0
            @variable(m_solve, x[1:n_x])
        end

        # Set the objective
        @objective(m_solve, Min, u)

        # Add the constraints
        exprs = AddSystemExpressions(m_solve, input_dict)
        @constraint(m_solve, f[j = 1:n_f], exprs[1][j] <= u)
        if n_h != 0
            @constraint(m_solve, h[i = 1:n_h], exprs[2][i] == 0)
        end
    end

    # Solve the model(s) and compute the SF
    if use_vulnerability_model
        for i = 1:n_θ
            for k = 1:num_test_pts
                JuMP.fix(θ[i, k], test_pts[i, k])
            end
        end
        status = solve(m_solve)
        if status == :Optimal
            return 1 - sum(getvalue(u) .>= toler) / num_pts
        else
            error("Problem not solved to optimality, cannot estimate SF.")
        end
    else
        infeasible_counter = 0
        feasible_results = zeros(Bool, num_test_pts)
        for k = 1:num_test_pts
            for i = 1:n_θ
                JuMP.fix(θ[i], test_pts[i, k])
            end
            status = solve(m_solve)
            if status == :Optimal
                if getvalue(u) <= toler
                    feasible_results[k] = true
                end
            else
                infeasible_counter += 1
            end
        end
        if infeasible_counter != 0
            @warn "Not all scenario subproblems not solved to optimality, estmation of SF might not be correct."
        end
        if use_flexibility_index
            return (sum(feasible_results) + num_pts - num_test_pts) / (num_pts - infeasible_counter)
        else
            return sum(feasible_results) / (num_pts - infeasible_counter)
        end
    end
end


text:
"""
    print_blocks(filename::String; verbose::Bool=true)

Reads the block names of blocks in a snapshot and returns them in an array.
Outputs them to console if `verbose=true`
"""
function print_blocks(filename::String; verbose::Bool=true)

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize != 8
        error("Block search not possible - use snap_format 2!")
    end

    p = position(f)
    seek(f,4)

    blocks = []

    while eof(f) != true

        name = Char.(read!(f, Array{Int8,1}(undef,4)))
        blockname = String(name)

        blockname = strip(blockname)

        push!(blocks, blockname)

        p = position(f)
        seek(f,p+8)

        skipsize = read(f, Int32)

        p = position(f)
        seek(f,p+skipsize+8)

    end

    if verbose
        println("Found blocks: ")
        for block ∈ blocks
            println(block)
        end
    end

    close(f)

    return String.(blocks)
end

"""
    read_info(filename; verbose::Bool=false)

Reads the info block of a snapshot and returns the information in an array of `Info_Line` types.
If `verbose=true` the blocknames are also printed to console.
"""
function read_info(filename::String; verbose::Bool=false)

    f = open(filename)
    seek(f,4)

    while eof(f) != true

        name = Char.(read!(f, Array{Int8,1}(undef,4)))
        blockname = String(name)

        p = position(f)
        seek(f,p+8)

        skipsize = read(f, Int32)

        if blockname == "INFO"
            n_blocks = Int(skipsize/40) # one info line is 40 bytes
            arr_info = Array{Info_Line,1}(undef,n_blocks)

            for i = 1:n_blocks
                arr_info[i] = read_info_line(f)
            end # for

            close(f)

            if verbose == true
                println("Found Info block.\nEntries are:\n")
                for i ∈ 1:length(arr_info)
                    println(i, " - ", arr_info[i].block_name)
                end # for
            end # verbose

            return arr_info

        else
            p = position(f)
            seek(f,p+skipsize+8)
        end # if blockname == "INFO"

    end # while eof(f) != true

    println("No info block present!")

    return 1
end

function read_info_line(f)

    # block name consists of 4 C Chars.
    name = Char.(read!(f, Array{Int8,1}(undef,4)))
    blockname = String(name)

    block_name = strip(blockname)

    # the datatype is stored as a 8 char word.
    letters = read!(f, Array{Int8,1}(undef,8))
    letters = Char.(letters)
    letters = string.(letters)

    data_type = ""
    for i = 1:length(letters)
        data_type *= letters[i]
    end

    # erase whitespace and make lower case for comparison
    data_type = lowercase(strip(data_type))

    # assign data types
    if data_type == "float" ||
       data_type == "floatn"
            dt = Float32
    elseif data_type == "long"
            dt = UInt32
    elseif data_type == "llong"
            dt = UInt64
    elseif data_type == "double" ||
           data_type == "doublen"
            dt = Float64
    end

    # array dimensions are stored as Int32
    n_dim = read(f,Int32)

    # read the information for which particles the block is relevant.
    # e.g. for gas particles: [ 1, 0, 0, 0, 0, 0 ]
    is_present = read!(f, Array{Int32,1}(undef,6))

    # construct the info line struct and return it.
    return Info_Line(block_name, dt, n_dim, is_present)
end

"""
    block_present(filename::String, blockname::String, blocks::Vector{String}=[""])

Checks if a given block is present in the snapshot file, or in the supplied `blocks` Vector.
"""
function block_present(filename::String, blockname::String, blocks::Vector{String}=[""])

    # if no blocks supplied, read in block information
    if blocks == [""]
        blocks = print_blocks(filename, verbose=false)
    end

    # loop over all blocks
    for block ∈ blocks
        if block == blockname
            return true
        end
    end

    return false
end

"""
    get_block_positions(filename::String)

Returns a dictionary with the starting positions of all blocks in a snapshot in bits.
"""
function get_block_positions(filename::String)

    f = open(filename)
    blocksize = read(f, Int32)

    # only works for snap format 2
    if blocksize != 8
        error("Block search not possible - use snap_format 2!")
    end

    p = position(f)
    seek(f,4)

    blocks = Vector{String}(undef, 0)
    pos    = Vector{Int64}(undef, 0)

    while eof(f) != true

        # read block name
        name = Char.(read!(f, Array{Int8,1}(undef,4)))
        blockname = String(name)

        blockname = strip(blockname)

        # store block name in array
        push!(blocks, blockname)

        read(f, Int32)
        read(f, Int32)

        skipsize = read(f, Int32)

        p = position(f)

        push!(pos, p)

        seek(f,p+skipsize+8)

    end

    close(f)

    d = Dict( blocks[i] => pos[i] for i = 1:length(blocks))

    return d
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
function get_api_key()
    "USDA_QUICK_SURVEY_KEY" in keys(ENV) ? (println("Your API key is: " * ENV["USDA_QUICK_SURVEY_KEY"])) : (println("You don't have a key in your environment. Use set_api_key"))
end

function set_api_key(api_key::String)
    length(api_key) != 36 ? (throw("Key must be 36 characters long")) : nothing

    if "USDA_QUICK_SURVEY_KEY" ∉ keys(ENV) || length(ENV["USDA_QUICK_SURVEY_KEY"]) != 36
        ENV["USDA_QUICK_SURVEY_KEY"] = api_key
        println("local API key is set.")
    else
        throw("You already have a registered key")
    end
end


text:
using ArgParse

include("./src/raft.jl")
using Raft

## I'm just a follower to start with
ps = Raft.PersistentState(0, 0, [])
vs = Raft.VolatileState(0, 0)

s = ArgParseSettings("Runs a Raft server")

@add_arg_table s begin
  "--server_id", "-s"
      help = "id of this server"
      arg_type = UInt64
      required = true
  "--port", "-p"
      help = "port to listen on"
      required = true
      arg_type = UInt64
  "--buddies", "-b"
      nargs = '+'
      help = "ports of buddies to talk to"
      arg_type = UInt64
      required = true
  "--amleader", "-l"
      action = :store_true
      help = "This server should be the leader"
end

## Read config:
parsed_args = parse_args(ARGS, s)
for key in parsed_args
  @show key
end

buddies = Dict{Raft.Server, Raft.Address}()
i = UInt64(1)
for b in parsed_args["buddies"]
  buddies[i] = Raft.Address(ip"127.0.0.1", b)
  i+=1
end

my_config = Raft.ServerConfig(parsed_args["server_id"],
  parsed_args["port"], buddies)
am_leader = parsed_args["amleader"]

## Wait for leader to contact me
if !am_leader
  Raft.run_follower(my_config)
else
  Raft.run_leader(my_config)
end


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug, @sync

using GPUArrays

using LinearAlgebra

using BFloat16s

using CEnum

using Memoize

using DataStructures

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")
include("libcublas_deprecated.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

# cache for created, but unused handles
const handle_cache_lock = ReentrantLock()
const idle_handles = DefaultDict{CuContext,Vector{cublasHandle_t}}(()->cublasHandle_t[])
const idle_xt_handles = DefaultDict{Any,Vector{cublasHandle_t}}(()->cublasHandle_t[])

function math_mode!(handle, mode)
    flags = 0

    # https://github.com/facebookresearch/faiss/issues/1385
    if version(handle) > v"11"
        flags = CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION
    end

    flags |= if mode == CUDA.PEDANTIC_MATH
        # prevent use of tensor cores
        if version(handle) < v"11"
            CUBLAS_DEFAULT_MATH
        else
            CUBLAS_PEDANTIC_MATH
        end
    elseif mode == CUDA.DEFAULT_MATH
        # use tensor cores, but don't reduce precision
        if version(handle) < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_DEFAULT_MATH
        end
    elseif mode == CUDA.FAST_MATH
        # we'll additionally select a compute-mode with reduced precision whenever possible
        if version(handle) < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_TF32_TENSOR_OP_MATH
        end
    end

    cublasSetMathMode(handle, flags)

    return
end

function handle()
    state = CUDA.active_state()
    handle, stream, math_mode = get!(task_local_storage(), (:CUBLAS, state.context)) do
        new_handle = @lock handle_cache_lock begin
            if isempty(idle_handles[state.context])
                cublasCreate()
                # FIXME: use cublasSetWorkspace? cublasSetStream reset it.
            else
                pop!(idle_handles[state.context])
            end
        end

        finalizer(current_task()) do task
            @spinlock handle_cache_lock begin
                push!(idle_handles[state.context], new_handle)
            end
        end
        # TODO: cublasDestroy to preserve memory, or at exit?

        cublasSetStream_v2(new_handle, state.stream)

        math_mode!(new_handle, state.math_mode)

        new_handle, state.stream, state.math_mode
    end::Tuple{cublasHandle_t,CuStream,CUDA.MathMode}

    if stream != state.stream
        cublasSetStream_v2(handle, state.stream)
        task_local_storage((:CUBLAS, state.context), (handle, state.stream, math_mode))
        stream = state.stream
    end

    if math_mode != state.math_mode
        math_mode!(handle, state.math_mode)
        task_local_storage((:CUBLAS, state.context), (handle, stream, state.math_mode))
        math_mode = state.math_mode
    end

    return handle
end

function xt_handle()
    ctxs = Tuple(context(dev) for dev in devices())
    get!(task_local_storage(), (:CUBLASxt, ctxs)) do
        handle = @lock handle_cache_lock begin
            if isempty(idle_xt_handles[ctxs])
                cublasXtCreate()
            else
                pop!(idle_xt_handles[ctxs])
            end
        end

        finalizer(current_task()) do task
            @spinlock handle_cache_lock begin
                push!(idle_xt_handles[ctxs], handle)
            end
        end
        # TODO: cublasXtDestroy to preserve memory, or at exit?

        # select all devices
        devs = convert.(Cint, devices())
        cublasXtDeviceSelect(handle, length(devs), devs)

        handle
    end::cublasHandle_t
end

function log_message(cstr)
    # NOTE: we can't `@debug` these messages, because the logging function is called for
    #       every line... also not sure what the i!/I! prefixes mean (info?)
    # NOTE: turns out we even can't `print` these messages, as cublasXt callbacks might
    #       happen from a different thread! we could strdup + uv_async_send, but I couldn't
    #       find an easy way to attach data to that
    # TODO: use a pre-allocated lock-free global message buffer?
    len = ccall(:strlen, Csize_t, (Cstring,), cstr)
    ccall(:write, Cint, (Cint, Cstring, Csize_t), 0, cstr, len)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUBLAS
    if isdebug(:init, CUBLAS)
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end


text:
__precompile__()

"""
# SeisRequests

Basic implementation of seismic web requests, returning only the raw results
of a request.

Currently implemented schemes are the FDSN Web Services specification, and IRIS's
web services specification for data (timeseries).  Further specifications will
be added in the future.

To request data, create a one of the [FDSNRequest](@ref) types and use the
[get_request](@ref) method on it.  This method returns a `HTTP.Message.Response`
containing the raw response in the field `body`.

### Specifications

FDSN Web Services request specification:
    http://www.fdsn.org/webservices/FDSN-WS-Specifications-1.1.pdf

Details on IRISWS timeseries can be found at:
    https://service.iris.edu/irisws/timeseries/1/
"""
module SeisRequests

export
    FDSNDataSelect,
    FDSNEvent,
    FDSNStation,
    IRISTimeSeries,
    add_server!,
    get_request,
    request_uri,
    server_list

using Compat
using Compat.Dates
import Compat: @info
@static if VERSION < v"0.7-"
    using Missings
end

using HTTP
using Parameters
import DataStructures: OrderedDict

"URIs of servers to which requests can be sent by key"
const SERVERS = Dict{String,String}(
    "IRIS" => "http://service.iris.edu",
    "INGV" => "http://webservices.ingv.it",
    "Orfeus" => "http://www.orfeus-eu.org",
    "NCEDC" => "http://service.ncedc.org")

const DEFAULT_SERVER = "IRIS"

"Dict describing the common HTTP status codes returned by FDSN services"
const STATUS_CODES = Dict{Int,String}(
    200 => "Successful request, results follow",
    204 => "Request was properly formatted and submitted but no data matches the selection",
    400 => "Bad request due to improper specification, unrecognized parameter, parameter value out of range, etc.",
    401 => "Unauthorized, authentication required",
    403 => "Authentication failed or access blocked to restricted data",
    404 => "Alternate to 204 (set via the ‘nodata’ parameter), normally used for results returned to a web browser.",
    413 => "Request would result in too much data being returned or the request itself is too large, returned error message should include the service limitations in the detailed description. Service limits should also be documented in the service WADL.",
    414 => "Request URI too large",
    500 => "Internal server error",
    503 => "Service temporarily unavailable, used in maintenance and error conditions")
const CODE_SUCCESS = 200
const CODES_FAILURE = (400, 401, 403, 404, 413, 414, 500, 503)

# Types which accept missing values
const MBool = Union{Missing,Bool}
const MDateTime = Union{Missing,DateTime}
const MFloat = Union{Missing,Float64}
const MInt = Union{Missing,Int}
const MString = Union{Missing,String}

"""
    SeisRequest

An abstract type representing any kind of web request for seismic data.

Current subtypes of `SeisRequest`:
- `FDSNRequest`
- `IRISRequest`
"""
abstract type SeisRequest end

include("fdsnws.jl")
include("irisws.jl")

"""Default version string for all requests."""
version_string(::SeisRequest) = "1"

"""
    request_uri(r::SeisRequest; server=$(DEFAULT_SERVER)) -> uri

Return a URI for the request `r`, which can then be obtained via HTTP GET.

`server` can either be a URI or one of the available servers.  (See [server_list](@ref).)
"""
function request_uri(r::SeisRequest; server=DEFAULT_SERVER)
    server = server in keys(SERVERS) ? SERVERS[server] : server
    protocol = protocol_string(r)
    service = service_string(r)
    version = version_string(r)
    uri = join((server, protocol, service, version, "query?"), "/")
    firstfield = true

    if server == "http://service.ncedc.org"
        for f in fieldnames(typeof(r))
            v = getfield(r, f)
            if !ismissing(v)
                if f == :starttime
                    firstfield ? (uri *= "$(f)"[1:5]*"=$(v)"; firstfield=false) : uri = join((uri, "$(f)"[1:5]*"=$(v)"), "&")

                elseif (f == :format) || (f == :nodata) #not print in uri

                else
                    firstfield ? (uri *= "$(f)"[1:3]*"=$(v)"; firstfield=false) : uri = join((uri, "$(f)"[1:3]*"=$(v)"), "&")
                end
            end
        end

    else
        for f in fieldnames(typeof(r))
            v = getfield(r, f)
            if !ismissing(v)
                if firstfield
                    uri *= "$(f)=$(v)"
                    firstfield = false
                else
                    uri = join((uri, "$(f)=$(v)"), "&")
                end
            end
        end
    end

    uri
end

"""
    add_server!(label, uri)

Add a server with key `label` to the global list of servers.

See also [server_list](@ref).
"""
function add_server!(label, uri)
    label ∈ keys(SERVERS) && throw(ArgumentError("A server with label `$label` already exists"))
    SERVERS[label] = uri
end

"""
    server_list() -> servers

Return a list of available servers for use with a `SeisRequest`.

See also [add_server!](@ref).
"""
server_list() = collect(keys(SERVERS))

"""
    get_request(r::SeisRequest; server=$(DEFAULT_SERVER)) -> response::HTTP.Message.Response

Return `response`, the result of requesting `r` via the HTTP GET command.  Optionally specify
the server either by URI or one of the available servers.  (See [server_list](@ref).)
"""
function get_request(r::SeisRequest; server=DEFAULT_SERVER, verbose=true)
    uri = request_uri(r; server=server)

    #Specify request headers
    #Wrong request headers cause server internal error (500)
    if server == "NCEDC"
        headers = ["Host" => "service.ncedc.org", "User-Agent" => "curl/7.60.0", "Accept" => "*/*"]
    elseif server == "IRIS"
        headers = ["Host" => "service.iris.edu", "User-Agent" => "curl/7.60.0", "Accept" => "*/*"]
    else
        headers = []
    end

    response = HTTP.request("GET", uri; headers=headers)
    status_text = STATUS_CODES[response.status]
    if verbose
        @info("Request status: " * status_text)
    else
        response.status in CODES_FAILURE && warn("Request status: " * status_text)
    end
    response
end

end # module


text:
using LidJul
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using GLMakie
using IterativeSolvers
using Preconditioners
using Random

function addmeasurements(msm,solvername,ti,ts,res,iterations=[])
    push!(msm,(solver=solvername,init_time=ti,solver_time=ts,residual=res,iterations=iterations))
end

function resnorm(splu,pref,b)
    ax=similar(b)
    LidJul.mul!(ax,splu,pref)
    norm(ax-b)
end

function test_solver!(msm,Solver,solver_args,splu,b)
    p=zeros(size(b))
    init_time=@elapsed solver=Solver(solver_args...)
    init_time=@elapsed solver=Solver(solver_args...)
    Random.seed!(1234)
    prand=rand(length(p))
    iterations=begin copyto!(p,prand) ; solve!(p,b,solver)  end
    solve_time=@elapsed  begin copyto!(p,prand) ; solve!(p,b,solver) end
    sname=string(Solver)
    length(solver_args)>1 && (sname=sname*"_"*string(solver_args[2]))

    addmeasurements(msm,sname,init_time,solve_time,resnorm(splu,p,b),iterations)
    p
end

function precoconstructor(A,preco)
    if preco==CholeskyPreconditioner
        return preco(A,4)
    elseif preco===nothing
        return IterativeSolvers.Identity()
    else
        return preco(A)
    end
end

function test_pcg(msm,preco,A,splu,Bxy)
    b=view(Bxy,1:length(Bxy))
    p=zeros(size(b))
    # x=view(Xxy,1:length(Xxy))
    sname="PCG"

    init_time=@elapsed pc=precoconstructor(A,preco)
    init_time=@elapsed pc=precoconstructor(A,preco)
    Random.seed!(1234)
    prand=rand(length(p))
    p,iterations = begin copyto!(p,prand) ; IterativeSolvers.cg!(p,A, b, Pl=pc;log=true,reltol = 1e-16,maxiter=2000) end
    solve_time = @elapsed begin copyto!(p,prand) ; IterativeSolvers.cg!(p,A, b, Pl=pc;log=true,reltol = 1e-16,maxiter=2000) end
    spreco=string(preco)[1:min(15,length(string(preco)))]
    sname="PCG"*spreco
    addmeasurements(msm,sname,init_time,solve_time,resnorm(splu,p,b),iterations)
    nothing
end
function test_stationary(msm,stm,A,splu,Bxy)
    spiter=PoissonSparseIterative(A)
    init_time=@elapsed PoissonSparseIterative(A)
    Xxy=similar(Bxy)
    Random.seed!(1234)
    Xxyrand=rand(size(Xxy))
    p,iterations = begin copyto!(Xxy,Xxyrand) ; solve!(Xxy,Bxy,spiter,stm) end
    solve_time = @elapsed begin copyto!(Xxy,Xxyrand) ; solve!(Xxy,Bxy,spiter,stm) end
    sname=string(stm)
    addmeasurements(msm,sname,init_time,solve_time,resnorm(splu,Xxy,Xxy),iterations)
    nothing
end

include("measureplots.jl")

function testpoisson(n,bc)
    msm=Vector{NamedTuple{(:solver,:init_time,:solver_time,:residual,:iterations),Tuple{String,Float64,Float64,Float64,Any}}}()

    p=ones(n,n)
    pref=ones(n,n)
    b=zeros(n,n)
    #LP: Create a Poisson's RHS compatible with a possible full Neumann boundary conditions
    #See bellow
    for i=1:n
        for j=1:n
            b[i,j]=sin(3π*(i-1)/(n-1))*sin(6π*(j-1)/(n-1))
        end
    end

    L=1.0
    dx=L/n
    dy=L/n
    @show dx,dy

    l2D=Laplacian2D(n,n,L,L,bc...)
    if bc==(neumann,neumann,neumann,neumann)
        #LP: for all neumann conditions the laplacian equation has no
        #solution if the rhs b does not belong to the range of A.
        #This impose the integral of b to be zero on the domain and
        # the sum of b on the mesh to be also zero...
        #See Mc Cormick and other for more details.
        @show sum(b)
        @assert(sum(b)<1.e-8)

        #Note that if this compatibility condition is ensured, and
        # x is a solution, then x+Cte is also a solution.
    end

    @time s2D=sparse_corr(l2D)
    splu,tlu=PoissonSparseLU(s2D)
    ts=@elapsed solve!(pref,b,splu)

    @time s2D_nocorr=SparseArrays.sparse(l2D)
    splu_nocorr,tlu_nocorr=PoissonSparseLU(s2D_nocorr)
    stationary_methods=["jacobi","gauss_seidel","sor","ssor"]
    for stm in stationary_methods
        test_stationary(msm,stm,s2D,splu,b)
    end

    test_pcg(msm,Preconditioners.DiagonalPreconditioner,s2D,splu,b)
    test_pcg(msm,nothing,s2D,splu,b)
    test_pcg(msm,Preconditioners.AMGPreconditioner{SmoothedAggregation},s2D,splu,b)
    test_pcg(msm,Preconditioners.AMGPreconditioner{RugeStuben},s2D,splu,b)
    test_solver!(msm,PoissonSparseAMG,(s2D,),splu,b)

    test_solver!(msm,PoissonTTSolver,(l2D,),splu,b)

    addmeasurements(msm,string(PoissonSparseLU),tlu,ts,resnorm(splu,pref,b))

    test_solver!(msm,PoissonGMG,(l2D,GSSmoother),splu_nocorr,b)

    test_solver!(msm,PoissonSparseCGILU,(s2D,),splu,b)

    sname=string(n)*tostring(bc)

    @save "msm"*sname*".jld2" msm

    measureplot(msm,sname)

    return p,pref,p,splu,b,msm

end

function tostring(bcs)
     r=string()
     for bc in bcs
         bc==dirichlet && (r*="D")
         bc==neumann && (r*="N")
    end
    r
end

function go()
    #Choose a power of two
    n=128
    #Choose boundary conditions
    # bc=(neumann,neumann,neumann,neumann)
    bc=(dirichlet,neumann,dirichlet,neumann)
    # bc=(dirichlet,dirichlet,dirichlet,dirichlet)

    p,pref,pbj,s,b,msm=testpoisson(n,bc)

   #*** Makie Plots
    x = range(0,1,length=n)
    y = range(0,1,length=n)

    p=pref

    pmax=maximum(p)-minimum(p)

    p./=pmax

    z = p # to match https://lazarusa.github.io/BeautifulMakie/surfWireLines/surfWireContour/ notation

    zmin, zmax = minimum(z), maximum(z)
    cmap = :viridis
    set_theme!(theme_dark())
    fig = Figure(resolution = (1200, 800), fontsize = 22)
    ax = Axis3(fig[1, 1], aspect = :data,
        xzpanelcolor = (:black, 0.75), yzpanelcolor = (:black, 0.75),
        zgridcolor = :grey, ygridcolor = :grey, xgridcolor = :grey)
        # ax = Axis3(fig[1, 1], aspect = :data, perspectiveness = 0.5, elevation = π / 9,
        # xzpanelcolor = (:black, 0.75), yzpanelcolor = (:black, 0.75),
        # zgridcolor = :grey, ygridcolor = :grey, xgridcolor = :grey)

    sm = surface!(ax, x, y, z; colormap = cmap, colorrange = (zmin, zmax),
        transparency = false)
    xm, ym, zm = minimum(ax.finallimits[])
    contour!(ax, x, y, z; levels = 20, colormap = cmap, linewidth = 2,
        colorrange = (zmin, zmax), transformation = (:xy, zm),
        transparency = true)
    wireframe!(ax, x, y, z; overdraw = true, transparency = true,
        color = (:black, 0.1))
    Colorbar(fig[1, 2], sm, height = Relative(0.5))
    colsize!(fig.layout, 1, Aspect(1, 1.0))
    # display(fig)
    set_theme!()

    # s=Makie.surface(x,y,p)
    # # xm, ym, zm = minimum(scene_limits(s))
    # # Makie.contour!(s,x,y,p, levels = 15, linewidth = 2, transformation = (:xy, zm))
    # Makie.contour!(s,x,y,p, levels = 15, linewidth = 2)
    # Makie.wireframe!(s,x,y,p, overdraw = true, transparency = true, color = (:black, 0.1))
    # display(AbstractPlotting.PlotDisplay(), s)
    # resize!(s,(1600,800))
    Makie.save("makie"*string(n)*tostring(bc)*".png",fig)
    nothing
end

go()


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

using Markdown
using InteractiveUtils

# ╔═╡ c7edc362-1d55-4556-94e2-5ccae0fce152
begin
	using LinearAlgebra
	using Statistics
	using PlutoUI
	using Plots
	gr()
end

# ╔═╡ e91395c4-e4dd-4bf7-8efa-c48d05418bb4
using Polynomials

# ╔═╡ f16d832c-67ad-11ec-3474-d9f0f2241ff9
html"""

Happy holiday! Remember to take care of yourself and your loved ones!

"""

# ╔═╡ f16d834c-67ad-11ec-05db-f726c8f8c643
md""" # 機械学習帳
ここを参考にした.[機械学習帳](https://chokkan.github.io/mlnote/regression/01sra.html)
## 1. 単回帰
### (続き)1.4. 平均二乗誤差は$a$,$b$の関数
"""

# ╔═╡ 6ce3e88a-e10d-45c2-9722-aeed933931d9
begin
	D = [
		1 3
		3 6
		6 5
		8 7]

	X = D[:,1]
	Y = D[:,2]
end;

# ╔═╡ bc133d95-c122-49be-8789-01dd2b87ac0a
function calc_msr(a, b, D)
	xm = mean(D[:,1])
	ym = mean(D[:,2])
	x2m = mean(D[:,1].^2)
	y2m = mean(D[:,2].^2)
	xym = mean(D[:,1] .* D[:,2])
	return @. x2m * a^2 + b^2 + 2*xm *a *b - 2*xym *a - 2*ym*b + y2m
end

# ╔═╡ 86db1e1c-5164-4e3f-b877-f9c6f1a535c7
begin
	N = 100
	aa = range(-10,10, N)
	bb = range(-10,10, N)
end

# ╔═╡ c47e13d4-f0e9-46b0-b746-73563344487e
J = calc_msr(aa, bb', D)

# ╔═╡ f1896c70-c8ce-4d7f-832e-11a3eede0883
let
	plot(xlim=(-10,10), ylim=(-10,10), framestyle=:box, aspect_ratio=1)
	plot!(xticks=-10:5:10, yticks=-10.0:2.5:10.0)
	plot!(xlabel="a", ylabel="b", zlabel="L̂_D(a,b)")
	heatmap!(aa, bb, J'; colorbar_scale=:log10, clim=(10^(-1),10^4),
		fillcolor=cgrad(:viridis) )
	J_levels = [1., 10.0, 100.0, 1000.0]
	contour!(aa, bb, J'; lc=:red, lw=0.5, ls=:dash, levels= J_levels)
end

# ╔═╡ 02c50675-5c7e-4662-84fd-56cf5e34e666
md""" ### 1.5. 平均二乗残差を最小化する
単回帰の目的関数(平均二乗残差)は
```math
L̂_D(a,b) = \frac{N}{1} \sum_{i=1}^{N} ϵ_i^2 = \frac{1}{N} \sum_{i=1}^{N} (y_i - ŷ_i)^2 = \frac{1}{N} \sum_{i=1}^{N} (y_i - ax_i - b)^2
```
である。
上式を $b$ で偏微分すると
```math
b = ȳ - ax̄
```
となる。 単回帰の傾き $a$ を求める式
```math
a = \frac{ \overline{xy} - (x̄)(ȳ)}{ \overline{x^2} - (x̄)^2}
	= \frac{  \rm{Cov[} \it{X, Y} \rm{]} } {\rm{Cov[} \it{X, X} \rm{]}}
    = \frac{  \rm{Cov[} \it{X, Y} \rm{]} } {\rm{Var[} \it{X} \rm{]}}
```

"""

# ╔═╡ 0367c8ba-1c50-43ea-828f-fb0cbaefe691
function set_plot()
	plot(size=(500,500), xlim=(0,10), ylim=(0,10), framestyle=:box, legend=:topleft)
	plot!(xlabel="x", ylabel="y", gridalpha=0.5)
end

# ╔═╡ 4765fa67-4e96-4cc0-85fa-266128fb8fae
function plot_y(X, Y, colormap)
	scatter!(X, Y; c=colormap,  markerstrokewidth=0, label="")
end

# ╔═╡ 8f419e32-93f9-4067-9461-e0e2f69a723a
function plot_ŷ(a, b, offset, colormap)

	offset_ŷ = offset

	a, b = a, b
	f(x, a, b) = a*x + b
	Ŷ = f.(X, a, b)
	x = [0, 10]

	a = round(a; digits=3)
	b = round(b; digits=3)

	label = (a == 1) ? "x + $(b))" : "$(a)x + $(b)"
	scatter!(D[:,1], Ŷ; c=colormap, shape=:star5, ms=6, markerstrokewidth=0, label="")
	plot!(x -> f(x, a, b); c=:black, lw=2, label=label)

	texts_ŷ = ["(x₁,ŷ₁)", "(x₂,ŷ₂)", "(x₃,ŷ₃)", "(x₄,ŷ₄)"]
	for i in 1:length(colormap)
		annotate!([(X[i],Ŷ[i]+offset_ŷ[i],text(texts_ŷ[i], 10,
			:center, colormap[i]))])
	end
	plot!()
end

# ╔═╡ 59bd7898-2dcf-4c23-815b-abf13a6954d0
function plot_error(a, b, offset, colormap)

	f(x, a, b) = a*x + b
	Ŷ = f.(X, a, b)

	offset_ϵ = offset
	x_ϵ = X + offset_ϵ
	y_ϵ = (Y + Ŷ) / 2
	texts_ϵ = ["ϵ₁", "ϵ₂", "ϵ₃", "ϵ₄"]

	for i in 1:length(colormap)
		plot!([(X[i], Ŷ[i]), (X[i], Y[i])]; c=colormap[i], ls=:dash, label="")
		annotate!([(x_ϵ[i],y_ϵ[i],text(texts_ϵ[i], 10, :center, colormap[i]))])
	end
	plot!()
end

# ╔═╡ e776c006-777f-4a6c-8a3e-a273820a0ad8
const colormap = [:red, :blue, :green, :orange]

# ╔═╡ 75364a95-c8ef-40d3-8373-6176540ceaf8
let
	# パラメータ a, b の計算
	x̄ = D[:,1] |> mean
	ȳ = D[:,2] |> mean

	Cov_XY = 1 / length(X) * sum( (X .- x̄) .* (Y .- ȳ) )
	Var_X =  1 / length(X) * sum( (X .- x̄).^2 )

	a = Cov_XY / Var_X
	b = ȳ - a* x̄

	# 可視化
	offset_y_new = [-1., 1.8, -1.2, 0.7]
	offset_ϵ_new = [0.3, -0.3, 0.3, -0.3]

	set_plot()
	plot_y(X, Y, colormap)
	plot_ŷ(a, b, offset_y_new, colormap)
	plot_error(a, b, offset_ϵ_new, colormap)

end

# ╔═╡ 792c62f0-222b-4048-ad6a-9655db74f080
md""" ### 1.6. 回帰直線の性質
"""

# ╔═╡ 0f1d5e27-f04c-47f4-bcff-bc9deea6eda2
function plot_means(X, Y)
	x̄ = D[:,1] |> mean
	ȳ = D[:,2] |> mean

	vline!([x̄]; c=:black, lw=1.5, ls=:dash, label="x̄")
	hline!([ȳ]; c=:black, lw=1.5, ls=:dot, label="ȳ")

end

# ╔═╡ 4c2841e2-2e5e-47ca-987a-0e4c5557f023
function plot_regions(X, Y)
	x̄ = D[:,1] |> mean
	ȳ = D[:,2] |> mean

	rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
	plot!(rectangle(0, 0, x̄, ȳ), color=:red, alpha=0.1, label="")
	plot!(rectangle(x̄, ȳ, 10-x̄, 10-ȳ), color=:red, alpha=0.1, label="")

	plot!(rectangle(0, ȳ, x̄, ȳ), color=:blue, alpha=0.1, label="")
	plot!(rectangle(x̄, 0, 10-x̄, ȳ), color=:blue, alpha=0.1, label="")

end

# ╔═╡ dde67d9a-5c75-4477-bbd0-3eefc7b5ef13
let
	set_plot()
	plot_y(X, Y, colormap)
	plot_means(X, Y)
	plot_regions(X, Y)
end

# ╔═╡ c0a1bf99-d350-47ee-9a01-763acdcdb63a
function plot_regions2(X, Y, colormap)
	x̄ = D[:,1] |> mean
	ȳ = D[:,2] |> mean

	rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

	xs = [D[1,1], D[2,1], x̄, x̄]
	ys = [D[1,2], D[2,2], D[3,2], ȳ]
	ws = [x̄-D[1,1], x̄-D[2,1], D[3,1]-x̄, D[4,1]-x̄]
	hs = [ȳ-D[1,2], ȳ-D[2,2], ȳ-D[3,2], D[4,2]-ȳ]

	for i in 1:length(xs)
		plot!(rectangle(xs[i], ys[i], ws[i], hs[i]), color=colormap[i],
			alpha=0.1, label="")
	end
	plot!()
end

# ╔═╡ 407a1c1b-bb2d-4188-997b-17565ef4431d
let
	set_plot()
	plot_y(X, Y, colormap)
	plot_means(X, Y)
	plot_regions2(X, Y, colormap)
end

# ╔═╡ cc5a8abc-ff3e-4693-94fe-79a73aaf2a19
md""" #### 1.6.3. 回帰直線は観測データの重心を通る
"""

# ╔═╡ 6da0bfa7-bd9e-4463-9057-bd40c56019eb
function plot_slope(x̄, ȳ, a, b)

	Cov_XY = 1 / length(X) * sum( (X .- x̄) .* (Y .- ȳ) )
	Var_X =  1 / length(X) * sum( (X .- x̄).^2 )

	a = Cov_XY / Var_X
	b = ȳ - a* x̄

	# horizontal
	plot!([0, x̄],[b, b]; lw=2, color=:brown, label="")
	annotate!([(x̄/2,b-0.3,text("x̄", 10, :center, color=:brown))])
	# vertical
	plot!([x̄, x̄], [b, ȳ]; lw=2, color=:brown, label="")
	annotate!([(x̄+0.3, b+(ȳ-b)/2,text("ax̄", 10, :center, color=:brown))])
end

# ╔═╡ f4215be2-c8af-4f7e-92cb-76b338714f1a
function plot_ŷ_only(a, b)

	a, b = a, b
	f(x, a, b) = a*x + b
	Ŷ = f.(X, a, b)
	x = [0, 10]

	a = round(a; digits=3)
	b = round(b; digits=3)

	label = (a == 1) ? "x + $(b))" : "$(a)x + $(b)"
	plot!(x -> f(x, a, b); c=:black, lw=2, label=label)

end

# ╔═╡ c1635018-bff4-4245-91ff-278fda46f3a0
let

	# パラメータ a, b の計算
	x̄ = D[:,1] |> mean
	ȳ = D[:,2] |> mean

	Cov_XY = 1 / length(X) * sum( (X .- x̄) .* (Y .- ȳ) )
	Var_X =  1 / length(X) * sum( (X .- x̄).^2 )

	a = Cov_XY / Var_X
	b = ȳ - a* x̄

	# ----- 可視化 -----
	# 初期化
	set_plot()
	# データ点をプロット
	plot_y(X, Y, colormap)
	# xとyの平均の直線をプロット
	plot_means(X, Y)
	# +(赤)、 -(青)になる領域をプロット
	plot_regions(X, Y)
	# スロープ部分(茶色)をプロット
	plot_slope(x̄, ȳ, a, b)
	# 推定値 ŷ をプロット
	plot_ŷ_only(a, b)
end

# ╔═╡ fcd91d11-e8e0-4419-a992-47f82b39b393
md""" #### 1.8. 単回帰の実施例
"""

# ╔═╡ e2ae904b-c4aa-4f3c-ad18-58b0fabb6877
begin
	mmt = [
		9.1, 11.2, 12.3, 18.9, 22.2, 26. , 30.9, 31.2, 28.8, 23. , 18.3,
       11.1,  8.3,  9.1, 12.5, 18.5, 23.6, 24.8, 30.1, 33.1, 29.8, 23. ,
       16.3, 11.2,  9.6, 10.3, 16.4, 19.2, 24.1, 26.5, 31.4, 33.2, 28.8,
	   23. , 17.4, 12.1, 10.6,  9.8, 14.5, 19.6, 24.7, 26.9, 30.5, 31.2,
	   26.9, 23. , 17.4, 11. , 10.4, 10.4, 15.5, 19.3, 26.4, 26.4, 30.1,
	   30.5, 26.4, 22.7, 17.8, 13.4, 10.6, 12.2, 14.9, 20.3, 25.2, 26.3,
	   29.7, 31.6, 27.7, 22.6, 15.5, 13.8, 10.8, 12.1, 13.4, 19.9, 25.1,
	   26.4, 31.8, 30.4, 26.8, 20.1, 16.6, 11.1,  9.4, 10.1, 16.9, 22.1,
	   24.6, 26.6, 32.7, 32.5, 26.6, 23. , 17.7, 12.1, 10.3, 11.6, 15.4,
	   19. , 25.3, 25.8, 27.5, 32.8, 29.4, 23.3, 17.7, 12.6, 11.1, 13.3,
	   16. , 18.2, 24. , 27.5, 27.7, 34.1, 28.1, 21.4, 18.6, 12.3
	]

	spend = [
	  463.,  360.,  380.,  584.,  763.,  886., 1168., 1325.,  847.,
      542.,  441.,  499.,  363.,  327.,  414.,  545.,  726.,  847.,
     1122., 1355.,  916.,  571.,  377.,  465.,  377.,  362.,  518.,
      683.,  838., 1012., 1267., 1464., 1000.,  629.,  448.,  466.,
      404.,  343.,  493.,  575.,  921., 1019., 1149., 1303.,  805.,
      739.,  587.,  561.,  486.,  470.,  564.,  609.,  899.,  946.,
     1295., 1325.,  760.,  667.,  564.,  633.,  478.,  450.,  567.,
      611.,  947.,  962., 1309., 1307.,  930.,  668.,  496.,  650.,
      506.,  423.,  531.,  672.,  871.,  986., 1368., 1319.,  924.,
      716.,  651.,  708.,  609.,  535.,  717.,  890., 1054., 1077.,
     1425., 1378.,  900.,  725.,  554.,  542.,  561.,  459.,  604.,
      745., 1105.,  973., 1263., 1533., 1044.,  821.,  621.,  601.,
      549.,  572.,  711.,  819., 1141., 1350., 1285., 1643., 1133.,
      784.,  682.,  587.
	]
end;

# ╔═╡ 37558fb1-0ea2-4aa5-bdae-d1b54d237812
function plot_temp_vs_spend()
	plot(title="MMT vs. Spend for Icecream", framestyle=:box)
	plot!(xlabel="month mean maximum temperature [̇Degree]", ylabel="spend [Yen]")
	plot!(xlim=(0,35), ylim=(-250,2000), yticks=-250:250:2000, gridalpha=0.5)
	scatter!(mmt, spend; marker=:dot, ms=3, label="")
end

# ╔═╡ 569e9185-fb79-4a0b-b32d-b6e20dee29e8
let
	ŷ = fit(mmt, spend, 1)
	x = range(0, 35, length=100)
	plot_temp_vs_spend()
	plot!(x, ŷ.(x); lw=2, color=:red, label="")
end

# ╔═╡ e802792b-052f-4bd5-976a-2f5983fe7e06

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
Plots = "~1.25.3"
PlutoUI = "~0.7.27"
Polynomials = "~2.0.22"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[deps.Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.2"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[deps.ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[deps.FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[deps.FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[deps.FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[deps.GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "0c603255764a1fa0b61752d2bec14cfbd18f7fe8"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+1"

[[deps.GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "30f2b340c2fff8410d89bfcdc9c0a6dd661ac5f7"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.62.1"

[[deps.GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f97acd98255568c3c9b416c5a3cf246c1315771b"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.63.0+0"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[deps.Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[deps.Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+2"

[[deps.Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[deps.Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "8d70835a3759cdd75881426fced1508bb7b7e1b6"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.1"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.Intervals]]
deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
git-tree-sha1 = "323a38ed1952d30586d0fe03412cde9399d3618b"
uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
version = "1.5.0"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[deps.LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[deps.LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a8f4f279b6fa3c3c4f1adadd78a621b13a506bce"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.9"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[deps.Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[deps.Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[deps.Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[deps.Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[deps.Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[deps.Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "73deac2cbae0820f43971fad6c08f6c4f2784ff2"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.2"

[[deps.NaNMath]]
git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.6"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[deps.Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.3"

[[deps.Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[deps.PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "68604313ed59f0408313228ba09e79252e4b2da8"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.1.2"

[[deps.Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"]
git-tree-sha1 = "7eda8e2a61e35b7f553172ef3d9eaa5e4e76d92e"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.25.3"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "fed057115644d04fba7f4d768faeeeff6ad11a60"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.27"

[[deps.Polynomials]]
deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"]
git-tree-sha1 = "a2515a5baa0b6954e96e10a8dc3f6855b4dedb97"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
version = "2.0.22"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.4.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.2.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "d88665adc9bcf45903013af0982e2fd05ae3d0a6"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2bb0cb32026a66037360606510fca5984ccc6b75"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.13"

[[deps.StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TimeZones]]
deps = ["Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "ce5aab0b0146b81efefae52f13002e19c2af57ac"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.7.0"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[deps.Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

[[deps.Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[deps.Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "66d72dc6fcc86352f01676e8f0f698562e60510f"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.23.0+0"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[deps.XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[deps.Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[deps.Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[deps.Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[deps.Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[deps.Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[deps.Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[deps.Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[deps.Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[deps.Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[deps.Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[deps.Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[deps.Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[deps.Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[deps.Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[deps.Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[deps.Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[deps.Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[deps.Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[deps.libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[deps.libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[deps.libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[deps.x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[deps.x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[deps.xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─f16d832c-67ad-11ec-3474-d9f0f2241ff9
# ╠═c7edc362-1d55-4556-94e2-5ccae0fce152
# ╠═f16d834c-67ad-11ec-05db-f726c8f8c643
# ╠═6ce3e88a-e10d-45c2-9722-aeed933931d9
# ╠═bc133d95-c122-49be-8789-01dd2b87ac0a
# ╠═86db1e1c-5164-4e3f-b877-f9c6f1a535c7
# ╠═c47e13d4-f0e9-46b0-b746-73563344487e
# ╠═f1896c70-c8ce-4d7f-832e-11a3eede0883
# ╠═02c50675-5c7e-4662-84fd-56cf5e34e666
# ╠═0367c8ba-1c50-43ea-828f-fb0cbaefe691
# ╠═4765fa67-4e96-4cc0-85fa-266128fb8fae
# ╠═8f419e32-93f9-4067-9461-e0e2f69a723a
# ╠═59bd7898-2dcf-4c23-815b-abf13a6954d0
# ╠═e776c006-777f-4a6c-8a3e-a273820a0ad8
# ╠═75364a95-c8ef-40d3-8373-6176540ceaf8
# ╠═792c62f0-222b-4048-ad6a-9655db74f080
# ╠═0f1d5e27-f04c-47f4-bcff-bc9deea6eda2
# ╠═4c2841e2-2e5e-47ca-987a-0e4c5557f023
# ╠═dde67d9a-5c75-4477-bbd0-3eefc7b5ef13
# ╠═c0a1bf99-d350-47ee-9a01-763acdcdb63a
# ╠═407a1c1b-bb2d-4188-997b-17565ef4431d
# ╠═cc5a8abc-ff3e-4693-94fe-79a73aaf2a19
# ╠═6da0bfa7-bd9e-4463-9057-bd40c56019eb
# ╠═f4215be2-c8af-4f7e-92cb-76b338714f1a
# ╠═c1635018-bff4-4245-91ff-278fda46f3a0
# ╠═fcd91d11-e8e0-4419-a992-47f82b39b393
# ╠═e2ae904b-c4aa-4f3c-ad18-58b0fabb6877
# ╠═37558fb1-0ea2-4aa5-bdae-d1b54d237812
# ╠═e91395c4-e4dd-4bf7-8efa-c48d05418bb4
# ╠═569e9185-fb79-4a0b-b32d-b6e20dee29e8
# ╠═e802792b-052f-4bd5-976a-2f5983fe7e06
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: harakiso/pluto

max_issues_repo_name: harakiso/pluto

max_forks_repo_name: harakiso/pluto


text:
module TileMaps

import Crayons

include("objects.jl")
include("object_indexable_array.jl")
include("tile_map.jl")
include("visualization.jl")

end


text:
if isdefined(:_HAVE_GSL)
N=5
A=randn(N,N)
B=randn(N,N)
Q=UniformHaar(2, N)

#Test case where there are no symbolic Haar matrices
@test_approx_eq eval(expectation(N, :Q, :(A*B))) A*B
#Test case where there is one pair of symbolic Haar matrices
@test_approx_eq eval(expectedtrace(N, :Q, :(A*Q*B*Q'))) trace(A)*trace(B)/N

println("Case 3")
println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))

end #_HAVE_GSL


text:
function paraview_collection(filename_noext::AbstractString)
    # Initialise Paraview collection file (extension .pvd).
    # filename_noext: filename without the extension (.pvd).
    xpvd = XMLDocument()
    xroot = create_root(xpvd, "VTKFile")
    set_attribute(xroot, "type",       "Collection")
    set_attribute(xroot, "version",    "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    new_child(xroot, "Collection")
    return CollectionFile(xpvd, string(filename_noext, ".pvd"))
end

function collection_add_timestep(pvd::CollectionFile, datfile::VTKFile,
                                 t::AbstractFloat)
  xroot = root(pvd.xdoc)
  xMBDS = find_element(xroot, "Collection")
  xDataSet = new_child(xMBDS, "DataSet")
  fname = splitdir(datfile.path)[2]
  set_attribute(xDataSet, "timestep", @sprintf("%f", t))
  set_attribute(xDataSet, "part", "0")
  set_attribute(xDataSet, "file", fname)
  push!(pvd.timeSteps, datfile)
  return
end

function vtk_save(pvd::CollectionFile)
    # Saves paraview collection file (.pvd).
    # Also saves the contained data files recursively.
    outfiles = [pvd.path]::Vector{UTF8String}
    for vtk in pvd.timeSteps
        append!(outfiles, vtk_save(vtk))
    end
    if isopen(pvd)
        save_file(pvd.xdoc, pvd.path)
        close(pvd)
    end
    return outfiles
end


text:
"""
```
segments                = meanshift(img, spatial_radius, range_radius; iter=50, eps=0.01))
```
Segments the image using meanshift clustering. Returns a `SegmentedImage`.

Parameters:
-    img                            = input grayscale image
-    spatial_radius/range_radius    = bandwidth parameters of truncated normal kernel. Controlling the size of the kernel determines the resolution of the mode detection.
-    iter/eps                       = stopping criterion for meanshift procedure. The algorithm stops after iter iterations or if kernel center moves less than eps distance in an update step, whichever comes first.
"""
function meanshift(img::Array{CT, 2}, spatial_radius::Real, range_radius::Real; iter::Int = 50, eps::Real = 0.01) where CT

    rows, cols = size(img)
    rowbins = Int(floor(rows/spatial_radius))
    colbins = Int(floor(cols/spatial_radius))
    colorbins = Int(floor(1/range_radius))

    colorbin(x)::Int = min(floor(x/range_radius) + 1, colorbins)
    rowbin(x)::Int = min(floor(x/spatial_radius) +1, rowbins)
    colbin(x)::Int = min(floor(x/spatial_radius) +1, colbins)

    buckets = Array{Vector{CartesianIndex{2}}}(undef, rowbins, colbins, colorbins)
    for i in CartesianIndices(size(buckets))
        buckets[i]=Array{CartesianIndex{2}}(undef, 0)
    end

    for i in CartesianIndices(size(img))
        push!( buckets[rowbin(i[1]), colbin(i[2]), colorbin(img[i])], i)
    end

    function dist2(a::SVector, b::SVector)::Float64
        return ((a[1]-b[1])/spatial_radius)^2 + ((a[2]-b[2])/spatial_radius)^2 + ((a[3]-b[3])/range_radius)^2
    end

    function neighborhood_mean(pt::SVector{3, Float64})::SVector{3, Float64}
        den = 0.0
        num = SVector(0.0, 0.0, 0.0)

        R = CartesianIndices(size(buckets))
        I1, Iend = first(R), last(R)
        I = CartesianIndex(rowbin(pt[1]), colbin(pt[2]), colorbin(pt[3]))

        for J in _colon(max(I1, I-_oneunit(I)), min(Iend, I+_oneunit(I)))
            for point in buckets[J]
                if dist2(pt, SVector(point[1], point[2], img[point])) <= 1
                    num += SVector(point[1], point[2], img[point])
                    den += 1
                end
            end
        end

        return den<=0 ? pt : num/den
    end

    modes = Array{Float64}(undef, 3, rows*cols)

    for i in CartesianIndices(size(img))
        pt::SVector{3, Float64} = SVector(i[1], i[2], img[i])
        for j in 1:iter
            nextpt = neighborhood_mean(pt)
            if dist2(pt, nextpt) < eps^2
                break
            end
            pt = nextpt
        end
        modes[:, i[1] + (i[2]-1)*rows] = [pt[1]/spatial_radius, pt[2]/spatial_radius, pt[3]/range_radius]
    end

    clusters = dbscan(modes, 1.414)
    num_segments = length(clusters)
    TM = meantype(CT)
    result              = similar(img, Int)
    labels              = Array(1:num_segments)
    region_means        = Dict{Int, TM}()
    region_pix_count    = Dict{Int, Int}()

    cluster_idx = 0
    for cluster in clusters
        cluster_idx += 1
        region_pix_count[cluster_idx] = cluster.size
        for index in [cluster.core_indices; cluster.boundary_indices]
            i, j = (index-1)%rows + 1, floor(Int, (index-1)/rows) + 1
            result[i, j] = cluster_idx
            region_means[cluster_idx] = get(region_means, cluster_idx, zero(TM)) + img[i, j]
        end
        region_means[cluster_idx] /= region_pix_count[cluster_idx]
    end

    return SegmentedImage(result, labels, region_means, region_pix_count)
end


text:
import FiniteMesh

cells, points = FiniteMesh.read_mesh("../mesh/SU2/naca0012.su2")

FiniteMesh.add_group!(cells, "../mesh/SU2/naca0012.su2")

cells, points = FiniteMesh.read_mesh("../test/square.msh")


text:
using Makie

fig = Figure(resolution = (1200, 900))

ax = Axis(fig[1, 1])
ax2 = Axis(fig[2, 1])

sl_x = Slider(fig[3, 1], range = 1:length(xs), startvalue = 1)

r = 16384-100:16384+100
line = @lift xs[$(sl_x.value)][r]
line2 = @lift xs2[$(sl_x.value)][r]

lines!(ax, r, line, color = :blue)
lines!(ax2, r, line2, color = :black, linewidth=0.5)

limits!(ax, extrema(r)..., -20, 20)
limits!(ax2, extrema(r)..., -20, 20)

fig


text:
export local_update!, SW_update!, Wolff_update!, loop_update!

@doc """
    local_update!(model, param)

Updates spin configuration by local spin flip and Metropolice algorithm
"""
@inline function local_update!(model::Model, param::Parameter)
    p = convert_parameter(model, param)
    return local_update!(model, p...)
end

@doc doc"""
Information of clusters in Swendsen-Wang algorithm.

# Fields
- `activated_bonds` : The number of activated (connected) bonds of each cluster.
- `clustersize` : The number of sites in each cluster.
- `clusterspin` : Spin variable of each cluster (e.g., 1 or -1 for `Ising`).
"""
mutable struct SWInfo
    activated_bonds :: Vector{Int}
    clustersize :: Vector{Int}
    clusterspin :: Vector{Int}
end
numclusters(sw::SWInfo) = length(sw.clustersize)

@doc """
    SW_update!(model, param::Parameter)

Updates spin configuration by Swendsen-Wang algorithm
"""
@inline function SW_update!(model::Model, param::Parameter)
    p = convert_parameter(model, param)
    return SW_update!(model, p...)
end

@doc """
    Wolff_update!(model, param::Parameter)

Updates spin configuration by Wolff algorithm
"""
@inline function Wolff_update!(model::Model, param::Parameter)
    p = convert_parameter(model, param)
    return Wolff_update!(model, p...)
end

@doc """
    loop_update!(model, param::Parameter)

Updates spin configuration by loop algorithm
"""
@inline function loop_update!(model::Model, param::Parameter)
    p = convert_parameter(model, param)
    return loop_update!(model, p...)
end


text:
samplingparameter(samplingstyle::SamplingStyle, ap::ApproximationProblem; samplingstrategy=DefaultSamplingStrategy(), options...)=
    samplingparameter(samplingstyle, samplingstrategy, ap; options...)

samplingparameter(samplingstyle::GenericSamplingStyle, ap::ApproximationProblem; options...) = -1

function samplingparameter(samplingstyle::SamplingStyle, samplingstrategy::SamplingStrategy, ap::ApproximationProblem;
        verbose=false, forget=false, options...)
    if haskey(options, :L)
        # The user specified L as an option
        L = options[:L]
        setsamplingparam!(ap, L)
    elseif samplingparam(ap) != nothing && !forget
        # Has it been computed before
        L = samplingparam(ap)
    else
        # It wasn't. We deduce its value from the options given and store the outcome.
        L = deduce_samplingparameter(samplingstyle, samplingstrategy, ap; verbose=verbose, options...)
        setsamplingparam!(ap, L)
    end
    verbose && println("Sampling parameter: using L = $L")
    return L
end

samplingparameter(samplingstyle::DiscreteGramStyle, samplingstrategy::SamplingStrategy, ap::ApproximationProblem; options...) =
    samplingparameter(OversamplingStyle(), samplingstrategy, ap; options...)

deduce_samplingparameter(ss::SamplingStyle, samplingstrategy::SamplingStrategy, ap::ApproximationProblem; options...) =
    deduce_samplingparameter(ss, samplingstrategy, platform(ap), parameter(ap); options...)

deduce_samplingparameter(ss::ProductSamplingStyle, samplingstrategy::SamplingStrategy, ap::ProductPlatformApproximation; options...) =
    map((x,style)->samplingparameter(style, samplingstrategy, x; options...), ApproximationProblems.unsafe_elements(ap), ss.styles)

function deduce_samplingparameter(ss::ProductSamplingStyle{NTuple{N,OversamplingStyle}}, samplingstrategy::SamplingStrategy, ap::ProductPlatformApproximation;
        dict=dictionary(ap), oversamplingfactor=2, verbose=false, options...) where N
    oversamplingfactor = oversamplingfactor_float2tuple(samplingstrategy, oversamplingfactor, dict)
    if haskey(options, :M)
        M = options[:M]
        if M isa Integer
            M = M_int2tuple(samplingstrategy, M, dict)
        end
    else
        M = oversamplingM(samplingstrategy, oversamplingfactor, dict)
    end
    L = tuple(map((api, Mi)->deduce_samplingparameter(OversamplingStyle(), samplingstrategy, api;
        oversamplingfactor=oversamplingfactor, M=Mi, verbose=verbose, options...), ApproximationProblems.unsafe_elements(ap), M)...)
    verbose && println("Sampling parameter: best match for M = $M is L = $L")
    L
end

M_int2tuple(::DefaultSamplingStrategy, M::Int, dict::TensorProductDict) =
    ntuple(k->round(Int,M*size(dict, k)/length(dict)))
oversamplingM(::DefaultSamplingStrategy, oversamplingfactor::Real, dict::TensorProductDict) =
    ntuple(k->round(Int,oversamplingfactor*size(dict, k)),Val(dimension(dict)))
oversamplingfactor_float2tuple(::DefaultSamplingStrategy, oversamplingfactor::Real, dict::TensorProductDict) =
    oversamplingfactor^(1/dimension(dict))

function deduce_oversampling_parameter(ss::OversamplingStyle, samplingstrategy::SamplingStrategy, args...; dict=dictionary(args...), verbose = false, oversamplingfactor = 2, options...)
    # In the absence of L, we deduce M and then find the best matching L
    # M is either supplied, or we compute it based on the (default) oversamplingfactor
    M = haskey(options, :M) ? options[:M] : round(Int, oversamplingfactor * length(dict))
    verbose && println("Sampling parameter: oversamplingfactor=$oversamplingfactor, length=$(length(dict)), M=$M")
    L = match_and_correct_sampling_parameter(samplingstrategy, args..., M; samplingstyle=ss, options...)
    verbose && println("Sampling parameter: best match for M = $M is L = $L")
    L
end

deduce_samplingparameter(ss::SamplingStyle, samplingstrategy::SamplingStrategy, platform::Platform, param; dict = dictionary(platform, param), options...) =
    deduce_samplingparameter(ss, samplingstrategy, dict; options...)
deduce_samplingparameter(ss::OversamplingStyle, samplingstrategy::SamplingStrategy, dict::Dictionary; options...) =
    deduce_oversampling_parameter(ss, samplingstrategy, dict; dict=dict, options...)
deduce_samplingparameter(ss::OversamplingStyle, samplingstrategy::SamplingStrategy, platform::Platform, param; dict=dictionary(platform, param), options...) =
    deduce_oversampling_parameter(ss, samplingstrategy, platform, param; dict=dict, options...)
deduce_samplingparameter(ss::InterpolationStyle, samplingstrategy::SamplingStrategy, dict::Dictionary; options...) =
    match_and_correct_sampling_parameter(samplingstrategy, dict, length(dict); samplingstyle=ss, options...)
# TODO: implement this one better (more general)
deduce_samplingparameter(::GramStyle, samplingstrategy::SamplingStrategy, dict::Dictionary; options...) = length(dict)

deduce_samplingparameter(::RectangularGramStyle, samplingstrategy::SamplingStrategy, dict::Dictionary; projectiondict, options...) = length(projectiondict)

deduce_samplingparameter(::GridStyle, samplingstrategy::SamplingStrategy, dict::Dictionary; options...) = param(dict)