text:
# corFun.jl

export CorParams, MaternParams,
    CorHypers, MaternHyper_ScInvChiSq,
    corr, corrMat, covMat, PDMat_adj;

abstract type CorParams end

mutable struct MaternParams <: CorParams
    ν::Float64
    lenscale::Float64
    MaternParams(ν, lenscale) = ν > 0.0 && lenscale > 0.0 ? new(ν, lenscale) : error("Both parameters must be positive.")
end
function MaternParams()
    MaternParams(2.5, 1.0)
end

abstract type CorHypers end

mutable struct MaternHyper_ScInvChiSq <: CorHypers
    lenscale_ν::Float64
    lenscale0::Float64
end
function MaternHyper_ScInvChiSq()
    MaternHyper_ScInvChiSq(5.0, 0.6)
end

function corr(d::T, params::MaternParams; logout::Bool=false) where T <: Real
    ## Matern

    ν = params.ν
    lenscale = params.lenscale

    d >= 0.0 || throw(error("Distance d must be non-negative."))
    lenscale > 0.0 || throw(error("Lengthscale parameter lenscale must be positive."))

    if ν == 0.5

        lcor = -d / lenscale

    elseif ν == 1.5

        s3ddl = sqrt(3.0) * d / lenscale
        lcor = log( s3ddl  + 1.0 ) - s3ddl

    elseif ν == 2.5

        s5ddl = sqrt(5.0) * d / lenscale
        lcor = log( s5ddl^2 / 3.0 + s5ddl + 1.0 ) - s5ddl

    elseif ν == Inf  # squared exponential correlation function

        lcor = -0.5 * ( d / lenscale )^2

    else
        throw(error("Unsupported value of ν."))
    end

    logout ? out = lcor : out = exp(lcor)
    return out
end

function corrMat(D::Matrix{T}, params::MaternParams; logout::Bool=false) where T <: Real
    if issymmetric(D)

        n = size(D,1)
        C = Matrix{T}(undef, n, n)

        for i = 1:n
            C[i,i] = corr(D[i,i], params, logout=logout)
            for j = (i+1):n
                C[i,j] = corr(D[i,j], params, logout=logout)
                C[j,i] = C[i,j]
            end
        end

    else

        n1, n2 = size(D)
        C = Matrix{T}(undef, n1, n2)

        for i = 1:n1
            for j = 1:n2
                C[i,j] = corr(D[i,j], params, logout=logout)
            end
        end

    end

    return Symmetric(C)
end

function covMat(D::Matrix{T}, variance::T, params::MaternParams; tol::Float64=1.0e-6) where T <: Real
    variance > 0.0 || throw("Variance must be positive. Its value is $(variance).")
    lv = log(variance)

    if issymmetric(D)

        n = size(D,1)
        C = Matrix{T}(undef, n, n)

        for i = 1:n
            C[i,i] = exp(lv + corr(D[i,i], params, logout=true))
            for j = (i+1):n
                C[i,j] =  exp(lv + corr(D[i,j], params, logout=true))
                C[j,i] = C[i,j]
            end
        end

        # C = PDMat_adj(C, tol) # do this manually where you need it.

    else

        n1, n2 = size(D)
        C = Matrix{T}(undef, n1, n2)

        for i = 1:n1
            for j = 1:n2
                C[i,j] =  exp(lv + corr(D[i,j], params, logout=true))
            end
        end

    end

    return C
end

function PDMat_adj(A::Union{Matrix{Float64}, Symmetric{Float64, Matrix{Float64}}, PDMat},
    maxadd::Float64=1.0e-6,
    epsfact::Float64=100.0, cumadd::Float64=0.0)

    try PDMat(A)
    catch excep
        if isa(excep, PosDefException) && cumadd <= maxadd
            a = epsfact * eps(Float64)
            A += a * I
            cumadd += a
            epsfactnext = 10.0 * epsfact
            PDMat_adj(A, maxadd, epsfactnext, cumadd)
        else
            PDMat(A) # just trigger original error
        end
    end
end


text:
struct ∇²SemWLS{Vt <: Union{AbstractArray, UniformScaling{Bool}},
        St <: AbstractArray}
    V::Vt
    s::St
end

function ∇²SemWLS(semwls::A) where {A <: SemWLS}
    return ∇²SemWLS(semwls.V, semwls.s)
end

function (diff::∇²SemWLS)(par, H::AbstractArray, model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    model.imply.imp_fun(model.imply.imp_cov, par)
    model.imply.gradient_fun(model.imply.∇Σ, par)
    J = (-2*(diff.s-model.imply.imp_cov)'*diff.V)'
    model.imply.hessian_fun(model.imply.∇²Σ, J, par)
    H .= 2*model.imply.∇Σ'*diff.V*model.imply.∇Σ + model.imply.∇²Σ
end

struct ∇²SemWLS_approx{Vt <: Union{AbstractArray, UniformScaling{Bool}},
        St <: AbstractArray}
    V::Vt
    s::St
end

function ∇²SemWLS_approx(semwls::A) where {A <: SemWLS}
    return ∇²SemWLS_approx(semwls.V, semwls.s)
end

function (diff::∇²SemWLS_approx)(par, H, model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    model.imply.imp_fun(model.imply.imp_cov, par)
    model.imply.gradient_fun(model.imply.∇Σ, par)
    H .= 2*model.imply.∇Σ'*diff.V*model.imply.∇Σ
end

ext: jl


text:
path = ["data"]

csvreading = CSVInput(name = "test_datastream.csv", descriptor = "1997");

xlsxreading = [
    XLSXInput(name = "test_datastream.xlsx", sheet = "1999", range = "A7:J14", descriptor = "1999"),
    XLSXInput(name = "test_datastream.xlsx", sheet = "2000", range = "A7:J14", descriptor = "2000")];

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

describing = Describe(col = :year);

ordering = Order(
    col = [:year, :from_code, :from_desc, :to_code, :to_desc, :value, :units],
    type = [Int, String, String, String, String, Int, String]);

renaming = [
    Rename(from = :IOCode, to = :from_code),
    Rename(from = :Name,   to = :from_desc)];

melting = Melt(
    on = [:from_desc],
    var = :to_desc,
    val = :value,
    type = Any);

adding = [
    Add(col = :region, val = "Colorado"),
    Add(col = :units, val = "millions of us dollars (USD)")];

mapping = Map(
        file = "regions.csv",
        from = [:from],
        to = [:to],
        input = [:region],
        output = [:region]);

replacing = [
    Replace(col = :value, from = "...",   to = 0),
    Replace(col = :value, from = missing, to = 0)];

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
struct Encoding{Data}
    data::Data
    title::String
end

struct Plot2D{Kind,X<:Encoding,Y<:Encoding,Color<:Union{Nothing,Encoding}}
    kind::Kind
    x::X
    y::Y
    color::Color
end

encoding(x) = Encoding(x.data, x.title)

plot2d(kind; x, y, color = nothing) =
    Plot2D(kind, encoding(x), encoding(y), color === nothing ? nothing : encoding(color))

struct LocationTimeKind end
const LocationTimePlot = Plot2D{LocationTimeKind}

function rescale(t::Period, unit::Type{<:Period})
    t, u = promote(t, unit(1))
    return t / u
end

function ascolumntable(table)
    if Tables.columnaccess(table)
        return table
    else
        return Tables.columntable(table)
    end
end

function MemoryAccesses.location_time(table; timeunit = Microsecond)::LocationTimePlot
    table = ascolumntable(table)

    unitstr = replace(string(timeunit(0)), r"^0 +" => "")
    xs = rescale(Nanosecond(1), timeunit) .* (table.time .- table.time[1])

    lb = minimum(table.access)
    ys = table.access .- lb

    zs = table.threadid

    return plot2d(
        LocationTimeKind();
        x = (data = xs, title = "Time [$unitstr]"),
        y = (data = ys, title = "Memory location"),
        color = (data = zs, title = "Thread ID"),
    )
end

function Base.show(io::IO, ::MIME"text/plain", plt::Plot2D)
    print(io, "Plot2D{", nameof(typeof(plt.kind)), ", …}")
    println(io)
    print(io, "  x: ", plt.x.title)
    println(io)
    print(io, "  y: ", plt.y.title)
    if plt.color !== nothing
        println(io)
        print(io, "  color: ", plt.color.title)
    end
    return
end


text:
function combine2tuples(vs::Vector...)
  combined::Vector = vec(collect(Base.product(vs...)))

  return combined
end

function combine2vectors(vs::T...) where T<:Vector
  combinedtuples::Vector = combine2tuples(vs...)

  N::Int = length(combinedtuples)
  M::Int = length(combinedtuples[1])

  #preallocate
  combined::Vector{T} = [T(undef, N) for i ∈ 1:M]

  for c ∈ 1:M, r ∈ 1:N
    combined[c][r] = combinedtuples[r][c]
  end

  return combined
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module LockUtils
export Lock, holding
type Lock
    sem::Base.Semaphore
    Lock() = new(Base.Semaphore(1))
end
function holding(f::Function, l::Lock, id)
    Base.acquire(l.sem)
    # print("$id HOLDING LOCK\n")
    f_ = f()
    Base.release(l.sem)
    # print("$id RELEASING LOCK\n")
    f_
end
function holding(f::Function, l::Lock)
    Base.acquire(l.sem)
    f_ = f()
    Base.release(l.sem)
    f_
end
export LevelTrigger, signal, waitSignal
type LevelTrigger
    condition::Lock
    val::Bool
    LevelTrigger() = new(Lock(), false)
end

function signal(l::LevelTrigger)
    holding(l.condition) do
        l.val = true
    end
end

function waitSignal(l::LevelTrigger)
    while (holding(l.condition) do
        l.val == false
        end)
        yield()

    end
end
function waitSignal(l::LevelTrigger, delay::Float64, totalDelay::Float64, timeout::Function)
    times = Int64(totalDelay/delay)
    for i =1:times
        if (holding(l.condition) do
            l.val
            end)
            return
        end
        wait(Timer(delay))
    end
    timeout()
    while !(holding(l.condition) do
            l.val
            end)
        wait(Timer(delay))
    end

end

end


text:
"""
Abstract type for a Robust MDP.

Should use the same interface as an MDP for things like actions, etc. If representative_mdp is implemented, then the interface methods for that MDP will be used when available.
"""
abstract RobustMDP{S,A}

POMDPs.isterminal(r::RobustMDP, s) = isterminal(representative_mdp(r), s)
POMDPs.discount(r::RobustMDP) = discount(representative_mdp(r))
POMDPs.create_action(r::RobustMDP) = create_action(representative_mdp(r))

next_action(gen::MCTS.RandomActionGenerator, rmdp::RobustMDP, s, ::RobustStateNode) = next_action(gen, rmdp, s)
function next_action(gen::MCTS.RandomActionGenerator, rmdp::RobustMDP, s)
    if isnull(gen.action_space)
        gen.action_space = Nullable{AbstractSpace}(actions(representative_mdp(rmdp)))
    end
    rand(gen.rng, actions(representative_mdp(rmdp), s, get(gen.action_space)))
end

type MDPCollection{S,A} <: RobustMDP{S,A}
    mdps::Vector{MDP{S,A}}
end

type MixedMDP{S,A} <: MDP{S,A}
    collection::MDPCollection{S,A}
    dist::Vector{Float64} # probability weights
end

function GenerativeModels.generate_sr{S,A}(mdp::MixedMDP{S,A}, s, a, rng::AbstractRNG)
    # select model
    r = rand(rng)
    psum = 0.0
    model = Nullable{MDP{S,A}}()
    for i in 1:length(dist)
        psum += mdp.dist[i]
        if r < psum
            model = Nullable{MDP{S,A}}(mdp.collection.mdps[i])
            break
        end
    end
    return generate_sr(get(model), s, a, rng)
end


text:
"""
      TrainTestSplits an Array\n
      --------------------\n
      x = [5,10,15,20]\n
      train, test = @tts x
       """
macro tts(x,at=.75)
    x = eval(x)
    TrainTestSplit(x,at)
end
"""
      Returns the normal distribution of an array.\n
      --------------------\n
      x = [5,10,15,20]\n
      norm = @norm x
       """
macro norm(x)
    x = eval(x)
    Normalizer(x).predict(x)
end
"""
      OneHotEncodes a dataframe\n
      Takes a symbol representing the column to one hot encode from and a DF.\n
      --------------------\n
      df = (:A => ["hello","world"], :B => ["Foo", "Bar"])\n
      encoded = @onehot df, :A
       """
macro onehot(df,symb)
    _onehotdf(df,symb)
end


text:
using FastTransforms
using Test

@testset "FastTransforms spin spherical harmonics" begin
    lmax = 4
    s = 1

    C1 = zeros(Complex{Float64}, lmax + 1, 2lmax + 1)
    C1[1, 2] = 1                 # l=1, m=-1

    C2 = zeros(Complex{Float64}, lmax + 1, 2lmax + 1)
    C2[2, 2] = 1                 # l=2, m=-1

    # Evalute
    P = plan_spinsph2fourier(C1, s)
    PS = plan_spinsph_synthesis(C1, s)
    F1 = PS * (P * C1)
    F2 = PS * (P * C2)

    # Calculate scalar product kernel
    F = conj(F1) .* F2

    # Integrate (by calculating l=0 mode)
    P = plan_sph2fourier(C1)
    PA = plan_sph_analysis(F)
    C = P \ (PA * F)

    # Should be zero
    @test abs(C[1, 1]) < 10eps()
end


text:
# This file was generated, do not modify it. # hide
Xc = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))
scitype(Xc)


text:
x = randn(100)
hist(x, nbins=10, norm="pdf", label="hist")
plot!(x->exp(-x^2/2)/√(2π), -3, 3, label="fit", smooth=true)
legend()


text:
# # Parallel Computing with Julia
#
# This course will cover:
#
# * Introduction to parallelism
#     * What is happening to our computers?
#
# * Parallelism strategies
#     * SIMD
#     * Multi-threading
#     * Tasks
#     * Multi-process
#         * Shared memory
#         * Distributed memory
#     * GPU programming
#
# * Challenges of parallel computing
#     * Order of execution
#         * execution of out order of Possibility
#         * simultaneous access and mutation
#     * Data access and movement
#     * Code access and movement
#     * Appropriately matching the parallelism strategy to your machine capabilities
#     * Appropriately matching the parallelism strategy with the problem at hand

#-

# ## What is happening to our computers!?
#
# ![](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/40-years-processor-trend.png)
#
# Not only have we gained multiple cores, but processors have become extremely
# complex, with multiple levels of caches, pipelines, predictions, speculations...
#
# ## What is hard about parallel computing
#   * We don't think in parallel
#   * We learn to write and reason about programs serially
#   * The desire for parallelism often comes _after_ you've written your algorithm (and found it too slow!)
#
# ## Summary:
#   * Current computer archetectures push us towards parallel programming for peak performance — even if we're not on a cluster!
#   * But it's hard to design good parallel algorithms
#   * And it's hard to express and reason about those algorithms

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule oxigraph_server_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("oxigraph_server")
JLLWrappers.@generate_main_file("oxigraph_server", UUID("ca7c0f55-ad95-5972-a3e2-8e5b4e3bed80"))
end  # module oxigraph_server_jll


text:
@testset "maximumlikelihood_gp.jl" begin
    n = 100

    x₁ = Variable("x₁", randn(n) / 3)
    x₂ = Variable("x₂", randn(n) / 10)

    ϕ = -.05 .+ x₁.value
    ξ = x₂.value

    σ = exp.(ϕ)
    θ = [-0.05; 1.0; 0.0; 1.0]

    pd = GeneralizedPareto.(σ, ξ)
    y = rand.(pd)

    @testset "gpfit(y; logscalecov, shapecov)" begin
        # model building with non-stationary logscale and shape
        fm = Extremes.gpfit(y,
            logscalecov = [x₁],
            shapecov = [x₂])

        # data is y
        @test fm.model.data.value ≈ y

        # logscale is x₁
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₁.value

        # shape is x₂
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₂.value

    end

    @testset "gpfit(y, initialvalues; logscalecov, shapecov)" begin
        # model building with non-stationary logscale and shape
        @test_logs Extremes.gpfit(y, zeros(Float64, 4),
            logscalecov = [x₁],
            shapecov = [x₂])

    end

    @testset "gpfit(df, datacol; logscalecovid, shapecovid)" begin
        # model building with non-stationary location, logscale and shape
        df = DataFrame(y = y, x1 = x₁.value, x2 = x₂.value)

        fm = Extremes.gpfit(df, :y, logscalecovid = [:x1], shapecovid = [:x2])

        # data is y
        @test fm.model.data.value ≈ y

        # logscale is x₁
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₁.value

        # shape is x₂
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₂.value

    end

    @testset "gpfit(df, datacol, initialvalues; logscalecovid, shapecovid)" begin
        # model building with non-stationary location, logscale and shape
        df = DataFrame(y = y, x1 = x₁.value, x2 = x₂.value)

        @test_logs Extremes.gpfit(df, :y, zeros(Float64, 4), logscalecovid = [:x1], shapecovid = [:x2])

    end

    @testset "gpfit(model, initialvalues)" begin
        # non-stationary location, logscale and shape
        model = Extremes.ThresholdExceedance(Variable("y", y),
            logscalecov = [x₁],
            shapecov = [x₂])

        fm = Extremes.gpfit(model, zeros(Float64, 4))

        # data is y
        @test fm.model.data.value ≈ y

        # logscale is x₁
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₁.value

        # shape is x₂
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₂.value

    end

end


text:
# -*- mode: julia; -*-
module WAVPlay
wavplay(data, fs) = warn("wavplay is not currently implemented on $(Sys.KERNEL)")
end # module


text:
function SeisBPFilter(d,dt,f1,f2,f3,f4)
    Ndims = ndims(d);
    N     = size(d);
    nt    = N[1];
    x     = reshape(d,nt,prod(N[2:end]));
    nf=4*nextpow(2,nt)

    i1 = floor(Int,nf*f1*dt)+1;
    i2 = floor(Int,nf*f2*dt)+1;
    i3 = floor(Int,nf*f3*dt)+1;
    i4 = floor(Int,nf*f4*dt)+1;
    up=collect(1:1:(i2-i1))/(i2-i1);
    down=collect((i4-i3):-1:1)/(i4-i3);
    aux1 = append!(vec(zeros(1,i1)),up,vec(ones(1,i3-i2)), down,vec(zeros(1,Int(nf/2)+1-i4)));
    aux2 = reverse(aux1[2:Int(nf/2)]);
    c = 0;
    F=append!(aux1,aux2);

    Phase = (pi/180.)*append!([0],-c*ones(1,Int(nf/2-1)),[0],c*ones(1,Int(nf/2-1)));
    Transfer = F.*exp.(-im .*Phase);

    xp=zeros(nf,size(x)[2]);
    xp[1:size(x)[1],:]=x;
    X = fft(xp,1);
    Y=zeros(ComplexF64,size(X));
    for k = 1:prod(N[2:end])
         Y[:,k] = Transfer .*X[:,k];
    end

    o = ifft(Y,1);

    o = real(o[1:nt,:]);

    o = reshape(o, (nt, N[end]));
    return o
end


text:
# ---
# title: 567. Permutation in String
# id: problem567
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/permutation-in-string/description/>
# hidden: true
# ---
#
# Given two strings **s1** and **s2** , write a function to return true if
# **s2** contains the permutation of **s1**. In other words, one of the first
# string's permutations is the **substring** of the second string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s1 = "ab" s2 = "eidbaooo"
#     Output: True
#     Explanation: s2 contains one permutation of s1 ("ba").
#
#
# **Example 2:**
#
#
#
#     Input: s1= "ab" s2 = "eidboaoo"
#     Output: False
#
#
#
#
# **Constraints:**
#
#   * The input strings only contain lower case letters.
#   * The length of both given strings is in range [1, 10,000].
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Flux
using Flux: mse, back!
using Flux.Optimise: _update_params!
using Statistics
using Formatting

"""
HillClimber

    Agent attempts to improve by random parameter search.

    Run `N` epsiodes using a random parameterization of `Q` and record the average score. Next, make a random update to the parameterization, and run another `N` episodes to calculate the fitness of the new `Q`. If the new `Q` is fitter than the old `Q`, switch to the new `Q`. Otherwise, generate a new random update and repeat.
```
"""
mutable struct HillClimber
    W
    σ
end

HillClimber(A, S) = HillClimber(randn(A, S), 0.1)

function action(agent::HillClimber, env, ϵ = 0.0)
    if rand() < ϵ
        return rand(env.action_space)  # explore
    else
        return argmax(agent.W * env.state)  # exploit
    end
end

function perturb!(agent)
    W_ = deepcopy(agent.W)
    agent.W .+= agent.σ * randn(size(agent.W))
    return W_
end

function run_episode!(agent::HillClimber, env; ϵ = 0.05)
    total_reward = 0.0
    steps = 0
    s′, done = reset!(env)
    while !done
        # render(env)
        a = action(agent, env, ϵ)
        s = deepcopy(s′)
        s′, r, done, info = step!(env, a)
        # println("s = $s, a = $a, r = $r, s′ = $s′, done = $done")
        total_reward += r
        steps += 1
        update!(agent)
    end
    return total_reward, steps
end

function evaluate(agent::HillClimber, env; n = 100)
    history = []
    for i in 1:n
        total_reward = 0.0
        steps = 0
        s′, done = reset!(env)
        while !done
            a = action(agent, env)  # no exploration here (ϵ = 0.0)
            s = deepcopy(s′)
            s′, r, done, info = step!(env, a)
            total_reward += r
            steps += 1
            # println("s = $s, a = $a, r = $r, s′ = $s′, done = $done")
        end
        # @info "Finished episode $i" total_reward steps
        push!(history, (total_reward=total_reward, steps=steps))
    end
    return mean([e.total_reward for e in history])
end

function train!(agent::HillClimber, env; max_episodes = Inf)
    init_memory!(agent, env)
    @info "Training agent..."
    ϵ = ϵ0
    history = []
    iter = 0
    improving = true
    fitness = evaluate(agent, env)
    while improving && episodes < max_episodes
        iter += 1
        W = perturb!(agent)
        fitness = evaluate(agent, env)
        fitness < agent.fitness && revert!(agent, W)
        @info format("iter: {}, score: {}, epsilon: {}", episodes, fitness)
    end
    return history
end


text:
# Sample script for plotting fieldlines with handpicked seeds, multi-processing version.
#
# To run on a single node,
# julia -p $ncores demo_fieldline_mp_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed, ParallelDataTransfer, Glob
@everywhere using Vlasiator, PyPlot, PyCall, Printf, LaTeXStrings, FieldTracer
@everywhere using Vlasiator: RE

function generate_seeds(coordmin, coordmax, dim_, nseeds)
   seeds = Matrix{Float64}(undef, 2, nseeds)
   for i in 1:nseeds
      seeds[1,i] = coordmin[dim_[1]] +
         (coordmax[dim_[1]] - coordmin[dim_[1]]) / nseeds * (i - 1)
      seeds[2,i] = -20RE
   end
   seeds
end

@everywhere function init_figure(pArgs, norm, ticks, seeds, extent)
   fig, ax = plt.subplots(1, 1; num=myid(),
      figsize=(6, 8), constrained_layout=true)

   fontsize = "x-large"

   ax.set_aspect("equal")
   # Set border line widths
   for loc in ("left", "bottom", "right", "top")
      edge = get(ax.spines, loc, nothing)
      edge.set_linewidth(2.0)
   end
   ax.xaxis.set_tick_params(width=2.0, length=3)
   ax.yaxis.set_tick_params(width=2.0, length=3)
   ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())

   ax.set_xlabel(pArgs.strx; fontsize)
   ax.set_ylabel(pArgs.stry; fontsize)

   ax.set_title("Density"; fontsize)

   x1, x2 = Vlasiator.get_axis(pArgs)

   range1 = searchsortedfirst(x1, extent[1]):searchsortedlast(x1, extent[2])
   range2 = searchsortedfirst(x2, extent[3]):searchsortedlast(x2, extent[4])

   fakedata = zeros(Float32, length(range2), length(range1))
   c = ax.pcolormesh(x1[range1], x2[range2], fakedata; norm, cmap=matplotlib.cm.turbo)

   format = matplotlib.ticker.FormatStrFormatter("%.1f")
   cb1 = colorbar(c; ax, ticks, format)
   cb1.ax.set_ylabel("[amu/cc]"; fontsize)
   cb1.outline.set_linewidth(1.0)

   fakeline = [0.0, 1.0]
   ls = [ax.plot(fakeline, fakeline,  color="w") for _ in 1:size(seeds,2)]

   return fig, ax, c, ls, range1, range2
end

@everywhere function update_plot!(ax, c, ls, range1, range2, dim_, seeds, grid1, grid2,
   outdir, file)
   isfile(outdir*file[end-8:end-5]*".png") && return

   println("file = $file")
   meta = load(file)

   data = Vlasiator.prep2d(meta, "proton/vg_rho", :mag)'
   c.set_array(data[range2,range1] ./ 1f6)

   str_title = @sprintf "Density pulse run, t= %4.1fs" meta.time
   ax.set_title(str_title; fontsize="x-large")

   b = meta["vg_b_vol"]
   b1 = reshape(b[dim_[1],:], meta.ncells[dim_[1]], meta.ncells[dim_[2]])
   b2 = reshape(b[dim_[2],:], meta.ncells[dim_[1]], meta.ncells[dim_[2]])
   # Find existing arrow annotations
   annotations = [child for child in ax.get_children() if
      pybuiltin(:isinstance)(child, matplotlib.text.Annotation)]
   # Remove existing arrows
   for a in annotations
      a.remove()
   end
   # Add new arrows along field lines
   for i in axes(seeds,2)
      startx, starty = seeds[:,i]
      x1, y1 = trace(b1, b2, startx, starty, grid1, grid2;
         ds=0.5, maxstep=4000, gridtype="ndgrid")
      x1 ./= RE
      y1 ./= RE
      if length(x1) < 5; continue; end
      ls[i][1].set_xdata(x1)
      ls[i][1].set_ydata(y1)
      add_arrow(ls[i][1])
   end

   savefig(outdir*file[end-8:end-5]*".png", bbox_inches="tight")
end

function make_jobs(files)
   for f in files
       put!(jobs, f)
   end
end

@everywhere function do_work(jobs, status,
   outdir, pArgs, norm, ticks, grid1, grid2, dim_, seeds, extent)

   fig, ax, c, ls, range1, range2 = init_figure(pArgs, norm, ticks, seeds, extent)

   while true
      file = take!(jobs)
      update_plot!(ax, c, ls, range1, range2, dim_, seeds, grid1, grid2, outdir, file)
      put!(status, true)
   end
   close(fig)
end

############################################################################################
files = glob("bulk*.vlsv", ".")

nfile = length(files)
# Set output directory
outdir = "out/"

const jobs   = RemoteChannel(()->Channel{String}(nfile))
const status = RemoteChannel(()->Channel{Bool}(nworkers()))

axisunit = EARTH # contour plot axes unit
extent = [0., 20., -20., 20.] # [RE], default full domain: [-Inf32, Inf32, -Inf32, Inf32]

# Upper/lower limits for each variable
ρmin, ρmax = 0.0, 11.0      # [amu/cc]

meta = load(files[1])
# Construct pieces for plotting
pArgs = Vlasiator.set_args(meta, "proton/vg_rho", axisunit; normal=:none)
norm, ticks = Vlasiator.set_colorbar(Linear, ρmin, ρmax)

# Mark spatial dimensions
dim_ = pArgs.stry[1] == 'Z' ? (1,3) : (1,2)

(;coordmin, coordmax, ncells) = meta

# Generate regular Cartesian range
grid1 = range(coordmin[dim_[1]], coordmax[dim_[1]], length=ncells[dim_[1]])
grid2 = range(coordmin[dim_[2]], coordmax[dim_[2]], length=ncells[dim_[2]])
# Generate seeds for in-plane field line tracing
nseeds = 10
seeds = generate_seeds(coordmin, coordmax, dim_, nseeds)

println("Total number of files: $nfile")
println("Running with $(nworkers()) workers...")

@async make_jobs(files) # Feed the jobs channel with all files to process.

@sync for p in workers()
   @async remote_do(do_work, p, jobs, status,
      outdir, pArgs, norm, ticks, grid1, grid2, dim_, seeds, extent)
end

let n = nfile
   t = @elapsed while n > 0 # wait for all jobs to complete
      take!(status)
      n -= 1
   end
   println("Finished in $(round(t, digits=2))s.")
end


text:
module OpenMPIData

using Reexport

include("download.jl")

end # module


text:
"""
*MultipleTesting* package

* Adjusting p-values for multiple testing
* Estimating the fraction π₀ of tests under the null hypothesis
* Combination of p-values
* Higher criticism
"""
module MultipleTesting

using StatsBase
import StatsBase: fit

using Distributions
import Distributions: estimate

import SpecialFunctions: digamma

export PValues,
    adjust,
    PValueAdjustment,
    Bonferroni,
    BenjaminiHochberg,
    BenjaminiHochbergAdaptive,
    BenjaminiYekutieli,
    BenjaminiLiu,
    Hochberg,
    Holm,
    Hommel,
    Sidak,
    ForwardStop,
    BarberCandes,
    estimate,
    Pi0Estimator,
    Storey,
    StoreyBootstrap,
    LeastSlope,
    Oracle,
    TwoStep,
    RightBoundary,
    CensoredBUM,
    CensoredBUMFit,
    BUM,
    BUMFit,
    FlatGrenander,
    ConvexDecreasing,
    ConvexDecreasingFit,
    isin,
    fit,
    BetaUniformMixtureModel,
    PValueCombination,
    combine,
    Fisher,
    Logit,
    Stouffer,
    Tippett,
    Simes,
    Wilkinson,
    Minimum,
    HigherCriticismScores,
    HigherCriticismThreshold

include("types.jl")
include("utils.jl")
include("pval-adjustment.jl")
include("pi0-estimators.jl")
include("model.jl")
include("combinations.jl")
include("higher-criticism.jl")

end


text:
#=

    This is a port of the SPPL model poisson_mean_inference.pynb

    Here is the PDF of Y = 0..20 (exact probabilities) from the SPPL model:
    (0, 0.19572329209779615),
    (1, 0.24541518485806976),
    (2, 0.15633363576277434),
    (3, 0.07275764007128807),
    (4, 0.03749784880330323),
    (5, 0.03225414602603766),
    (6, 0.036272539178041456),
    (7, 0.03919655164760693),
    (8, 0.03782130413444809),
    (9, 0.032547494609049904),
    (10, 0.025221805385485633),
    (11, 0.017769927037907036),
    (12, 0.011477227802063514),
    (13, 0.006844292955361095),
    (14, 0.00379364785941744),
    (15, 0.0019704640279916506),
    (16, 0.000975431521715664),
    (17, 0.0004844975064651556),
    (18, 0.0002799391781190319),
    (19, 0.00023526933546341676)
    (20, 0.0002864142016352986)

    This Turing.jl model:

    Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            X    0.0256    1.0030     0.0100    0.0124   6401.6166    0.9999     2199.1125
            Y    4.2629    7.6756     0.0768    0.0975   6534.2074    1.0000     2244.6607

    Quantiles
    parameters      2.5%     25.0%     50.0%     75.0%     97.5%
        Symbol   Float64   Float64   Float64   Float64   Float64

            X   -1.9514   -0.6385    0.0195    0.7112    2.0227
            Y    0.0000    1.0000    2.0000    4.0000   26.0000

    Distributions of variable Y:
    1.00000 =>    2448  (0.244800)
    0.00000 =>    1827  (0.182700)
    2.00000 =>    1776  (0.177600)
    3.00000 =>     987  (0.098700)
    4.00000 =>     614  (0.061400)
    5.00000 =>     382  (0.038200)
    6.00000 =>     287  (0.028700)
    7.00000 =>     237  (0.023700)
    8.00000 =>     192  (0.019200)
    9.00000 =>     145  (0.014500)
    11.00000 =>     113  (0.011300)
    10.00000 =>     106  (0.010600)
    13.00000 =>      94  (0.009400)
    12.00000 =>      80  (0.008000)
    14.00000 =>      75  (0.007500)
    15.00000 =>      59  (0.005900)
    16.00000 =>      55  (0.005500)
    20.00000 =>      55  (0.005500)
    17.00000 =>      42  (0.004200)
    18.00000 =>      36  (0.003600)
    24.00000 =>      30  (0.003000)
    21.00000 =>      27  (0.002700)
    23.00000 =>      25  (0.002500)
    19.00000 =>      23  (0.002300)
    26.00000 =>      21  (0.002100)
    27.00000 =>      19  (0.001900)
    25.00000 =>      15  (0.001500)
    22.00000 =>      15  (0.001500)
    30.00000 =>      13  (0.001300)
    34.00000 =>      13  (0.001300)
    36.00000 =>      12  (0.001200)
    32.00000 =>      11  (0.001100)
    28.00000 =>      11  (0.001100)
    29.00000 =>      11  (0.001100)
    46.00000 =>      10  (0.001000)
    52.00000 =>       9  (0.000900)
    31.00000 =>       9  (0.000900)
    35.00000 =>       8  (0.000800)
    39.00000 =>       8  (0.000800)
    33.00000 =>       8  (0.000800)
    40.00000 =>       8  (0.000800)
    38.00000 =>       7  (0.000700)
    37.00000 =>       7  (0.000700)
    42.00000 =>       6  (0.000600)
    57.00000 =>       5  (0.000500)
    64.00000 =>       4  (0.000400)
    49.00000 =>       4  (0.000400)
    47.00000 =>       4  (0.000400)
    60.00000 =>       3  (0.000300)
    45.00000 =>       3  (0.000300)
    41.00000 =>       3  (0.000300)
    61.00000 =>       3  (0.000300)
    53.00000 =>       3  (0.000300)
    48.00000 =>       3  (0.000300)
    56.00000 =>       2  (0.000200)
    55.00000 =>       2  (0.000200)
    110.00000 =>       2  (0.000200)
    67.00000 =>       2  (0.000200)
    73.00000 =>       2  (0.000200)
    43.00000 =>       2  (0.000200)
    50.00000 =>       2  (0.000200)
    51.00000 =>       2  (0.000200)
    139.00000 =>       1  (0.000100)
    63.00000 =>       1  (0.000100)
    54.00000 =>       1  (0.000100)
    58.00000 =>       1  (0.000100)
    69.00000 =>       1  (0.000100)
    95.00000 =>       1  (0.000100)
    66.00000 =>       1  (0.000100)
    93.00000 =>       1  (0.000100)
    44.00000 =>       1  (0.000100)
    70.00000 =>       1  (0.000100)
    80.00000 =>       1  (0.000100)
    65.00000 =>       1  (0.000100)
    97.00000 =>       1  (0.000100)

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function poisson_mean_inference(n=20)
    X ~ truncated(Normal(),-4, 4)
    Y ~ Poisson(1+2*abs(X^3))
    return Y
end

n=20
model = poisson_mean_inference(n)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:Y)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export marginalrule

@marginalrule GCV(:y_x) (m_y::Any, m_x::Any, q_z::Any, q_κ::Any, q_ω::Any) = begin

    ksi = mean(q_κ) ^ 2 * var(q_z) + mean(q_z) ^ 2 * var(q_κ) + var(q_z) * var(q_κ)
    A = exp(-mean(q_ω) + var(q_ω) / 2)
    B = exp(-mean(q_κ) * mean(q_z) + ksi / 2)
    W = [ precision(m_y) + A * B -A * B; -A * B precision(m_x) + A * B ]
    m = cholinv(W) * [ mean(m_y) * precision(m_y); mean(m_x) * precision(m_x) ]

    return MvNormalMeanPrecision(m, W)
end


text:
"""
    unsafe_decompress!([ctx = DecompressionContext(;kwargs...)], dest::Ptr{T}, dest_size, src::Ptr{UInt8}, src_size; kwargs...)

Decompress data from `src` buffer, into `dest` buffer using context `ctx` with no checks, return count of elements written to dest or negative value if an error.
`dest_size` and `src_size` are the byte sizes of buffers `dest` and `src` respectively

Return count of the elements written to dest

The `unsafe` prefix on this function indicates that no validation is performed on the pointers dest and src to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_decompress!(ctx::DecompressionContext, dest::Ptr{T}, dest_size, src::Ptr{UInt8}, src_size) where {T}
    res = Lib.blosc2_decompress_ctx(ctx.context, src, src_size, dest, dest_size)
    return res >= 0 ?
    ceil(Int64, res / sizeof(T)) :
    res
end

unsafe_decompress!(dest::Ptr{T}, dest_size, src::Ptr{UInt8}, src_size; kwargs...) where {T} =
    unsafe_decompress!(DecompressionContext(;kwargs...), dest, dest_size, src, src_size)

"""
    decompress!([ctx = DecompressionContext(;kwargs...)], dest::Vector{T}, src::Vector{UInt8}, dest_offset = 1, src_offset = 1; kwargs...)

Decompress data from buffer `src` at offset `src_offset` (1-indexed) into vector `dest` at offset `dest_offset` (1-indexed).
Return count of the elements written to `dest`
"""
function decompress!(ctx::DecompressionContext, dest::Vector{T},  src::Vector{UInt8}, dest_offset = 1, src_offset = 1) where {T}
    @boundscheck begin
        ulen = uncompressed_length(T, src, src_offset)
        checkbounds(dest, dest_offset:dest_offset + ulen - 1)
    end
    src_size = length(src) - src_offset + 1
    dest_size = (length(dest) - dest_offset + 1) * sizeof(T)
    result = GC.@preserve dest src begin
        unsafe_decompress!(ctx, pointer(dest, dest_offset), dest_size, pointer(src, src_offset), src_size)
    end
    result < 0 && error("Decompression error. Error code: $result")
    return  result
end

decompress!(dest::Vector{T}, src::Vector{UInt8}, dest_offset = 1, src_offset = 1;kwargs...) where {T} =
        decompress!(DecompressionContext(;kwargs...), dest, src, dest_offset, src_offset)

"""
    decompress([ctx = DecompressionContext(kwargs...)], ::Type{T}, src::Array{UInt8}, [src_offset = 1]; kwargs...)::Vector{T}

Return `Vector{T}` consisting of decompression data from `src` starting at `src_offset` using context `ctx`
"""
function decompress(ctx::DecompressionContext, ::Type{T}, src::Array{UInt8}, src_offset = 1) where {T}
    result = Vector{T}(undef, uncompressed_length(T, src))
    @inbounds r = decompress!(ctx, result, src, 1, src_offset)
    r < 0 && error("Decompression error. Error code: $r")
    return result
end
decompress(::Type{T}, src::Array{UInt8}, src_offset = 1; kwargs...) where {T} = decompress(DecompressionContext(;kwargs...), T, src, src_offset)

"""
    decompress_items!(ctx::DecompressionContext, dest::Vector{T},  src::Vector{UInt8}, range::UnitRange{<:Integer}, dest_offset, src_offset)

Decompress items in `range` from buffer `src` at offset `src_offset` (1-indexed) into vector `dest` at offset `dest_offset` (1-indexed).
Return count of the elements written to `dest`
"""
function decompress_items!(ctx::DecompressionContext, dest::Vector{T},
            src::Vector{UInt8}, range::UnitRange{<:Integer}, dest_offset = 1, src_offset = 1) where {T}

    @boundscheck begin
        checkbounds(src, src_offset)
        checkbounds(dest, dest_offset:dest_offset + length(range) - 1)
    end
    src_size = length(src) - src_offset + 1
    dest_size = (length(dest) - dest_offset + 1) * sizeof(T)

    result = GC.@preserve dest src begin
        Lib.blosc2_getitem_ctx(ctx.context,
            pointer(src, src_offset), src_size,
            first(range) - 1, length(range),
            pointer(dest, dest_offset), dest_size
        )
    end
    result < 0 && error("Decompression error. Error code: $result")
    return ceil(Int, result / sizeof(T))
end

decompress_items!(dest::Vector{T}, src::Vector{UInt8}, range::UnitRange{<:Integer}, dest_offset = 1, src_offset = 1; kwargs...) where {T} =
            decompress_items!(DecompressionContext(;kwargs...), dest, src, range, dest_offset, src_offset)


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

symbol() = :testing
emit(symbol)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.symbol)}, svec(), symbol() in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/lower/symbol.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.symbol)}"(%arg0: !jlir<"typeof(Main.symbol)">) -> !jlir.Symbol attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<":testing">} : () -> !jlir.Symbol
# CHECK-NEXT:     "jlir.return"(%0) : (!jlir.Symbol) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.symbol)}"(%arg0: !jlir<"typeof(Main.symbol)">) -> !jlir.Symbol attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<":testing">} : () -> !jlir.Symbol
# CHECK-NEXT:     return %0 : !jlir.Symbol
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.symbol)}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>) -> !llvm.ptr<struct<"struct_jl_value_type", opaque>> attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %1 = llvm.inttoptr %0 : i64 to !llvm.ptr<struct<"struct_jl_value_type", opaque>>
# CHECK-NEXT:     llvm.return %1 : !llvm.ptr<struct<"struct_jl_value_type", opaque>>
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.symbol)}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>) -> !llvm.ptr<struct<"struct_jl_value_type", opaque>> attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.symbol)}"(%arg0) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>) -> !llvm.ptr<struct<"struct_jl_value_type", opaque>>
# CHECK-NEXT:     llvm.return %0 : !llvm.ptr<struct<"struct_jl_value_type", opaque>>
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
#=
 Methods and types associated with single and double layer potentials
=#

module ImmersedLayers

using Reexport
using DocStringExtensions
@reexport using CartesianGrids
@reexport using RigidBodyTools
@reexport using ConstrainedSystems

import CartesianGrids: Laplacian

#using CartesianGrids
#using RigidBodyTools

using LinearAlgebra

using UnPack

export DoubleLayer, SingleLayer, MaskType, Mask, ComplementaryMask,
        DoubleLayer!, SingleLayer!, Mask!,
        SurfaceCache,SurfaceScalarCache,SurfaceVectorCache,
        @ilmproblem,construct_system,update_system,update_system!,
        BasicScalarILMProblem,BasicVectorILMProblem,prob_cache,
        AbstractScalingType,GridScaling,IndexScaling,
        BasicILMCache,ILMSystem,
        ODEFunctionList,zeros_sol,init_sol,
        similar_grid,similar_gridgrad,similar_gridcurl,similar_griddiv,similar_gridgradcurl,
        similar_surface,similar_surfacescalar,
        zeros_grid,zeros_gridgrad,zeros_gridcurl,zeros_griddiv,zeros_gridgradcurl,
        zeros_surface,zeros_surfacescalar,
        ones_grid,ones_gridgrad,ones_gridcurl,ones_griddiv,ones_gridgradcurl,
        ones_surface,ones_surfacescalar,
        x_grid,y_grid,x_gridcurl,y_gridcurl,x_griddiv,y_griddiv,x_gridgrad,y_gridgrad,
        AbstractExtraILMCache,AbstractScalarILMProblem,AbstractVectorILMProblem,
        ConvectiveDerivativeCache,convective_derivative, convective_derivative!,
        regularize!, interpolate!,
        regularize_normal!,normal_interpolate!,
        regularize_normal_symm!,normal_interpolate_symm!,
        regularize_normal_cross!,normal_cross_interpolate!,
        regularize_normal_dot!,normal_dot_interpolate!,
        surface_curl!,surface_divergence!,surface_grad!,inverse_laplacian!,
        laplacian!,Laplacian,
        surface_curl_cross!,surface_divergence_cross!,surface_grad_cross!,
        surface_divergence_symm!,surface_grad_symm!,
        mask,mask!,complementary_mask,complementary_mask!,
        ScalarPotentialCache,VectorPotentialCache,VectorFieldCache,
        vectorpotential_from_masked_curlv!,scalarpotential_from_masked_divv!,
        vectorpotential_from_curlv!,vecfield_from_vectorpotential!,masked_curlv_from_curlv_masked!,
        scalarpotential_from_divv!,masked_divv_from_divv_masked!,vecfield_from_scalarpotential!,
        divv_masked_from_masked_divv!,curlv_masked_from_masked_curlv!,
        vecfield_helmholtz!,vectorpotential_uniformvecfield!,scalarpotential_uniformvecfield!,vecfield_uniformvecfield!,
        create_CLinvCT,create_CL2invCT,create_GLinvD,create_nRTRn,create_RTLinvR,
        create_GLinvD_cross,create_CLinvCT_scalar,create_surface_filter,
        AreaRegionCache,LineRegionCache,PointRegionCache,arccoord,
        ForcingModelAndRegion,apply_forcing!,
        solve, @snapshotoutput,@scalarsurfacemetric,@vectorsurfacemetric,
        prescribed_surface_jump!,prescribed_surface_average!

abstract type AbstractBasicCache{N,GCT} end

abstract type AbstractScalingType end

abstract type GridScaling <: AbstractScalingType end
abstract type IndexScaling <: AbstractScalingType end

abstract type LayerType{N} end

const DEFAULT_SCALING = GridScaling
const DEFAULT_DATA_TYPE = Float64
const DEFAULT_DDF = CartesianGrids.Yang3

"""
$(TYPEDEF)

A system of operators and caches for immersed layer problems. This is constructed
by [`construct_system`](@ref)
"""
mutable struct ILMSystem{static,PT,N,PHT,BCF,FF,DTF,MTF,BCT,ECT}

  phys_params :: PHT
  bc :: BCF
  forcing :: FF
  timestep_func :: DTF
  motions :: MTF
  base_cache :: BCT
  extra_cache :: ECT

end

include("cartesian_extensions.jl")
include("tools.jl")
include("cache.jl")
include("problem.jl")
include("system.jl")
include("layers.jl")
include("grid_operators.jl")
include("surface_operators.jl")
include("forcing.jl")
include("matrix_operators.jl")
include("helmholtz.jl")
include("timemarching.jl")
include("output.jl")
include("bc.jl")

include("plot_recipes.jl")

@deprecate DoubleLayer surface_divergence!
@deprecate SingleLayer regularize!

end


text:
# These functions are for testing purposes only
# They need to be tested and abstractized

# remember that the psoils are corrected by surface tension, psoil = p_25 * f_st
#=
function test_soil_from_psoil(
            node::SPACMono{FT},
            psoils::Array{FT,1},
            Δt::FT=FT(10)
) where {FT<:AbstractFloat}
    # unpack the values
    @unpack n_canopy = node;
    n_sl = node.canopy_rt.nAzi * node.canopy_rt.nIncl;

    # update the soil water contents and potential in each layer
    for i_root in 1:node.n_root
        node.swc[i_root] = max(node.mswc[i_root], soil_swc(node.plant_hs.roots[i_root].sh, psoils[i_root]));
        node.plant_hs.roots[i_root].p_ups = psoils[i_root];
    end

    e_sum = FT(0);

    # loop through the canopy layers
    for i_can in 1:node.n_canopy
        rt_layer = node.n_canopy + 1 - i_can;
        f_view = (node.can_opt.Ps[rt_layer] + node.can_opt.Ps[rt_layer+1]) / 2;

        e_i_can = FT(0);
        for i_leaf in 1:(n_sl+1)
            e_i_can += node.plant_ps[i_can].g_lw[i_leaf] *
                       (node.plant_ps[i_can].p_sat - node.envirs[i_can].p_H₂O) /
                       node.envirs[i_can].p_atm *
                       node.plant_ps[i_can].LAIx[i_leaf] *
                       node.plant_ps[i_can].LA;
        end

        # calculate the photosynthetic rates
        gas_exchange!(node.photo_set, node.plant_ps[i_can], node.envirs[i_can], GswDrive());
        gsw_control!(node.photo_set, node.plant_ps[i_can], node.envirs[i_can]);

        # use the ball-berry model here for now as the ∂A/∂E and ∂A/∂Θ functions are not yet ready
        gsw_ss = stomatal_conductance(node.stomata_model, node.plant_ps[i_can], node.envirs[i_can], FT(1));

        # assume τ = 10 minutes
        for i_leaf in 1:(n_sl+1)
            node.plant_ps[i_can].g_sw[i_leaf] += (gsw_ss[i_leaf] - node.plant_ps[i_can].g_sw[i_leaf]) / 600 * Δt;
        end

        e_sum += e_i_can;
    end

    # update root flow rates
    roots_flow!(node.plant_hs.roots,
                node.plant_hs.container_k,
                node.plant_hs.container_p,
                node.plant_hs.container_q,
                e_sum);

    return nothing
end

function test_soil_from_swc(
            node::SPACMono{FT},
            swcs::Array{FT,1},
            Δt::FT=FT(10)
) where {FT<:AbstractFloat}
    # unpack the values
    @unpack n_canopy = node;
    n_sl = node.canopy_rt.nAzi * node.canopy_rt.nIncl;

    # update the soil water contents and potential in each layer
    for i_root in 1:node.n_root
        node.swc[i_root] = max(node.mswc[i_root], swcs[i_root]);
        node.p_soil[i_root] = soil_p_25_swc(node.plant_hs.roots[i_root].sh,
                                            node.swc[i_root]) *
                              node.plant_hs.roots[i_root].f_st;
        # pass the pressure to root_hs
        node.plant_hs.roots[i_root].p_ups = node.p_soil[i_root];
    end

    e_sum = FT(0);

    # loop through the canopy layers
    for i_can in 1:node.n_canopy
        rt_layer = node.n_canopy + 1 - i_can;
        f_view = (node.can_opt.Ps[rt_layer] + node.can_opt.Ps[rt_layer+1]) / 2;

        e_i_can = FT(0);
        for i_leaf in 1:(n_sl+1)
            e_i_can += node.plant_ps[i_can].g_lw[i_leaf] *
                       (node.plant_ps[i_can].p_sat - node.envirs[i_can].p_H₂O) /
                       node.envirs[i_can].p_atm *
                       node.plant_ps[i_can].LAIx[i_leaf] *
                       node.plant_ps[i_can].LA;
        end

        # calculate the photosynthetic rates
        gas_exchange!(node.photo_set, node.plant_ps[i_can], node.envirs[i_can], GswDrive());
        gsw_control!(node.photo_set, node.plant_ps[i_can], node.envirs[i_can]);

        # use the ball-berry model here for now as the ∂A/∂E and ∂A/∂Θ functions are not yet ready
        gsw_ss = stomatal_conductance(node.stomata_model, node.plant_ps[i_can], node.envirs[i_can], FT(1));

        # assume τ = 10 minutes
        for i_leaf in 1:(n_sl+1)
            node.plant_ps[i_can].g_sw[i_leaf] += (gsw_ss[i_leaf] - node.plant_ps[i_can].g_sw[i_leaf]) / 600 * Δt;
        end

        e_sum += e_i_can;
    end

    # update root flow rates
    roots_flow!(node.plant_hs.roots,
                node.plant_hs.container_k,
                node.plant_hs.container_p,
                node.plant_hs.container_q,
                e_sum);

    return nothing
end
=#

function test_soil_from_psoil!(
            node::SPACMono{FT},
            psoils::Array{FT,1},
            e_sum::FT
) where {FT<:AbstractFloat}
    # update the soil water contents and potential in each layer
    for i_root in 1:node.n_root
        node.swc[i_root] = max(node.mswc[i_root],
                               soil_swc(node.plant_hs.roots[i_root].sh,
                                        psoils[i_root])
                              );
        node.plant_hs.roots[i_root].p_ups = psoils[i_root];
    end

    # update root flow rates
    roots_flow!(node.plant_hs, e_sum);

    return nothing
end

function test_soil_from_swc!(
            node::SPACMono{FT},
            swcs::Array{FT,1},
            e_sum::FT
) where {FT<:AbstractFloat}
    # update the soil water contents and potential in each layer
    for i_root in 1:node.n_root
        node.swc[i_root] = max(node.mswc[i_root], swcs[i_root]);
        node.p_soil[i_root] = soil_p_25_swc(node.plant_hs.roots[i_root].sh,
                                            node.swc[i_root]) *
                              node.plant_hs.roots[i_root].f_st;
        # pass the pressure to root_hs
        node.plant_hs.roots[i_root].p_ups = node.p_soil[i_root];
    end

    # update root flow rates
    roots_flow!(node.plant_hs, e_sum);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Wayland_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Wayland")
JLLWrappers.@generate_main_file("Wayland", UUID("a2964d1f-97da-50d4-b82a-358c7fce9d89"))
end  # module Wayland_jll


text:
# script to plot model predictions for the interleaved dataset
#
# The script can be called in multiple ways:
#
# (i) discover the model parameters
#
# plotmodelitlvd.jl [model]
#
# (ii) plot model fits for previously fitted model
#
# plotmodelitlvd.jl [model] [objective]
#
# (iii) plot model fits for given parameter values
#
# plotmodelitlvd.jl [model] [p1] [p2] ...
#
# The plots are written to figs/[model]_[objective]_itlvd.pdf

using Random, Printf, CSV
import Fontconfig, Cairo

include("common.jl")
include("simmodel.jl")
include("behaviorstats.jl")
include("models.jl")

const nrats = 1
const t = Theme(minor_label_font_size=6pt, major_label_font_size=8pt)
const fitpath = "fits"
const datapath = "../fitdata/$(fitpath)"
const figpath = "../fitdata/$(fitpath)"

# returns model for given name
function createmodel(modelname)
    modeldict = Dict("opt" => OptimalModel,
                     "optlapse" => OptimalLapseModel,
                     "optcollapse" => OptimalColLapseModel,
                     "optbiaslapse" => OptimalBiasLapseModel,
                     "optcolbiaslapse" => OptimalColBiasLapseModel,
                     "adflapse" => ADFLapseModel,
                     "adfcollapse" => ADFColLapseModel,
                     "adfbiaslapse" => ADFBiasLapseModel,
                     "adfcolbiaslapse" => ADFColBiasLapseModel,
                     "gammalapse" => GammaLapseModel,
                     "gammacollapse" => GammaColLapseModel,
                     "gammabiaslapse" => GammaBiasLapseModel,
                     "gammacolbiaslapse" => GammaColBiasLapseModel,
                     "deltalapse" => DeltaLapseModel,
                     "deltacollapse" => DeltaColLapseModel,
                     "deltabiaslapse" => DeltaBiasLapseModel,
                     "deltacolbiaslapse" => DeltaColBiasLapseModel,
                     "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                     "lcalapse" => LCALapseModel,
                     "lcacollapse" => LCAColLapseModel)
    if modelname in keys(modeldict)
        return modeldict[modelname]()
    else
        error("Uknown model identifier $modelname")
    end
end

function plotitlvddataset(m_itlvd, plotfilename)
    mixctrst = [0.12, 0.36, 0.6, 1.0]
    nconcentr = 4
    d = load_data(:itlvd)
    # generate plot layers, iterating over concetrations
    lchron = Array{Gadfly.Layer}(undef, nconcentr, 4)
    lpsych = Array{Gadfly.Layer}(undef, nconcentr, 4)
    for i = 1:nconcentr  # iterate over different concentrations
        # stats for given concentration
        _, d_bymag = getpsychchronstats(d[d[:concentr] .== i,:], nrats)
        _, m_bymag = getpsychchronstats(m_itlvd[m_itlvd[:concentr] .== i,:], nrats)
        d_bymag[:tmin] = d_bymag[:tμ] - d_bymag[:tsem]
        d_bymag[:tmax] = d_bymag[:tμ] + d_bymag[:tsem]
        m_bymag[:tmin] = m_bymag[:tμ] - m_bymag[:tsem]
        m_bymag[:tmax] = m_bymag[:tμ] + m_bymag[:tsem]
        d_bymag[:pcorrmin] = d_bymag[:pcorrμ] - d_bymag[:pcorrsem]
        d_bymag[:pcorrmax] = d_bymag[:pcorrμ] + d_bymag[:pcorrsem]
        m_bymag[:pcorrmin] = m_bymag[:pcorrμ] - m_bymag[:pcorrsem]
        m_bymag[:pcorrmax] = m_bymag[:pcorrμ] + m_bymag[:pcorrsem]
        # the 8i-36 turns the stimmag into the (1-4) range
        d_bymag[:stimctrst] = mixctrst[d_bymag[:stimmag] .+ (8i - 36)]
        m_bymag[:stimctrst] = mixctrst[m_bymag[:stimmag] .+ (8i - 36)]
        ti = copy(t)
        ti.default_color = RGB(([1.0, 1.0, 1.0] .- 0.2i)...)
        tr = copy(ti)
        tr.lowlight_color = c->RGBA([[1.0, 1.0, 1.0] .- 0.2i; [0.2]]...)  # shaded color
        lchron[i,1] = layer(d_bymag, x="stimctrst", y="tμ", Geom.point, ti)[1]
        lchron[i,2] = layer(d_bymag, x="stimctrst", y="tμ", ymin="tmin", ymax="tmax",
                            Geom.errorbar, ti)[1]
        lchron[i,3] = layer(m_bymag, x="stimctrst", y="tμ", Geom.line, ti)[1]
        lchron[i,4] = layer(m_bymag, x="stimctrst", ymin="tmin", ymax="tmax",
                            Geom.ribbon, tr)[1]
        lpsych[i,1] = layer(d_bymag, x="stimctrst", y="pcorrμ", Geom.point, ti)[1]
        lpsych[i,2] = layer(d_bymag, x="stimctrst", y="pcorrμ", ymin="pcorrmin", ymax="pcorrmax",
                            Geom.errorbar, ti)[1]
        lpsych[i,3] = layer(m_bymag, x="stimctrst", y="pcorrμ", Geom.line, ti)[1]
        lpsych[i,4] = layer(m_bymag, x="stimctrst", ymin="pcorrmin", ymax="pcorrmax",
                            Geom.ribbon, tr)[1]
    end
    # assemble plots
    pchronid = plot(lchron[:,1], lchron[:,2], lchron[:,3], lchron[:,4],
        Guide.xlabel("concentration contrast", orientation=:horizontal),
        Guide.ylabel("reaction time [s]", orientation=:vertical),
        Guide.title("Chronometric curve, interleaved"),
        Scale.x_continuous(minvalue=0.0, maxvalue=1.0),
        Scale.y_continuous(minvalue=0.26, maxvalue=0.44))
    ppsychid = plot(lpsych[:,1], lpsych[:,2], lpsych[:,3], lpsych[:,4],
        Guide.xlabel("concentration contrast", orientation=:horizontal),
        Guide.ylabel("correct choiced (fraction)", orientation=:vertical),
        Guide.title("Psychometric curve, interleaved"),
        Scale.x_continuous(minvalue=0.0, maxvalue=1.0),
        Scale.y_continuous(minvalue=0.5, maxvalue=1.0))
    p = vstack(pchronid, ppsychid)
    draw(PDF(plotfilename, 8.5inch, 11inch), p)
end

# process arguments
if isempty(ARGS)
    error("Plotting script needs to be called with at least on argument")
end
const modelname = lowercase(ARGS[1])
const m = createmodel(modelname)
const ϕns = ϕnames(m)
const ϕn = length(ϕns)
if length(ARGS) == 1
    println("Model $modelname has $ϕn parameters:")
    println(foldl((a, b) -> a * ", " * b, [string(i) for i in ϕns]))
    quit()

elseif length(ARGS) == 2
    basefilename = "$(modelname)_$(ARGS[2])"
    fitfilename = "$(datapath)/fit_$(basefilename).csv"
    if !isfile(fitfilename)
        error("Fitting data $fitfilename not found")
    end
    println("Loading fitting data from $fitfilename")
    d = CSV.read(fitfilename)
    (maxllh, maxi) = findmax(d[:llh])
    if maxi != size(d,1) && maxllh > d[end,end]
        warn("Maximum llh $maxllh found in line $maxi instead of last line")
    end
    ϕ = Vector(d[maxi,1:end-1])
else
    basefilename = modelname
    if length(ARGS) != ϕn+1
        error("Expected $(ϕn+1) parameters for model $modelname")
    end
    # process and print arguments
    ϕ = [float(s) for s in ARGS[2:end]]
end
println("Plotting model $modelname with parameters")
println(foldl((a, b) -> a * ", " * b,
        [@sprintf("%s=%.2f", ϕns[i], ϕ[i]) for i in 1:ϕn]))

# simulate model
Random.seed!(0)
const evidence, ids = catitlvdevidence(itlvd_cs, 10, 10^5)
const m_itlvd = simitlvdmodel(m, ϕ, evidence, ids)

# plot itlvd condition
plotitlvddataset(m_itlvd, "$(figpath)/$(basefilename)_itlvd.pdf")
println("Plot written to $(figpath)/$(basefilename)_itlvd.pdf")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   Pre-defined formatters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export ft_printf, ft_round, ft_latex_sn

"""
    function ft_printf(ftv_str, [columns])

Apply the formats `ftv_str` (see the function `sprintf1` of the package
**Formatting.jl**) to the elements in the columns `columns`.

If `ftv_str` is a `Vector`, then `columns` must be also be a `Vector` with the
same number of elements. If `ftv_str` is a `String`, and `columns` is not
specified (or is empty), then the format will be applied to the entire table.
Otherwise, if `ftv_str` is a `String` and `columns` is a `Vector`, then the
format will be applied only to the columns in `columns`.

# Remarks

This formatter will be applied only to the cells that are of type `Number`.

"""
ft_printf(ftv_str::String) = ft_printf([ftv_str])
ft_printf(ftv_str::String, column::Int) = ft_printf(ftv_str, [column])
ft_printf(ftv_str::String, columns::AbstractVector{Int}) =
    ft_printf([ftv_str for i = 1:length(columns)], columns)

function ft_printf(ftv_str::Vector{String}, columns::AbstractVector{Int} = Int[])
    lc = length(columns)

    lc == 0 && (length(ftv_str) != 1) &&
    error("If columns is empty, then ftv_str must have only one element.")

    lc > 0 && (length(ftv_str) != lc) &&
    error("The vector columns must have the same number of elements of the vector ftv_str.")

    # By using the function `sprintf1` from the package Formatting.jl, it was
    # possible to reduce the time to print a 100x5 table of `Float64`s from 25s
    # to 0.04s. Thanks to @RalphAS for the tip!

    if lc == 0
        return Dict{Int,Function}(0 => (v,i) -> begin
            return typeof(v) <: Number ? sprintf1(ftv_str[1], v) : v
        end)
    else
        ft = Dict{Int,Function}()

        for i = 1:length(columns)
            push!(ft, columns[i] => (v,j) -> begin
                return typeof(v) <: Number ? sprintf1(ftv_str[i], v) : v
            end)
        end

        return ft
    end
end

"""
    function ft_round(digits, [columns])

Round the elements in the columns `columns` to the number of digits in `digits`.

If `digits` is a `Vector`, then `columns` must be also be a `Vector` with the
same number of elements. If `digits` is a `Number`, and `columns` is not
specified (or is empty), then the rounding will be applied to the entire table.
Otherwise, if `digits` is a `Number` and `columns` is a `Vector`, then the
elements in the columns `columns` will be rounded to the number of digits
`digits`.

"""
ft_round(digits::Int) = ft_round([digits])
ft_round(digits::Int, columns::AbstractVector{Int}) =
    ft_round([digits for i = 1:length(columns)], columns)

function ft_round(digits::AbstractVector{Int}, columns::AbstractVector{Int} = Int[])
    lc = length(columns)

    lc == 0 && (length(digits) != 1) &&
    error("If columns is empty, then digits must have only one element.")

    lc > 0 && (length(digits) != lc) &&
    error("The vector columns must have the same number of elements of the vector digits.")

    if lc == 0
        return Dict{Int,Function}(0 => (v,i) -> begin
                                      if applicable(round,v)
                                          return round(v; digits = digits[1])
                                      else
                                          return v
                                      end
                                  end)
    else
        ft = Dict{Int,Function}()

        for i = 1:length(columns)
            push!(ft, columns[i] => (v,j) -> begin
                      if applicable(round,v)
                          return round(v; digits = digits[i])
                      else
                          return v
                      end
                  end)
        end

        return ft
    end
end

"""
    function ft_latex_sn(m_digits, [columns])

Format the numbers of the elements in the columns `columns` to a scientific
notation using LaTeX. The number is first printed using `sprintf1` functions
with the `g` modifier and then converted to the LaTeX format. The number of
digits in the mantissa can be selected by the argument `m_digits`.

If `m_digits` is a `Vector`, then `columns` must be also be a `Vector` with the
same number of elements. If `m_digits` is a `Integer`, and `columns` is not
specified (or is empty), then the format will be applied to the entire table.
Otherwise, if `m_digits` is a `String` and `columns` is a `Vector`, then the
format will be applied only to the columns in `columns`.

# Remarks

This formatter will be applied only to the cells that are of type `Number`.

"""
ft_latex_sn(m_digits::Int) = ft_latex_sn([m_digits])
ft_latex_sn(m_digits::Int, columns::AbstractVector{Int}) =
    ft_latex_sn([m_digits for i = 1:length(columns)], columns)

function ft_latex_sn(m_digits::AbstractVector{Int}, columns::AbstractVector{Int} = Int[])
    lc = length(columns)

    lc == 0 && (length(m_digits) != 1) &&
    error("If columns is empty, then `m_digits` must have only one element.")

    lc > 0 && (length(m_digits) != lc) &&
    error("The vector columns must have the same number of elements of the vector `m_digits`.")

    if lc == 0
        return Dict{Int,Function}(0 => (v,i) -> begin
            str::String = ""
            if typeof(v) <: Number
                str = sprintf1("%." * string(m_digits[1]) * "g", v)

                # Check if we have scientific notation.
                aux = match(r"e[+,-][0-9]+", str)
                if aux != nothing
                    exp_str = " \\cdot 10^{" * string(parse(Int,aux.match[2:end])) * "}"
                    str = replace(str, r"e.*" => exp_str)
                    str = "\$" * str * "\$"
                end
            else
                str = sprint(print, v)
            end

            return str
        end)
    else
        ft = Dict{Int,Function}()

        for i = 1:length(columns)
            push!(ft, columns[i] => (v,j) -> begin
                str::String = ""

                if typeof(v) <: Number
                    str = sprintf1("%." * string(m_digits[i]) * "g", v)

                    # Check if we have scientific notation.
                    aux = match(r"e[+,-][0-9]+", str)
                    if aux != nothing
                        exp_str = " \\cdot 10^{" * string(parse(Int,aux.match[2:end])) * "}"
                        str = replace(str, r"e.*" => exp_str)
                        str = "\$" * str * "\$"
                    end
                else
                    str = sprint(print, v)
                end

                return str
            end)
        end

        return ft
    end
end


text:
#========================================================================================#
"""
	SimpleGAs

A simple Genetic Algorithm implementation.

Author: Niall Palfreyman, 13/3/2022.
"""
module SimpleGAs

export SimpleGA, mu!, px!, elite!, temperature!, init!, depict, evolve!

include("./Objectives.jl")
include("./Decoders.jl")
include("../Casinos/Casinos.jl")

using .Objectives, .Decoders, .Casinos
using Statistics: std, mean
using Random: shuffle!

#-----------------------------------------------------------------------------------------
# Module types:

"""
	SimpleGA

A very basic GA with the capability of optimising binary-encoded functions. This implementation
uses mutation, single-point crossover, elitism and sigma-scaling.
"""
mutable struct SimpleGA
	ngenes				# Number of genes in each individual
	npop				# Number of individuals in population
	population			# Population of individual binary chromosomes
	mu					# Mutation rate
	px					# Crossover rate
	temperature			# Sigma-scaling coeff: -ve values maximise objective function
	nNonelitist			# Number of non-elite recombinable individuals
	objective			# Objective function
	decoder				# Binary decoding for evaluation
	casino				# Mutation repository (Casino)

	"The one and only SimpleGA constructor"
	function SimpleGA( obj, npop, nbits::Int=-1)
		if rem(npop,2)==1
			# Warning: npop must be even for our crossover algorithm!
			npop = npop + 1
		end
		ngenes = nbits * dim(obj)
		decoder = nbits > 0 ? Decoder(obj.domain,nbits) : nothing

		new(
			ngenes,									# ngenes
			npop,									# npop
			[rand(Bool,ngenes) for _ in 1:npop],	# population
			2/(ngenes*npop),						# mu: About 2 per generation
			1,										# px
			1,										# temperature
			npop,									# nNonelitist: No elitism
			obj,									# objective
			decoder,								# decoder
			Casino(ngenes,npop)						# casino
		)
	end
end

#-----------------------------------------------------------------------------------------
# Module methods:

"""
	mu!(sga,mu)

Set the mutation probability.
"""
function mu!( sga, mu)
	sga.mu = max(min(mu,1),0)
end

#-----------------------------------------------------------------------------------------
"""
	px!(sga,px)

Set the crossover probability.
"""
function px!( sga, px)
	sga.px = max(min(px,1),0)
end

#-----------------------------------------------------------------------------------------
"""
	elite!(sga,elite)

Set number of elite individuals that are to be retained into next generation. If elite<1,
it describes a proportion of elite individuals; if elite>=1, it is the number of elites.
nNonelitist is the number of non-elite individuals - implemented this way for performance
reasons.
"""
function elite!( sga, elite)
	if elite < 1
		elite = ceil( Int, sga.npop * max(min(elite,1),0))
	else
		elite = floor( Int, elite)
	end
	if rem(elite,2) == 1
		# Recombined population must remain even!
		elite = elite + 1
	end
	sga.nNonelitist = sga.npop - elite
end

#-----------------------------------------------------------------------------------------
"""
	temperature!(sga,sigma)

Set the temperature of the algorithm. High temperature means the GA considers many poor
solutions, while low temperature means the GA only accepts high-quality solutions.
"""
function temperature!( sga, theta)
	sga.temperature = theta
end;

#-----------------------------------------------------------------------------------------
"""
	init!(sga)

Reinitialise SimpleGA population.
"""
function init!( sga)
	sga.population = [rand(0:1,ngenes) for _ in npop]
end

#-----------------------------------------------------------------------------------------
"""
	depict(sga)

Display current status of population.
"""
function depict( sga)
	fit,evals = fitness(sga)
	_,maxi = findmax(fit)
	println( "Avg value:       ", mean(evals))
	println( "Best value:      ", evals[maxi])
	println( "Best individual: ", individuals[maxi])
end

#-----------------------------------------------------------------------------------------
"""
	evolve!(sga)

Evolve the population through ngen generations.
"""
function evolve!( sga, ngen::Int=1)
	for _ in 1:ngen
		fitnesses,_ = fitness(sga)

		if sga.nNonelitist < sga.npop
			# Apply elitism - sort elite individuals to
			# end of population list:
			sp = sortperm(fitnesses)
			fitnesses = fitnesses[sp]
			sga.population = sga.population[sp]
		end

		if sga.px > 0
			# Perform crossover:
			recombine!(sga,fitnesses)
		end

		if sga.mu > 0
			# Perform mutation:
			mutate!(sga)
		end
	end

end

#-----------------------------------------------------------------------------------------
"""
	fitness(sga)

Calculate normalised fitness of the population. fit is a column vector of normalised
fitnesses of sga population, minus all sub-sigma-scaled individuals (see Mitchell p.168).
Negative sigma-scaling maximises the objective function; higher magnitudes raise the
fitness pressure. objeval is a colum vector of evaluations of the population. indiv is a
matrix of decoded values of population individuals.
"""
function fitness( sga)
	if sga.decoder === nothing
		# No decoding required
		individuals = sga.population
	else
		individuals = map(sga.population) do x
			decode(sga.decoder,x)
		end
	end

	# Calculate the current objective evaluations of the population:
	evaluations = evaluate(sga.objective,individuals)
	# Normalise the evaluations into frequencies:
	sigma = std(evaluations)					# Standard deviation
	if sigma == 0
		# Singular case: all evaluations were equal to the mean:
		fitness = ones(sga.npop)
	else
		# Exorcise all evaluations worse than 1/temperature standard
		# deviations above mean value:
		fitness = 1 .+ (mean(evaluations) .- evaluations) ./ (sga.temperature*sigma)
		fitness[fitness .<= 0] .= 0
	end

	(fitness,evaluations,individuals)					# Return value
end

#-----------------------------------------------------------------------------------------
"""
	recombine!(sga)

Select and recombine members of the sga population. We use Stochastic Universal Sampling
(SUS - see Mitchell p. 168) to select the prospective parent indices. SUS is really cool
and clever, and almost always improves performance.
"""
function recombine!( sga::SimpleGA, fitnesses=nothing)
	# Only recalculate fitnesses if they aren't provided:
	if fitnesses === nothing
		fitnesses,_ = fitness(sga)
	end

	# Set up a roulette wheel of fitness-biased slots:
	roulette = cumsum( fitnesses)/sum(fitnesses)
	# Set up npop (rigidly) uniformly distributed ball throws (SUS):
	throws = rem.( rand() .+ (1:sga.npop)./sga.npop, 1)

	# Throw ball onto roulette wheel to create randomised list of
	# fitness-selected parent indices:
	parents = zeros(Integer,sga.npop)
	for parent in 1:sga.npop
		for slot in 1:sga.npop
			if throws[parent] <= roulette[slot]
				parents[parent] = slot
				break
			end
		end
	end
	shuffle!(parents)				# Shuffle the order of parents

	# 1st half of parents are Mummies; 2nd half are Daddies:
	nMatings = sga.nNonelitist ÷ 2
	mummy = sga.population[parents[1:nMatings]]
	daddy = sga.population[parents[nMatings+1:sga.nNonelitist]]

	# Create next generation:
	progeny = similar(sga.population)
	for m in 1:nMatings
		xpt = rand(1:length(mummy[m])-1)			# Crossover point

		sally = deepcopy.(mummy[m])					# mummy-based progeny
		billy = deepcopy.(daddy[m])					# daddy-based progeny

		# To do: Now transfer the second part of daddy to sally, and the second
		# part of mummy to billy:
		# ???										# Sally also takes after daddy
		# ???										# Billy also takes after mummy

		# To do: Now store sally as m-th progeny and billy as (m+nMatings)-th progeny:
		# ???
		# ???
	end

	# To do: ... and finally replace the old fogeys by the new progeny:
	# ???
	sga
end

#-----------------------------------------------------------------------------------------
"""
	mutate!(sga)

Mutate the (non-elite) SimpleGA population Draw mutation masks from a pregenerated casino.
This significantly accelerates code by avoiding many random number generations, with no
apparent impact on GA behaviour.
"""
function mutate!( sga)
	if sga.mu > 0
		# Only expend effort on mutating if it is really wanted:

		# You need to implement the new Casinos.draw() method to run this next line of code:
		# Find loci to mutate (but ignore the elitists):
		loci = draw( sga.casino, sga.ngenes, sga.nNonelitist, sga.mu)

		# Now mutate everyone except the elitists:
		for i ∈ 1:sga.nNonelitist
			# Negate the mutated genes:
			sga.population[i][loci[:,i]] = .!(sga.population[i][loci[:,i]])
		end
	end

	sga
end

#-----------------------------------------------------------------------------------------
"""
	demo()

Demo the SimpleGA type.
"""
function demo()
	obj = Objective(7)
	sga = SimpleGA(obj,128,20)
	temperature!(sga,2)
	depict(sga)

	nIter = 5000
	println( "\nAfter $nIter iterations ...")
	evolve!(sga,nIter)
	depict(sga)
end

end		# of SimpleGAs


text:
# extensions to JuMP for multiobjective optimisation

module MultiJuMP

using JuMP
import JuMP: JuMPTypes, getvalue
import MathProgBase
using RecipesBase, LaTeXStrings
import Combinatorics
using LinearAlgebra
export multi_model, SingleObjective, get_multidata
export getutopia, getnadir
export WeightedSum, EpsilonCons, NBI

include("types.jl")
include("methods.jl")
include("linear.jl")
include("nonlinear.jl")
include("plot_recipe.jl")

end


text:
function optiEdificio(dcn, dca, dcp, dcc, dcu, dcf, dcr, alturaPso, ps_base, superficieTerreno, superficieTerrenoBruta)

    areaBasalPso = polyShape.polyArea_v2(ps_base)
    numTiposDepto = length(dcc.SUPDEPTOUTIL);
    superficieDensidad = dcn.FLAGDENSIDADBRUTA ? superficieTerrenoBruta : superficieTerreno
    maxDeptos = dcn.DENSIDADMAX / 4 * superficieDensidad / 10000;

    ##############################################
    # PARTE "3": DEFINICIÓN SOLVER               #
    ##############################################

    m = Model(Cbc.Optimizer)
    set_optimizer_attribute(m, "ratioGap", 0.001)
    set_optimizer_attribute(m, "threads", 3)

    #m = Model(GAMS.Optimizer)
    #set_optimizer_attribute(m, "OptCR", 0.001)
    #set_optimizer_attribute(m, "Threads", 3)
    #set_optimizer_attribute(m, "logOption", 0)

    ##############################################
    # PARTE "4": VARIABLES DE DECISION           #
    ##############################################

    @variables(m, begin
        0 <= numPisos <= dcn.MAXPISOS
        0 <= numDeptosTipo[u = 1:numTiposDepto]
        0 <= CostoUnitTerreno
    end)

    ##############################################
    # PARTE "5": EXPRESIONES AUXILIARES          #
    ##############################################

    # Cálculo Número de Estacionamientos
    estacionamientosViviendas = sum(dcn.ESTACIONAMIENTOSPORVIV .* numDeptosTipo)
    estacionamientosVisitas = estacionamientosViviendas * dcn.PORCADICESTACVISITAS;
    estacionamientosNormales = estacionamientosViviendas + estacionamientosVisitas
    estacionamientosDiscapacitados = (maxDeptos <= 20) ? 1 :
                                     ((maxDeptos <= 50) ? 2 :
                                     ((maxDeptos <= 200) ? 3 :
                                     ((maxDeptos <= 400) ? 4 :
                                     ((maxDeptos <= 500) ? 5 : (0.01*maxDeptos)))))
    estacionamientosBicicletas = estacionamientosNormales * dcn.ESTBICICLETAPOREST;
    descuentoEstCercaniaMetro = estacionamientosNormales * 0.5*dcn.REDUCCIONESTPORDISTMETRO
    descuentoEstBicicletas = estacionamientosBicicletas/dcn.BICICLETASPOREST
    estacionamientosNormales = estacionamientosNormales - descuentoEstCercaniaMetro - descuentoEstBicicletas
    cambioEstBicicletas = (dcn.FLAGCAMBIOESTPORBICICLETA) ? estacionamientosNormales / 3 : 0
    estacionamientosVendibles = estacionamientosViviendas - descuentoEstCercaniaMetro - descuentoEstBicicletas - cambioEstBicicletas
    estacionamientosBicicletas = estacionamientosBicicletas + cambioEstBicicletas * dcn.BICICLETASPOREST
    estacionamientos = estacionamientosVendibles + estacionamientosVisitas + estacionamientosDiscapacitados + estacionamientosBicicletas / dcn.BICICLETASPOREST;

    # Cálculo de Superficies
    superficieUtilDepto = numDeptosTipo .* dcc.SUPDEPTOUTIL;
    superficieUtil = sum(superficieUtilDepto);
    superficieTerrazaDepto = 2 .* superficieUtilDepto .* dca.PORCTERRAZA;
    superficieInteriorDepto = superficieUtilDepto .* (1 - dca.PORCTERRAZA);
    superficieVendibleDepto = superficieInteriorDepto + .5 * superficieTerrazaDepto; # = superficieUtilDepto
    superficieVendible = sum(superficieVendibleDepto);
    superficieComun = dca.PORCSUPCOMUN * superficieUtil;
    superficieLosaSNT = superficieVendible + superficieComun;
    superficieLosaBNT = dcn.SUPPORESTACIONAMIENTO * estacionamientosVendibles;

    # Cálculo de Ingresos por Ventas
    ingresosVentaDeptos = sum(superficieVendibleDepto .* dcc.PRECIOVENTA);
    ingresosVentaEstacionamientos = estacionamientosVendibles * dcc.PRECIOVENTAEST;
    IngresosVentas = ingresosVentaDeptos + ingresosVentaEstacionamientos;

    # Cálculo de Costos de Terreno y Construcción
    CostoTerreno = CostoUnitTerreno * superficieTerreno * (1 + dcu.ComisionCorredor);
    CostoConstruccionSNT = dcu.LosaSNT * superficieLosaSNT; # Costos con IVA
    CostoConstruccionBNT = dcu.LosaBNT * superficieLosaBNT; # Costos con IVA
    CostoEstacVisita = dcu.EstacVisita * dcn.SUPPORESTACIONAMIENTO * estacionamientosVisitas; # Costos con IVA
    devolucionIVA = (CostoConstruccionSNT + CostoConstruccionBNT + CostoEstacVisita) / 1.19 * 0.19 * 0.65; # Divide por 1.19 para dejar en base sin IVA
    CostoConstruccion = dcu.Demolicion * superficieTerreno + CostoConstruccionSNT + CostoConstruccionBNT + CostoEstacVisita +
                    dcu.ExtrasPostVenta * (CostoConstruccionSNT + CostoConstruccionBNT) - devolucionIVA;

    # Cálculo de Costos de Permisos, Habilitación, Puesta en Marcha y Atención de Clientes
    CostoInmobiliarioObra = dcu.EmpalmesAportes * sum(numDeptosTipo) + dcu.DerechosPermisosMunicipales +
                    dcu.duracionProyecto * dcu.InspeccionCertificacionTecnica + dcu.Contribuciones * (dcu.duracionProyecto / 12) * 4 +
                    dcu.CopiaPlanosArquitctura + dcu.CopiaPlanosObra;
    CostosHabilitacion = dcu.HallAcceso + dcu.JardinesPaisajismo + dcu.HabilitacionSalaUsosMultiples + dcu.HabilitacionPiscinaEnfermeria;
    CostosPuestaEnMarcha = dcu.duracionProyecto * dcu.CuentasServiciosDeptos + dcu.duracionProyecto * dcu.GastosComunes +
                    dcu.ContribucionesDeptos + dcu.Seguros;
    CostosAtencionCliente = dcu.ResponsabilidadCivil * CostoConstruccion + dcu.PostVentaInmobiliaria * CostoConstruccion +
                    dcu.SeguroVentaEnVerde * CostoConstruccion;

    # Cálculo de Honorarios de Proyecto, Administración y Gestión de Venta
    CostoHonorariosProyectos = dcu.Arquitectura * ingresosVentaDeptos + dcu.GestionAdministracion * ingresosVentaDeptos +
                    dcu.duracionProyecto * (dcu.ITO + dcu.Contabilidad) + dcu.Legales * CostoConstruccion +
                    dcu.Calculo * (superficieLosaSNT + superficieLosaBNT) + dcu.MecanicaSuelo + dcu.Topografia +
                    dcu.ProyectoElectrico * CostoConstruccion + dcu.ProyectoSanitario * CostoConstruccion +
                    dcu.ProyectoSeguridad * CostoConstruccion + dcu.ProyectoPavimentacion + dcu.ProyectoBasura + dcu.ProyectoCalefaccion;
    CostoPromocionVenta = dcu.Marketing * IngresosVentas + dcu.PilotSalaVentas * IngresosVentas;
    CostoOtrosVenta = dcu.NotarialesLegales * sum(numDeptosTipo) + dcu.CBR * sum(numDeptosTipo);
    CostoVenta = CostoPromocionVenta + CostoOtrosVenta;

    # Cálculo de Otros Costos
    Imprevistos = dcu.Imprevistos * IngresosVentas;
    DebitoIVA = (IngresosVentas - CostoTerreno) / 1.19 * 0.19;
    CreditoIVA = (CostoConstruccion + CostoVenta + CostosHabilitacion) / 1.19 * 0.19;
    CostoPagoIVA = DebitoIVA - CreditoIVA;

    # Cálculo Costo Total
    CostoTotal = CostoTerreno + CostoConstruccion + CostoHonorariosProyectos +
                 CostoVenta + CostoInmobiliarioObra + CostosHabilitacion +
                 CostosPuestaEnMarcha + CostosAtencionCliente + Imprevistos + CostoPagoIVA;

    ##############################################
    # PARTE "6": RESTRICCIONES DEL MIP           #
    ##############################################

    @constraints(m, begin
    # Restricción de Altura Máxima y Área Basal Máxima (Coeficiente de Ocupación)
        numPisos * dca.ALTURAPISO <= alturaPso

    # Restricciones que establecen relaciones entre areaBasal, numDeptos, largo y ancho
        superficieLosaSNT <= areaBasalPso * numPisos
        maxCoefConstructibilidad, superficieUtil <= superficieTerreno * dcn.COEFCONSTRUCTIBILIDAD * (1 + 0.3 * dcp.FUSIONTERRENOS)

    # Restricciones de Rentabilidad Mínima
        restRentabilidadMin, IngresosVentas >= dcr.RetornoExigido * CostoTotal

    # Restricción Suma de Departamentos menor a Máximo Número de Departamentos
        RestDensidadMax, sum(numDeptosTipo) <= maxDeptos

    # Restricciones de Demanda
        RestDemanda, numDeptosTipo .<= maxDeptos .* dcc.MAXPORCTIPODEPTO

    end)

    ##############################################
    # PARTE "7": FUNCIÓN OBJETIVO Y EJECUCIÓN    #
    ##############################################

    @objective(m, Max, CostoUnitTerreno)
    # @objective(m, Max, sum(numDeptosTipo))
    # @objective(m, Max, IngresosVentas-CostoTotal)

    # Resuelve el problema de optimización
    JuMP.optimize!(m)

    if termination_status(m) == MOI.OPTIMAL
        status = true

    ##############################################
    # PARTE "8": PRESENTACION DE RESULTADOS      #
    ##############################################

        # Cálculo del Flujo de Caja y la Tir
        fl_IngresosVentas = dcf.IngresosVentas .* JuMP.value(IngresosVentas);
        fl_CostoTerreno = dcf.CostoTerreno .* JuMP.value(CostoTerreno);
        fl_CostoConstruccion = dcf.CostoConstruccion .* JuMP.value(CostoConstruccion);
        fl_CostoHonorariosProyectos = dcf.CostoHonorariosProyectos .* JuMP.value(CostoHonorariosProyectos);
        fl_CostoVenta = dcf.CostoVenta .* JuMP.value(CostoVenta);
        fl_CostoInmobiliarioObra = dcf.CostoInmobiliarioObra .* JuMP.value(CostoInmobiliarioObra);
        fl_CostosHabilitacion = dcf.CostosHabilitacion .* CostosHabilitacion;
        fl_CostosPuestaEnMarcha = dcf.CostosPuestaEnMarcha .* CostosPuestaEnMarcha;
        fl_CostosAtencionCliente = dcf.CostosAtencionCliente .* JuMP.value(CostosAtencionCliente);
        fl_Imprevistos = dcf.Imprevistos .* JuMP.value(Imprevistos);
        fl_CostoPagoIVA = fl_IngresosVentas .* (1 - JuMP.value(CostoTerreno) / sum(fl_IngresosVentas)) .* (0.19 / 1.19) -
                    (fl_CostoConstruccion + fl_CostoVenta + fl_CostosHabilitacion) .* (0.19 / 1.19);
        fl_FlujoCajaNetoAntesImpuesto = fl_IngresosVentas - fl_CostoTerreno - fl_CostoConstruccion - fl_CostoHonorariosProyectos - fl_CostoVenta -
                                fl_CostoInmobiliarioObra - fl_CostosHabilitacion - fl_CostosPuestaEnMarcha - fl_CostosAtencionCliente -
                                fl_Imprevistos - fl_CostoPagoIVA;
        tirAntesImpuesto = (NPFinancial.irr(fl_FlujoCajaNetoAntesImpuesto) + 1)^2 - 1;

        numFlujos = length(fl_FlujoCajaNetoAntesImpuesto);
        fl_FlujoCajaNetoAntesImpuestos_acum = cumsum(fl_FlujoCajaNetoAntesImpuesto);
        vecCorrel = 1:numFlujos;
        flag = cumsum(fl_FlujoCajaNetoAntesImpuestos_acum .> 0) .== 1;
        flagCorrel = vecCorrel .* flag;
        posNoNeg = findall(x->x > 0, flagCorrel)
        fl_ImpuestoRenta = zeros(Float64, (numFlujos,));
        fl_ImpuestoRenta[posNoNeg] = fl_FlujoCajaNetoAntesImpuestos_acum[posNoNeg] * dcf.TasaImpuestoRenta;
        fl_FlujoCajaNetoDespuesImpuesto = fl_FlujoCajaNetoAntesImpuesto - fl_ImpuestoRenta;
        tirDespuesImpuesto = (NPFinancial.irr(fl_FlujoCajaNetoDespuesImpuesto) + 1)^2 - 1;

        fl_LineaCredito = (dcf.ingresoLineaCredito - dcf.pagoLineaCredito) * JuMP.value(CostoConstruccion);
        fl_MontoUtilizadoLineaCredito = cumsum(fl_LineaCredito);
        fl_InteresLineaCredito = fl_MontoUtilizadoLineaCredito * dcf.TasaInteresLineaCredito / 2;
        fl_CostoOperacionalCredito = fl_CostoConstruccion ./ sum(fl_CostoConstruccion) * dcf.DuracionLineaCredito * dcf.Tasaciones;
        fl_FlujoCajaApalancadoNetoAntesImpuesto = fl_FlujoCajaNetoAntesImpuesto + fl_LineaCredito - fl_InteresLineaCredito - fl_CostoOperacionalCredito;
        tirApalancadoAntesImpuestos = (NPFinancial.irr(fl_FlujoCajaApalancadoNetoAntesImpuesto) + 1)^2 - 1;

        # Construye estructura con los resultados de la optimización
        sn = salidaNormativa(
            maxDeptos, # maxNumDeptos
            dcn.COEFOCUPACION * superficieTerreno, # maxOcupacion
            superficieTerreno * dcn.COEFCONSTRUCTIBILIDAD * (1 + 0.3 * dcp.FUSIONTERRENOS), # maxConstructibilidad
            dcn.MAXPISOS, # maxPisos
            dcn.ALTURAMAX, # maxAltura
            JuMP.value(estacionamientosVendibles), # minEstacionamientosVendible
            JuMP.value(estacionamientosVisitas), # minEstacionamientosVisita
            estacionamientosDiscapacitados # minEstacionamientosDiscapacitados
        )

        sa = salidaArquitectonica(
            JuMP.value.(numDeptosTipo), # numDeptosTipo
            sum(JuMP.value.(numDeptosTipo)), # numDeptos
            areaBasalPso, # ocupacion
            JuMP.value(superficieUtil), # constructibilidad
            JuMP.value(numPisos), # numPisos
            JuMP.value(numPisos)*dca.ALTURAPISO, # altura
            JuMP.value(superficieVendible), # superficieUtilSNT
            JuMP.value(superficieComun), # superficieComunSNT
            JuMP.value(superficieVendible) + JuMP.value(superficieComun), # superficieEdificadaSNT
            areaBasalPso, # superficiePorPiso
            JuMP.value(estacionamientosVendibles), # estacionamientosVendibles
            JuMP.value(estacionamientosVisitas), # estacionamientosVisita
            JuMP.value(estacionamientosVendibles) + JuMP.value(estacionamientosVisitas), # numEstacionamientos Totales
            JuMP.value(estacionamientosBicicletas) # numBicicleteros
            )

        si = salidaIndicadores(
            JuMP.value(IngresosVentas), # IngresosVentas
            JuMP.value(CostoTotal), # CostoTotal
            JuMP.value(IngresosVentas) - JuMP.value(CostoTotal), # MargenAntesImpuesto
            tirAntesImpuesto, # TirAntesImpuestos
            dcf.TasaImpuestoRenta * (JuMP.value(IngresosVentas) - JuMP.value(CostoTotal)), # ImpuestoRenta
            JuMP.value(IngresosVentas) - JuMP.value(CostoTotal) - dcf.TasaImpuestoRenta * (JuMP.value(IngresosVentas) - JuMP.value(CostoTotal)), # UtilidadDespuesImpuesto
            (JuMP.value(IngresosVentas) - JuMP.value(CostoTotal)) / JuMP.value(CostoTotal), # RentabilidadTotalBruta
            (JuMP.value(IngresosVentas) - JuMP.value(CostoTotal) - dcf.TasaImpuestoRenta * (JuMP.value(IngresosVentas) - JuMP.value(CostoTotal))) / JuMP.value(CostoTotal), # RentabilidadTotalNeta
            JuMP.value(CostoTerreno) / JuMP.value(IngresosVentas), # IncidenciaTerreno
        )

        st = salidaTerreno(
            superficieTerreno, # superficieTerreno
            superficieTerrenoBruta, # superficieBruta
            0, # superficieEdificacion
            JuMP.value(CostoUnitTerreno * superficieTerreno), # costoTerrenoAntesCorr
            JuMP.value(CostoUnitTerreno), # costoUnitTerrenoAntesCorr
            JuMP.value(CostoUnitTerreno * superficieTerreno * dcu.ComisionCorredor), # CostoCorredor
            JuMP.value(CostoTerreno), # costoTerreno
            JuMP.value(CostoTerreno) / superficieTerreno # costoUnitTerreno
        )

        sm = salidaMonetaria(
            JuMP.value(IngresosVentas), # IngresosVentas
            JuMP.value(CostoTotal), # CostoTotal
            JuMP.value(CostoTerreno), # CostoTerreno
            JuMP.value(CostoTerreno) / superficieTerreno, # CostoUnitarioTerreno
            JuMP.value(CostoConstruccion), # CostoConstruccion
            JuMP.value(CostoHonorariosProyectos), # CostoHonorariosProyectos
            JuMP.value(CostoVenta), # CostoVenta
            JuMP.value(CostoInmobiliarioObra), # CostoInmobiliarioObra
            (CostosHabilitacion), # CostosHabilitacion
            (CostosPuestaEnMarcha), # CostosPuestaEnMarcha
            JuMP.value(CostosAtencionCliente), # CostosAtencionCliente
            JuMP.value(Imprevistos), # Imprevistos
            JuMP.value(CostoPagoIVA) # CostoPagoIVA
        )

        sf = salidaFlujoCaja(
            fl_IngresosVentas, # IngresosVentas
            fl_CostoTerreno, # CostoTerreno
            fl_CostoConstruccion, # CostoConstruccion
            fl_CostoHonorariosProyectos, # CostoHonorariosProyectos
            fl_CostoVenta, # CostoVenta
            fl_CostoInmobiliarioObra, # CostoInmobiliarioObra
            fl_CostosHabilitacion, # CostosHabilitacion
            fl_CostosPuestaEnMarcha, # CostosPuestaEnMarcha
            fl_CostosAtencionCliente, # CostosAtencionCliente
            fl_Imprevistos, # Imprevistos
            fl_CostoPagoIVA, # CostoPagoIVA
            fl_FlujoCajaNetoAntesImpuesto, # flujoCajaNetoAntesImpuestos
            tirAntesImpuesto, # TirAntesImpuestos
            fl_ImpuestoRenta, # ImpuestoRenta
            fl_FlujoCajaNetoDespuesImpuesto, # NetoDespuesImpuestos
            tirDespuesImpuesto, # TirDespuesImpuestos
            fl_LineaCredito, # LineaCredito
            fl_MontoUtilizadoLineaCredito, # MontoUtilizadoLineaCredito
            fl_InteresLineaCredito, # InteresLineaCredito
            fl_CostoOperacionalCredito, # CostoOperacionalCredito
            fl_FlujoCajaApalancadoNetoAntesImpuesto, # ApalancadoNetoAntesImpuesto
            tirApalancadoAntesImpuestos # TirApalancadoAntesImpuestos
        )

    else
        status = false
        sn = nothing
        sa = nothing
        si = nothing
        st = nothing
        sm = nothing
        sf = nothing
    end

    return sn, sa, si, st, sm, sf, status

end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, LinearAlgebra, Plots

begin
    x0 = -1
    x1 = 1
    ncell = 16
    dx = (x1 - x0) / ncell
    deg = 3
    nsp = deg + 1
    cfl = 0.1
    dt = cfl * dx
    t = 0.0
    a = 1.0
    tspan = (0.0, 2.0)
    nt = tspan[2] / dt |> Int
    bc = :period
end
ps = FRPSpace1D(x0, x1, ncell, deg)

u = zeros(ncell, deg + 1)
for i = 1:ncell, j = 1:deg+1
    u[i, j] = sin(π * ps.xpg[i, j])
end

prob = FR.FRAdvectionProblem(u, tspan, ps, a, bc)
itg = init(prob, Tsit5(), saveat = tspan[2], adaptive = false, dt = dt)

for iter = 1:nt
    step!(itg)
end

begin
    x = zeros(ncell * nsp)
    sol = zeros(ncell * nsp)
    for i = 1:ncell
        idx0 = (i - 1) * nsp
        idx = idx0+1:idx0+nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]
            sol[idx] = itg.u[i, j]
        end
    end
    ref = @. sin(π * x)
end

plot(x, sol, label = "t=1")
plot!(x, ref, line = :dash, label = "t=0")

begin
    dx |> println
    FR.L1_error(sol, ref, dx) |> println
    FR.L2_error(sol, ref, dx) |> println
    FR.L∞_error(sol, ref, dx) |> println
end


text:
using DashBootstrapComponents

spinners = html_div([
    dbc_button(
        dbc_spinner(size = "sm"),
        color = "primary",
        disabled = true,
        className = "me-1",
    ),
    dbc_button(
        [dbc_spinner(size = "sm"), " Loading..."],
        color = "primary",
        disabled = true,
    ),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
## permkronecker ##

# Kronecker product for matrices induces a Kronecker product on permutations.
"""
    permkron(p, q)

Compute the Kronecker product of permutations `p` and `q`
induced by their matrix representations.
"""
function permkron(p::PList{T}, q::PList{S}) where {T, S}
    np = length(p); nq = length(q)
    dc = Array{promote_type(T, S)}(undef, np*nq)
@inbounds for i in 1:np
    for k in 1:nq
            dc[nq*(i-1) + k] = nq*(p[i]-1)+q[k]
        end
    end
    dc
end

# This is partly broken
function permkron(pt::PDict{T}, qt::PDict) where {T<:Real}
    dout = Dict{T, T}()
    p = getdata(pt); q = getdata(qt)
    (isempty(p) || isempty(q)) && error("Can't yet do sparse kronecker product with identity permutations")
    np = convert(T,maximum(p)[1])
    nq = convert(T,maximum(q)[1])
    maxk = zero(T)
    #    for (i,pi) in p
    for i in 1:np
        pi = get(p,i,zero(T))
        pi = pi == 0 ? i : pi
        #        for (k,qk) in q
        for k in 1:nq
            qk = get(q,k,zero(T))
            qk = qk == 0 ? k : qk
            ind1 = nq *(pi-1) + qk
            ind2 = nq * (i-1) + k
            ind1 == ind2 ? continue : nothing
            z = get(dout,ind2,zero(T))
            z != zero(T) && continue
            dout[ind2] = ind1
            ind2 > maxk ? maxk = ind2 : nothing
            ind1 > maxk ? maxk = ind1 : nothing  # should be superfluous
        end
    end
    return PDict(dout, permlength(pt) * permlength(qt))
end

# This a bit more efficient than full matrix multiplication. Cuts out one loop.
function permkron(a::PList{T}, b::PMatrix{S}) where {T<:Real,S<:Real}
    (nrowa, ncola) = (length(a),length(a))
    (nrowb, ncolb) = size(b)
    R = zeros(promote_type(T,S), nrowa * nrowb, ncola * ncolb)
    d = invperm(a)
    for j = 1:ncola, l = 1:ncolb
        soff = ncola * (j-1)
        i = d[j]
        roff = ncola * (i-1)
        for k = 1:nrowb
            R[roff+k,soff+l] = b[k,l]
        end
    end
    return R
end

# For testing. If there is a difference, it is small
# function permkron2{T<:Real,S<:Real}(a::AbstractVector{T}, b::AbstractMatrix{S})
#     (nrowa, ncola) = size(a)
#     (nrowb, ncolb) = size(a)
#     R = zeros(promote_type(T,S), nrowa * nrowb, ncola * ncolb)
#     d = a.data
#     for j = 1:ncola, l = 1:ncolb
#         soff = ncola * (j-1)
#         i = d[j]
#         roff = ncola * (i-1)
#         for k = 1:nrowb
#             R[soff+l,roff+k] = b[l,k]
#         end
#     end
#     R
# end


text:
# Autogenerated wrapper script for SuiteSparse_LBT_jll for aarch64-apple-darwin
export libamd, libbtf, libcamd, libccolamd, libcholmod, libcolamd, libklu, libldl, librbio, libspqr, libsuitesparse_wrapper, libsuitesparseconfig, libumfpack

using libblastrampoline_jll
JLLWrappers.@generate_wrapper_header("SuiteSparse_LBT")
JLLWrappers.@declare_library_product(libamd, "@rpath/libamd.2.4.6.dylib")
JLLWrappers.@declare_library_product(libbtf, "@rpath/libbtf.1.2.6.dylib")
JLLWrappers.@declare_library_product(libcamd, "@rpath/libcamd.2.4.6.dylib")
JLLWrappers.@declare_library_product(libccolamd, "@rpath/libccolamd.2.9.6.dylib")
JLLWrappers.@declare_library_product(libcholmod, "@rpath/libcholmod.3.0.13.dylib")
JLLWrappers.@declare_library_product(libcolamd, "@rpath/libcolamd.2.9.6.dylib")
JLLWrappers.@declare_library_product(libklu, "@rpath/libklu.1.3.8.dylib")
JLLWrappers.@declare_library_product(libldl, "@rpath/libldl.2.2.6.dylib")
JLLWrappers.@declare_library_product(librbio, "@rpath/librbio.2.2.6.dylib")
JLLWrappers.@declare_library_product(libspqr, "@rpath/libspqr.2.0.9.dylib")
JLLWrappers.@declare_library_product(libsuitesparse_wrapper, "@rpath/libsuitesparse_wrapper.dylib")
JLLWrappers.@declare_library_product(libsuitesparseconfig, "@rpath/libsuitesparseconfig.5.4.0.dylib")
JLLWrappers.@declare_library_product(libumfpack, "@rpath/libumfpack.5.7.8.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libblastrampoline_jll)
    JLLWrappers.@init_library_product(
        libamd,
        "lib/libamd.2.4.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbtf,
        "lib/libbtf.1.2.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcamd,
        "lib/libcamd.2.4.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libccolamd,
        "lib/libccolamd.2.9.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcholmod,
        "lib/libcholmod.3.0.13.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolamd,
        "lib/libcolamd.2.9.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libklu,
        "lib/libklu.1.3.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libldl,
        "lib/libldl.2.2.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librbio,
        "lib/librbio.2.2.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspqr,
        "lib/libspqr.2.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsuitesparse_wrapper,
        "lib/libsuitesparse_wrapper.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsuitesparseconfig,
        "lib/libsuitesparseconfig.5.4.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libumfpack,
        "lib/libumfpack.5.7.8.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()