text:
"""
    test_criteria StoppingCriterion()

Runs a series of tests to check functionality of a StoppingCriterion

- `update(::StoppingCriterion, loss, state=nothing)` is defined
- Checks that `done(::StoppingCriterion, nothing)` is `false`
- `done(::StoppingCriterion, state)` can handle state after `update` or `update_training`
- `message(::StoppingCriterion, state)` can handle state after `update` or `update_training`

"""
macro test_criteria(criteria)
    quote
        @test $criteria isa StoppingCriterion
        @testset "state is Nothing" begin
            loss = rand()
            using EarlyStopping: update, update_training

            # Check update for nothing state
            c1 = update($criteria, loss)
            c2 = update($criteria, loss, nothing)
            @test compare_state(c1, c2)

            # Check update_training for nothing state
            c1 = update_training($criteria, loss)
            c2 = update_training($criteria, loss, nothing)
            @test compare_state(c1, c2)
        end

        # Check that done can be called after `update` or `update_training`
        @testset "done" begin
            loss = rand()
            using EarlyStopping: done, update, update_training
            @test done($criteria, nothing) == false
            @test done($criteria, update($criteria, loss)) isa Bool
            @test done($criteria, update_training($criteria, loss)) isa Bool

            # Training then out-of-sample
            loss2 = rand()
            s = update_training($criteria, loss2, update($criteria, loss))
            @test done($criteria, s) isa Bool

            # Out of sample then
            @test done($criteria, update($criteria, loss2, update_training($criteria, loss))) isa Bool
        end

        # Check that `message` can handle state after `update` or `update_training`
        @testset "message" begin
            loss = rand()
            using EarlyStopping: update, update_training
            @test message($criteria, update($criteria, loss)) isa String
            @test message($criteria, update_training($criteria, loss)) isa String
        end
    end
end

"""
    compare_state(s1, s2)

Helper method to check if StoppingCriterion states `s1` and `s2` are effectively
equivalent.
"""
compare_state(s1, s2) = s1 == s2
compare_state(s1::DateTime, s2::DateTime) = abs(s1 - s2) < Millisecond(10)

# Extend == for PQState
import Base: ==
function ==(s1::EarlyStopping.PQState, s2::EarlyStopping.PQState)
    length(s1.training_losses) == length(s2.training_losses) &&
    all(s1.training_losses .== s2.training_losses) &&
    s1.waiting_for_out_of_sample == s2.waiting_for_out_of_sample &&
    s1.loss == s2.loss &&
    s1.min_loss == s2.min_loss
end


text:
mutable struct CachedFactor{T<:Unsigned} <: AbstractFactor{T}
    cache::Dict{Tuple{Integer, Integer}, AbstractVector{T}}
    basefactor::AbstractFactor{T}
    lockobj::Threads.TatasLock
end

Base.length(var::CachedFactor) = length(var.basefactor)

getname(var::CachedFactor) = getname(var.basefactor)

getlevels(var::CachedFactor) = getlevels(var.basefactor)

function CachedFactor(basefactor::AbstractFactor{T}) where {T<:Unsigned}
    CachedFactor{T}(Dict{Tuple{Integer, Integer}, AbstractVector{T}}(), basefactor, Threads.TatasLock())
end

function cache(basefactor::AbstractFactor{T}) where {T<:Unsigned}
    CachedFactor{T}(Dict{Tuple{Integer, Integer}, AbstractVector{T}}(), basefactor, Threads.TatasLock())
end

function slice(factor::CachedFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    basefactor = factor.basefactor
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if isa(basefactor, CachedFactor) || isa(basefactor, Factor)
        slice(basefactor, fromobs, toobs, slicelength)
    else
        lockobj = factor.lockobj
        lock(lockobj)
        try
            if !((fromobs, toobs) in keys(factor.cache))
                v, _ = tryread(slice(basefactor, fromobs, toobs, toobs - fromobs + 1))
                factor.cache[(fromobs, toobs)] = get(v)
            end
        finally
            unlock(lockobj)
        end
        slice(factor.cache[(fromobs, toobs)], 1, toobs - fromobs + 1, slicelength)
    end
end

function Base.map(factor::CachedFactor{T}, dataframe::AbstractDataFrame) where {T<:Unsigned}
    CachedFactor(map(factor.basefactor, dataframe))
end

function isordinal(factor::CachedFactor{T}) where {T<:Unsigned}
    isordinal(factor.basefactor)
end

ext: jl

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using Mustache

tpl = mt"""
Hello {{name}}
You have just won {{value}} dollars!
{{#in_ca}}
Well, {{taxed_value}} dollars, after taxes.
{{/in_ca}}
"""

d = Dict(
"name" => "Chris",
"value" => 10000,
"taxed_value" => 10000 - (10000 * 0.4),
#"in_ca" => true
)

Ans = render(tpl, d);
print(Ans)

tpl = mt"""

<title>{{:TITLE}}</title>

name	summary

{{#:D}}

{{:names}}	{{:summs}}

{{/:D}}

"""
_names = String[]
_summaries = String[]
for s in sort(map(string, names(Main)))
    v = Symbol(s)
    if isdefined(Main,v)
        push!(_names, s)
        push!(_summaries, summary(eval(v)))
    end
end

using DataFrames
d = DataFrame(names=_names, summs=_summaries)

out = render(tpl, TITLE="A quick table", D=d)
print(out)


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash(assets_folder="hr_assets")
app.index_string = """<!DOCTYPE html>

        {%metas%}
        <title>{%title%}</title>
        {%css%}

        {%app_entry%}

            {%config%}
            {%scripts%}
            {%renderer%}

"""
app.layout = html_div(id="hot-reload-content") do
    html_h3("Hot reload")
end

run_server(app,
        dev_tools_hot_reload=true,
        dev_tools_hot_reload_interval=0.1,
        dev_tools_hot_reload_max_retry=100,
)


text:
# * 5%, ** 1%, *** 0.1%
function default_ascii_estim_decoration(s::String, pval::Float64)
  if pval<0.0
      error("p value needs to be nonnegative.")
  end
  if (pval > 0.1)
      return "$s"
  elseif (pval > 0.05)
      return "$s"
  elseif (pval > 0.01)
      return "$s*"
  elseif (pval > 0.001)
      return "$s**"
  else
      return "$s***"
  end
end
function default_latex_estim_decoration(s::String, pval::Float64)
  if pval<0.0
      error("p value needs to be nonnegative.")
  end
  if (pval > 0.1)
      return "$s"
  elseif (pval > 0.05)
      return "$s"
  elseif (pval > 0.01)
      return "$s\\sym{*}"
  elseif (pval > 0.001)
      return "$s\\sym{**}"
  else
      return "$s\\sym{***}"
  end
end


text:
using Images, IJulia, Plots
gr() # FASTER!

#using PyPlot, PyCall,
# using Plots
#@pyimport matplotlib.cm as cm
#@pyimport matplotlib.animation as anim

function makeMask(r_out, r_in=0.0, AA=0.0; normalize=true)
    @assert(AA>=0.0, "Antialiasing border size cannot be negative!")
    mask = zeros(Float64, 2*r_out+1, 2*r_out+1)
    cen = (r_out+1, r_out+1)
    for x in -r_out:r_out
        for y in -r_out:r_out
            l = sqrt(x^2 + y^2)
            if r_in <= l < r_out + AA/2
                if l < r_out - AA/2
                    mask[x+cen[1], y+cen[2]] = 1.0
                else
                    mask[x+cen[1], y+cen[2]] = (r_out + AA/2 - l)/AA
                end
            end
        end
    end

    return normalize ? mask/sum(mask) : mask
end

"""
    makeInnerMask(r_in[, AA=0.0])

Create a matrix filled with a ring of non-zero floats.
The sum of matrix is normalized to be equal to 1.0.

# Arguments
* `r_in::Int`: the inner radius of the ring. Must be smaller than r_out.
* `AA::Float64`: Thickness of antialiasing. 0.0 means no antialiasing.
"""
function makeInnerMask(r_in, AA=0.0; verbose=false)::Array{Float64,2}
    mask = makeMask(r_in, 0.0, AA)
    if verbose; plt[:figure](); plt[:axis]("off"); plt[:imshow](mask) end
    return mask
end

"""
    makeOuterMask(r_out, r_in[, AA=false])

Create a matrix filled with a ring of non-zero floats.
The sum of matrix is normalized to be equal to 1.0.

# Arguments
* `r_out::Int64`: the outer radius of the ring. Must be greater than r_in.
* `r_in::Int64`: the inner radius of the ring. Must be smaller than r_out.
* `AA`: Thickness of antialiasing. 0.0 means no antialiasing.
"""
function makeOuterMask(r_out, r_in, AA=0.0; verbose=false, subtract=false)::Array{Float64,2}
    mask = nothing
    if subtract
        mask = makeOuterMasksubtracted(r_out, r_in, AA)
    else
        mask = makeMask(r_out, r_in, AA)
    end

    if verbose; plt[:figure](); plt[:axis]("off"); plt[:imshow](mask) end
    return mask
end

"""
In this version, the inner mask is subtracted. Using antialiasing, this will lead to an continuous transition from inner to outer.
For further information, see 'makeOuterMask'.
"""
function makeOuterMasksubtracted(r_out, r_in, AA=0.0; verbose=false)::Array{Float64,2}
    mask = makeMask(r_out, 0.0, AA; normalize=false)

    # make the hole of the ring
    #mask_sub = zeros(Float64, 2*r_out+1, 2*r_out+1)
    mask_in = makeMask(r_in, 0.0, AA, normalize=false)
    println(size(mask_in))
    println(size(mask))
    println(size(mask[r_out-r_in+1:r_out+r_in+1, r_out-r_in+1:r_out+r_in+1]))
    flush(stdout)
    mask[r_out-r_in+1:r_out+r_in+1,r_out-r_in+1:r_out+r_in+1] -= mask_in

    #normalize mask, so that sum(mask) == 1.0
    mask /= sum(mask)

    if verbose; plt[:figure](); plt[:axis]("off"); plt[:imshow](mask) end
    return mask
end

"""
Add a splat to the given grid. Can be a square or circle.
Note: In a stable world, splats should become rings.

# Arguments
* `grid`:
* `height`: grid height
* `width`: grid width
* `r_out`: outer radius of conv mask
* `at`   : position as 2D array. If 'nothing', a random position is chosen instead.
* `val`  : Value of splat. Must be within [0,1]
* `round`: if true, make round splats. Make squares otherwise
"""
function splat!(grid::Array{Float64,2}, r_out, at = nothing, val = nothing; round=true)
    width = size(grid)[1]
    height = size(grid)[2]

    if at == nothing
        x = floor(Int64, rand()*width)+1
        y = floor(Int64, rand()*height)+1
    else
        if (x < 0 || x > width || y < 0 || y > height)
            throw(ArgumentError("Given position is outside of the grid boundaries!"))
        end
        x = at[1]
        y = at[2]
    end

    if val == nothing
        c = max(0.5, rand()) # val ∈ [0.5, 1.0]
    else
        if (val < 0 | val > 1)
            throw(ArgumentError("val must be within [0,1]!"))
        end
        c = val
    end

    for dx in -r_out:r_out
        for dy in -r_out:r_out
            # make a splat around [x, y]
            ix = x+dx
            iy = y+dy
            if 1<=ix<=width && 1<=iy<=height && (!round || round && sqrt(dx^2+dy^2) <= r_out)
                grid[ix,iy] = c
            end
        end
    end

    return grid
end

"""
Return a grid initiated with splats.

# Arguments
* `height`: grid height
* `width`: grid width
* `r_out`: outer radius of conv mask
"""
function initGrid(height, width, r_out)::Array{Float64,2}
	grid = zeros(Float64, height,width)
	for t in 0:((width/r_out)*(height/r_out))
		splat!(grid, r_out)
	end
    return grid
end

function sigma1(x::Float64, a::Float64, alpha::Float64)::Float64
    return 1.0 / ( 1.0 + exp(-(x-a)*4.0/alpha));
end

"""
Corrected rules according to ref. implementation
"""
function sigma2(x::Float64, a::Float64, b::Float64)::Float64
    αn = 0.028
    return sigma1(x,a,αn) * ( 1.0 - sigma1(x,b,αn));
end

function sigmam(x::Float64, y::Float64, m::Float64)::Float64
    αm = 0.147
    return x * ( 1.0 - sigma1(m,0.5,αm)) + y*sigma1(m, 0.5, αm);
end

"""
Return the new filling based of the area of the `outer` and the `inner` ring.
"""
function discrete_step(outer::Float64, inner::Float64)::Float64
    # Static Rule Set 1
    # Tested with euler stepping
    # Result: results in static, non moving pattern for dt < 0.35
    # Result: for dt = 0.5, incremental filling can be observed
    #b1 = 0.278
    #b2 = 0.365
    #d1 = 0.267
    #d2 = 0.445
    b1 = 0.257
    b2 = 0.336
    d1 = 0.365
    d2 = 0.549

    return sigmam(sigma2(outer,b1,b2), sigma2(outer,d1,d2),inner);
end

function discrete_as_euler(outer::Float64, inner::Float64)::Float64
    return 2.0 * discrete_step(outer, inner) - 1.0;
end

function next_step_as_euler(f::Float64, outer::Float64, inner::Float64, dt::Float64)::Float64
    return max(0.0, min(1.0, f + dt * discrete_as_euler(outer,inner)));
end

"""
Create two grids. Either empty or randomized. Randomized version uses splats.
Return a tuple. The first entry is a new grid, the second is an uninitialized, 2d array of the same size (for swapping)

Splats will be of size `r_out`.
"""
function createGrid(width=100, height=100, splat=true; r_out=21)::Tuple{Array{Float64,2}, Array{Float64,2}}
    if splat
        return (initGrid(height, width, r_out), Array{Float64, 2}(undef, width, height))
    else
        return (zeros(Float64, width, height), Array{Float64,2}(undef, width, height))
    end
end

"""
Start the smoothlife simulation. All results are estimates. Smaller time steps and greater masks improve resolution
but increase computation time also.

#Arguments
* `runs`: The number of total iterations (convolutions)
* `grid_width`:
* `grid_height`:
* `r_in::Int64`: Size of the inner ring. Greater radius takes more computation time, but improves precision.
* `r_out::Int64`: Size of outer ring. Greater radius takes more computation time, but improves precision/stability. Must be greater than radius_inner. Use 3*radius_inner for best results.
* `dt`: delta time. How much the time passes between two calculations. The greater, the rougher the estimation.
* `sleepTime`: Number of seconds the cpu waits after one complete convolution
* `showOnConsole`: if true, output grid in a python frame. Updates the same frame (video functionality)
* `verbose`: additional information output (e.g. masks)
"""
function simulate(
        ;
        runs = 10,
        start_grid = nothing,
        grid_width = 512,
        grid_height = 512,
        r_in = 7,
        r_out = 3*r_in,
        dt = 0.05,
        sleepTime=0.1,
        save_all_steps = true,
        showOnConsole=true,
        verbose=false
    )::Dict{}

    if start_grid == nothing
        (curGrid, newGrid) = createGrid(grid_width, grid_height; r_out=r_out)
    else
        (curGrid, newGrid) = (start_grid, zeros(Float64, grid_width, grid_height))
    end

    if save_all_steps == true
        # prepare full array for storage
        all_steps = Array{Float64}(undef, grid_width, grid_height, runs)
    end

    # prevent Atom from capturing the plot
    # that will not work!
    #ENV["PLOTS_USE_ATOM_PLOTPANE"] = "false"

    mask_in = makeInnerMask(r_in, 1.0; verbose=verbose)
    mask_out = makeOuterMask(r_out, r_in, 1.0; subtract=true, verbose=verbose)

    #plt[:axis]("off")
    #im = plt[:imshow](curGrid, cmap=cm.Greys_r, vmin = 0.0, vmax = 1.0)

    for run in 1:runs
        # convolution -> get fillings
        convInner = imfilter(curGrid, mask_in, "circular")  # HACK: fft is faster for masks greater 20x20
        convOuter = imfilter(curGrid, mask_out, "circular")

        # calc new frame
        newGrid = zeros(Float64, grid_width, grid_height) # just for debugging.
        for x in 1:grid_width
            for y in 1:grid_height
                #newGrid[x,y] = smoothStep(curGrid[x,y], convOuter[x,y], convInner[x,y], dt)
                newGrid[x,y] = next_step_as_euler(curGrid[x,y], convOuter[x,y], convInner[x,y], dt)
            end
        end

        # swap
        curGrid = newGrid

        if save_all_steps
            all_steps[:,:,run] = curGrid
        end

        # more colorful version, but may need a gradient to control the fill scale
        grad = ColorGradient(:grays)
        # grad = ColorGradient(RGBA{Float64}[RGBA{Float64}(0.05,0.05,0.05,1.0), RGBA{Float64}(0.95,0.95,0.95,1.0)], [0.0, 1.0])
        heatmap(curGrid,
            c = grad, # seems to be the same as :grays, tbh.
            legend = nothing,
            xticks = nothing,
            yticks = nothing,
            aspect_ratio = 1.0)
        gui()
        sleep(sleepTime)
    end

    return Dict("last" => curGrid, "all" => all_steps);
end

function save_runs_as_video(all_steps, file_name, save_type = ".mp4")
    if !(save_type in (".gif", ".mp4", ".mov"))
        throw(ArgumentError("Invalid save type! Please use either '.gif', '.mp4' or '.mov'!"))
    end

    dim = size(res["all"])
    nruns = dim[3]
    width = dim[1]
    height = dim[2]
    anim = @animate for i in 1:nruns
        heatmap(res["all"][:,:,i],
            c = grad, # seems to be the same as :grays, tbh.
            legend = nothing,
            xticks = nothing,
            yticks = nothing,
            aspect_ratio = 1.0)
    end

    # this function will also produce other file types based on the file ending!
    gif(anim, file_name * save_type, fps = 15);
    return anim;
end

max_stars_repo_name: Xentrics/SmoothLifeJulia

max_issues_repo_name: Xentrics/SmoothLifeJulia

max_forks_repo_name: Xentrics/SmoothLifeJulia


text:
function holm(pValues)
  @check_pValues pValues
  n = length(pValues)
  if n <= 1
    return pValues
  end
  sortedIndexes = order(pValues)
  originalOrder = order(sortedIndexes)
  sortedPValues = pValues[sortedIndexes]
  @step_down_adjustment (n + 1 - i) sortedPValues
  min(1, sortedPValues[originalOrder])
end


text:
export RadiationModel, NoRadiation

abstract type RadiationModel end

vars_state(::RadiationModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::RadiationModel,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function integral_set_auxiliary_state!(::RadiationModel, integ::Vars, aux::Vars) end
function integral_load_auxiliary_state!(
    ::RadiationModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) end
function reverse_integral_set_auxiliary_state!(
    ::RadiationModel,
    integ::Vars,
    aux::Vars,
) end
function reverse_integral_load_auxiliary_state!(
    ::RadiationModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) end
function flux_radiation!(
    ::RadiationModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end

struct NoRadiation <: RadiationModel end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
@testset "regression" begin
    Random.seed!(14)
    beta = [10; 20; 30]
    n = 1000
    X = rand(n, 3)
    y = X*beta + randn(n)

    model = regression(y, X)

    ss = statespace(model)

    @test ss.smoother.alpha[end, :] ≈ [10; 20; 30] rtol = 1e-2
    @test ss.model.H[1] ≈ 1.0 rtol = 1e-2
end


text:
## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want

    using StatGeochem
    using Plots; gr();
    using Statistics, DelimitedFiles, SpecialFunctions

## --- Input dataset

    # We'll store the data in a data structure called a "dictionary"
    test = Dict()

    # Make fields called "Latitude" and "Longitude" with locations
    test["Latitude"] = [-26.18, -21.55, 36.6, 54.4, 36.59, 32.5, 49.85, 49.58, 50.2, 42.4725, 48.47, 25.7725, 12.1851, 17.9, -18.1378, 16.1, 67.11, 48.25, 24.42, 23.2, 40.78, -15.38, 21.94, 61.2296, 49.0, 57.03, 63.6, 62.0, 69.47, 54.8189, 22.39, -17.9317, 48.12, 46.47, 49.0, 13.5, 13.5, 33.0, 20.536, 41.5783, 46.5045, 44.5633, 37.7, 43.7506, 37.2202, 58.07, 34.75, 38.25, 40.42, 43.96, 45.38, 37.78, 35.62, 65.25, 46.23, 53.68, 46.38, 45.11, 42.25, 60.25,46.2, 43.31, 39.81, 35.36, 34.88, 34.75, 62.76, 46.45, 54.75, 44.12, 33.61, 38.0, 36.62, 41.11, 47.95, 37.55, 36.5, 36.75, 48.1732, 22.045, 38.45, 61.4597, 36.0, 66.0, 37.5, 23.544, 69.45, 70.0, 6.25, 52.4, 4.58, -32.0, 64.45, 63.89, 63.95, 64.45, 57.03, 8.11, 7.35, 5.14]
    test["Longitude"] = [27.84, 119.9, -118.3, -67.0, 27.17, 103.0, 7.88, 7.15, -60.6, 11.9731, -81.4, 100.32, -84.1862, -65.75, -69.1483, 76.6, 28.85, -78.25, 33.43, 35.1, 14.05, 167.83, 100.92, -131.514, -85.5, 120.72, 36.3, -94.5, 30.6, -100.804, -3.76, -177.567, -77.77, 89.58, -78.0, 78.0, 78.0, 102.0, -104.7, -121.658, -114.802, -114.282, -119.3, -114.641, -119.35, -135.63, -118.75, -91.41, -75.45, -123.66, -109.77, -107.66, -115.16, -147.0, -122.25,-166.4, -122.06, -109.9, -72.12, -152.25, -122.19, -117.31, -120.0, -112.7, -118.88, -118.86, -150.93, -122.03, -131.99, -121.75, -113.09, -72.87, -116.5, -106.25, -91.95, -90.7, -106.12, -106.25, 99.8333, -160.223, 27.3, -139.599, -78.9, -30.0, 140.5, 141.495, 86.22, -52.0, 10.5, -92.75, 9.66, 147.0, 29.05, 28.95, 29.08, 29.05, 120.72, 38.37, 38.42, 123.67]

    # Fill in some other fields with age [Ma], major element [wt. %], and trace elements [ppm] data
    # Notice there are some NaNs, but this is OK. This is a small sampling of real data from EarthChem
    test["Age"] = [3210.0, 3300.0, 85.0, 2100.0, 0.0, 218.0, 275.0, 275.0, 1500.0, 0.254, 2650.0, 44.25, 19.5, 105.8, 1.305, 2600.0, 2415.0, 2100.0, 90.0, 975.0, 32.7685, 0.1, 249.0, 401.5, 2750.0, 3250.0, 2750.0, 2700.0, 2050.0, 1900.0, 188.8, 19.6, 1400.0, 407.0, 2720.0, 2520.0, 2520.0, 748.0, 2.0, 0.149, 46.0, 49.0, 123.0, 46.0, 114.0, 385.5, 33.4, 1400.0, 175.0, 33.4, 2700.0, 33.4, 2100.0, 395.5, 0.005, 33.4, 19.5, 2700.0, 395.5, 175.0, 0.005, 14.55, 33.4, 12.8, 14.55, 156.5, 44.4, 19.5, 429.5, 0.905, 33.4, 397.0, 14.55, NaN, 2750.0, 1400.0, 1400.0, 2100.0, 12.8, 0.5, 218.3, 226.3, 173.55, 44.25, 22.0, 12.8, 158.25, 34.05, 5.0e-6, 2750.0, 44.85, 47.5, 2750.0,2750.0, 2750.0, 2750.0, 3250.0, 1.31, 1.31, 44.85]
    test["SiO2"] = [73.94, 68.07, 67.53, 58.4, 58.8485, 72.53, 72.31, 76.65, 64.37, 50.82, 59.7, 50.87, 55.55, 51.71, 66.01, 50.8, 50.51, 72.11, 58.1, 52.3, 56.69, 50.78, 55.05, 53.2, 58.9, 50.3, 77.9, 57.2, 61.03, 50.68, 54.48, 56.83, 50.34, 50.33, 69.84, 69.57, 73.2, 68.2, 50.26, 52.1, 53.06, 54.01, 69.3, 72.42, 61.47, 50.6, 51.1, 51.6, 52.2, 54.6, 55.3, 57.1, 57.2, 59.9, 60.7, 60.9, 61.25, 61.7, 62.1, 62.5, 63.1, 63.8, 64.9, 66.9, 67.5, 67.8, 68.0, 68.6, 68.7, 68.9, 69.7, 71.3, 73.1, 73.2, 74.3, 75.8, 76.0, 76.9, 47.53, 45.81, 46.21, 48.99, 46.9, 43.86, 48.73, 48.56, 48.3, 45.6, 49.25, 45.74, 47.29, 44.93, 41.7, 44.9, 47.8, 48.2, 49.6, 48.1, 46.35, 49.47]
    test["TiO2"] = [0.28, 0.32, 0.517, 0.72, 1.10718, 0.29, 0.18, 0.07, 0.59, 0.79, 1.31, 0.53, 0.58, 0.7, 0.657, 0.98, 0.93, 0.66, 0.46, 1.3, 0.39, 1.19, 1.19, 0.29, 0.29, 0.83, 0.235, 0.78, 1.29, 0.9, 1.4, 1.2, 0.46, 2.17, 0.34, 0.42,0.2, 0.3, 1.6, 0.93, 1.88, 0.83, 0.48, 0.219, 0.89, 1.65, 1.8, 0.79, 1.16, 2.2, 0.7, 0.84, 1.63, 0.53, 0.7, 0.42, 1.07, 5.00543e-5, 1.6, 0.57, 0.67, 0.25, 0.14, 0.31, 0.34, 0.54, 0.54, 0.54, 0.41, 0.36, 0.37, 0.26, 0.2, 0.25, 0.04, 0.15, 0.18, 0.08, 2.23, 1.31, 3.13, 0.64, 0.55, 3.72, 0.6, 0.77, 1.33, 1.35, 2.93, 0.6, 2.71, 2.2, 0.59, 0.24, 0.48, 0.58, 0.94, 2.37, 2.55, 1.22]
    test["Al2O3"] = [14.19, 16.55, 15.61, 18.5, 17.3459, 14.15, 14.23, 12.6, 14.64, 19.29, 10.42, 9.89, 15.32, 17.14, 16.33, 12.37, 15.94, 14.83, 16.99, 16.41, 18.83, 15.57, 16.85, 13.1, 18.7, 8.48, 12.4, 14.7, 9.2, 15.41, 14.02, 14.16, 18.04, 16.84, 15.11, 15.49, 13.75, 16.83, 13.26, 18.8, 15.82, 10.08, 14.9, 13.27, 16.29, 14.2, 17.1, 13.0, 14.1, 12.3, 14.2, 16.0, 14.2, 15.1, 17.8, 17.4, 15.64, 17.9, 15.4, 15.7, 17.8, 11.8, 18.2, 15.0, 14.7, 15.7, 16.2, 13.8, 13.8, 14.7, 15.2, 15.1, 13.1, 13.3, 13.5, 11.1, 12.1, 13.1, 14.8, 15.13, 14.54, 15.48, 17.8, 9.32, 14.67, 16.29, 16.42, 13.5, 17.88, 11.58, 15.57, 13.53, 4.0, 6.82, 11.3, 14.0, 15.4, 15.17, 14.88, 18.13]
    test["FeOT"] = [2.12355, 2.69043, 3.36529, 6.67659, 6.40665, 1.64665, 1.43683, 0.82, 4.58, 6.87, 6.25366, 7.08, 6.45455, 8.22, 3.63, 11.19, 10.71, 0.57, 8.77, 8.37855, 3.42, NaN, 11.2868, 8.1, 7.56, 12.5, 2.37, 7.74, 15.394, 11.14, 11.31, 9.77, 10.66, 10.55, 2.7, 3.03, 1.61, 2.36, 7.16506, 7.66, 8.82, 8.7, 3.18979, 1.55, NaN, 11.0677, 4.02967, 12.6873, 10.2578, 13.2096, 7.90033, 6.80257, 10.1679, 4.32809, 4.82298, 4.79599, 6.10971, 5.37187, 4.55973, 6.10971, 4.56204, 2.13255, 1.06178, 2.3755, 3.32968, 3.17894, 3.62087, 4.41807, 3.00537, 2.3575, 1.74563,2.59693, 1.1498, 2.77977, 0.584877, 2.45648, 1.51168, 1.36091, 9.84955, 12.28, 12.07, 7.5, 9.21966, 14.1458,NaN, 8.98, 12.02, 10.89, 10.03, 10.84, 11.29, 11.75, 11.4, 9.71, 10.5, 9.28, 13.0, 11.27, 11.56, 6.46]
    test["MgO"] = [0.66, 0.62, 1.44, 3.48, 2.71457, 0.58, 0.71, 0.55, 1.99, 4.32, 1.15, 16.22, 5.75, 5.53, 1.35, 9.97, 6.64, 0.22, 0.44, 4.62, 0.39, 6.35, 3.68, 11.01, 3.4, 13.3, 0.34, 6.1, 1.83, 7.99, 5.36, 2.09, 10.91, 5.62, 1.74, 0.88, 0.16, 1.42, 9.79, 7.16, 7.3, 10.63, 1.1, 0.31, 3.29, 0.98, 2.4, 10.3, 7.51, 2.4, 5.55, 3.11, 3.25, 1.11, 2.32, 2.1, 3.3, 3.12, 1.7, 2.77, 2.2, 0.18, NaN, 1.2, 0.5, 1.6, 1.3, 0.59, 1.11, 0.53, 0.32, 0.55, 0.16, 0.6,0.19, NaN, 0.05, 0.19, 8.05, 11.09, 7.13, 8.82, 9.9, 8.86, 12.71, 8.76, 6.83, 9.57, 3.66, 4.36, 7.13, 11.44,21.5, 22.0, 15.3, 5.45, 7.3, 8.26, 7.57, 7.42]
    test["CaO"] = [0.58, 2.64, 3.65, 3.41, 6.51955, 2.1, 0.34, 0.35, 4.78, 8.6, 6.84, 6.11, 8.68, 7.33, 3.57, 9.88, 9.63, 1.97, 2.97, 7.45, 1.94, 10.56, 0.92, 8.12, 4.3, 13.4, 0.9, 10.36, 2.49, 10.52, 9.74, 8.68, 6.46, 6.56, 3.28, 3.45, 1.21, 2.66, 7.57, 9.23, 7.37, 9.78, 2.8, 3.38, 6.07, 6.95, 10.3, 0.43, 11.0, 5.9, 7.51, 4.2, 6.05, 3.29, 5.57, 4.79, 5.64, 1.73, 4.4, 6.25, 5.45, 1.21, 0.49, 2.9, 2.1, 3.3, 3.7, 2.3, 1.71, 1.83, 1.54, 2.2, 2.0, 2.3, 0.83, 0.03, 0.6, 0.95, 6.59, 10.66, 5.61, 7.69, 11.0, 12.29, 10.2, 11.46, 11.02, 11.1, 9.21, 10.47, 9.23, 10.16, 8.63, 7.57, 9.3, 9.7, 11.5, 9.78, 10.94, 12.45]
    test["Na2O"] = [4.15, 4.89, 3.77, 3.8, 4.15415, 3.22, 2.94, 1.59, 4.31, 1.96, 3.81, 1.46, 2.63, 4.24, 4.17, 1.98, 3.44, 5.67, 6.34, 4.31, 5.69, 3.03, 4.92, 3.3, 2.8, 0.6, 3.36, 1.5, 1.52, 2.95, 2.4, 2.87, 0.05, 4.25, 5.24, 4.48, 3.42, 5.07, 2.98, 3.42, 3.49, 2.57, 3.8, 4.5, 3.23, 3.38, 3.7, NaN, 2.21, 3.4, 3.24, 3.83, 2.49, 2.71, 4.32, 3.85, 4.32, 4.35, 3.0, 2.75, 4.5, 3.92, 3.75, 3.5, 3.7, 3.2, 3.9, 4.2, 3.51, 4.45, 4.9, 3.6, 3.4, 4.0, 3.48, 1.34, 3.1, 3.8, 4.65, 2.52, 3.64, 2.88, 2.02, 1.15, 1.88, 2.14, 2.17, 1.58, 3.79, 2.35, 3.28, 3.46, 0.02, 0.12, 1.42, 2.55, 1.81, 2.98, 2.6, 3.16]
    test["K2O"] = [4.02, 3.1, 3.63, 4.19, 1.81769, 4.48, 5.36, 5.02, 1.1, 6.06, 0.47, 4.95, 1.5, 0.4, 3.82, 0.26, 0.69, 1.55, 3.48, 1.02, 7.33, 1.23, 0.81, 0.76, 3.1, 0.22, 2.04, 0.44, 2.76, 0.13, 0.91, 0.71, 1.82, 1.57, 1.0, 1.25, 5.23, 1.92, 4.66, 0.561, 1.67, 2.91, 3.3, 4.2, 2.6, 2.23, 1.9, 3.21, 0.61, 0.65, 1.39, 3.85, 2.45, 10.3, 1.23, 2.0, 1.79, 4.98, 3.3, 1.59, 1.27, 2.75, 8.98, 2.59, 3.8, 3.2, 2.7, 2.81, 4.58, 2.77, 4.9, 3.1, 4.5, 1.3, 5.27,7.54, 5.1, 2.7, 4.26, 0.38, 0.7, 1.63, 0.24, 0.7, 0.22, 1.09, 0.23, 0.07, 1.86, 0.45, 1.06, 1.48, 0.53, 0.02, 0.05, 0.19, 0.13, 0.97, 0.49, 0.16]
    test["La"] = [15.78, 46.0, 28.8, 35.6, 27.2, 30.18, 50.0, 49.0, 31.1, 91.5, 9.87, 21.4, 9.17, 8.4, 28.6113, 4.34, 8.0, 45.52, 80.7, NaN, 108.0, 15.2, 7.67, 0.8, 36.3, NaN, 46.7, 21.0, 34.2, 2.1, 14.22, 8.81, 1.34, 25.8, 9.0, NaN, 85.1, 7.0, 48.1, 8.76, 36.0, 28.6, 30.0, 44.0, 24.1, 44.7, 50.0, 12.0, 10.7, 50.0, NaN, 43.0, 75.0, 91.0, 13.0, 17.5, 27.0, NaN, 122.0, 8.45, 9.0, 73.0, 53.5, 27.5, 50.0, 50.0, 93.0, 32.0, 28.9, 23.0, 83.0, 13.0, 100.0, 22.6, NaN, 20.2, 200.0, NaN, 67.3, 8.16, 0.0, 2.16, 4.8, 17.707, 3.06, 31.5, 5.54, 2.95, 61.05, 2.18, 31.35, 50.0, 10.0, 20.0, 0.0, 20.0, 3.54, 25.7, 18.39, 2.85]
    test["Yb"] = [NaN, 2.4, NaN, 0.72, 3.1, 0.69, NaN, NaN, 3.28, 2.61, 6.19, 1.46, 1.48, 1.64, 1.47864, 2.2, NaN, 2.71, 5.72, NaN, 4.7, 2.84, 2.81, 1.3, 1.76, NaN, 5.66, 1.71, 5.65, 2.03, 2.48, 5.5, 1.78, 4.1, NaN, NaN, 0.94, 0.36, 1.61, 1.3, 1.65, 1.86, 2.5, 2.49, 2.6, 4.1, 5.0, 2.7, 2.2, 7.0, NaN, 2.6, 5.0, 3.16, 2.0, 1.48, 2.0, 1.5, 5.8,2.74, 1.0, 9.36, 3.0, 1.12, 1.0, 1.5, 2.2, 4.0, 1.64, 3.0, 3.83, 0.77, 3.0, 4.81, NaN, 8.34, 15.0, NaN, 1.36, 1.59, NaN, 1.37, 2.3, 1.435, 1.76, 2.17, 2.65, 1.96, 2.53, 1.31, 1.88, 1.98, NaN, NaN, NaN, NaN, 2.1, 2.18,2.15, 2.26]

    # We're also going to want a list called "elements" with names of all the fields we want to resample
    elements = ["Latitude","Longitude","Age","SiO2","TiO2","Al2O3","FeOT","MgO","CaO","Na2O","K2O","La","Yb"]

    # Now let's add uncertainties, starting with age uncertainty
    test["Age_sigma"] = [10.0, 10.0, 10.0, 400.0, 10.0, 10.0, 24.0, 24.0, 10.0, 0.027, 150.0, 21.25, 3.5, 6.2, 1.295, 10.0, 35.0, 400.0, 10.0, 125.0, 32.7315, 10.0, 10.0, 42.5, 250.0, 250.0, 250.0, 10.0, 450.0, 10.0, 12.8, 14.3, 300.0, 9.0, 10.0, 10.0, 10.0, 7.0, 1.0, 0.26, 4.0, 1.0, 23.0, 3.7, 2.0, 31.5, 31.6, 300.0, 31.0, 31.6, 150.0, 31.6, 400.0, 147.5, 0.005, 31.6, 14.2, 150.0, 147.5, 31.0, 0.005, 9.25, 31.6, 11.0, 9.25, 91.5, 20.6, 14.2, 13.5, 0.895,31.6, 146.0, 9.25, 10.0, 250.0, 300.0, 300.0, 400.0, 10.2, 0.5, 16.7, 24.7, 28.05, 21.25, 10.0, 10.2, 92.75,31.45, 10.0, 250.0, 10.95, 2.5, 250.0, 250.0, 250.0, 250.0, 250.0, 1.295, 1.295, 10.95]

    # For this dataset, lat and lon are good to 0.01 degrees (absolute)
    test["Latitude_sigma"] = 0.01 * ones(size(test["Latitude"]))
    test["Longitude_sigma"] = 0.01 * ones(size(test["Longitude"]))

    # We'll use a 1% relative (1-sigma) default analytical uncertainty for the rest of the elements
    for i=4:length(elements)
        test[elements[i]*"_sigma"] = test[elements[i]] * 0.01
    end

## --- Resample

    # Compute proximity coefficients (inverse weights)
    k = invweight(test["Latitude"], test["Longitude"], test["Age"])

    # # Alternatively, we could weight only by location or only by age (though won't make much difference with this dataset)
    # k = invweight_location(test["Latitude"], test["Longitude"])
    # k = invweight_age(test["Age"])

    # Probability of keeping a given data point when sampling
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the data in each resampling

    # Resample a few hundred times (all elements!)
    nresamplings = 200
    mctest = bsresample(test, nresamplings*length(test["SiO2"]), elements, p)

## --- Approach 1: use the bulk-resampled dataset we just created

    # Calculate mean MgO for 8 bins between 40% SiO2 and 80% SiO2 from resampled dataset
    # (c = bin centers, m = means, e = 1-sigma S.E.M)
    (c,m,e) = binmeans(mctest["SiO2"],mctest["MgO"],40,80,8; resamplingratio=nresamplings)

    # Plot results
    plot(c,m,yerror=2*e,label="",xlabel="SiO2", ylabel="MgO",xlims=(40,80),framestyle=:box)

## --- Approach 2: resample the binned means for one element at a time (Can resample many times)

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_means(test["SiO2"],test["MgO"],40,80,8, p=p, x_sigma=test["SiO2_sigma"], nresamplings=10000)

    # Plot results
    plot(c,m,yerror=(el,eu),label="",xlabel="SiO2", ylabel="MgO",xlims=(40,80), framestyle=:box)

## --- Download and unzip Keller and Schoene (2012) dataset
    if ~isfile("ign.h5") # Unless it already exists
        download("https://storage.googleapis.com/statgeochem/ign.h5.gz","./ign.h5.gz")
        run(`gunzip -f ign.h5.gz`) # Unzip file
    end

    # Read HDF5 file
    using HDF5
    ign = h5read("ign.h5","vars")

## --- Compute proximity coefficients (inverse weights)

    # # Compute inverse weights
    # k = invweight(ign["Latitude"] .|> Float32, ign["Longitude"] .|> Float32, ign["Age"] .|> Float32)

    # Since this is somewhat computatually intensive, let's load a precomputed version instead
    k = ign["k"]

    # Probability of keeping a given data point when sampling
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the data in each resampling
    p[vec(ign["Elevation"].<-100)] .= 0 # Consider only continental crust

    # Age uncertainty
    ign["Age_sigma"] = (ign["Age_Max"]-ign["Age_Min"])/2;
    t = (ign["Age_sigma"] .< 50) .| isnan.(ign["Age_sigma"]) # Find points with < 50 Ma absolute uncertainty
    ign["Age_sigma"][t] .= 50 # Set 50 Ma minimum age uncertainty (1-sigma)

    # Location uncertainty
    ign["Latitude_sigma"] = ign["Loc_Prec"]
    ign["Longitude_sigma"] = ign["Loc_Prec"]

## --- Try resampling a single variable to reproduce the MgO trend from K&S 2012
    xmin = 0
    xmax = 3900
    nbins = 39
    elem = "SiO2"
    unit = "wt. %"

    # Look only at samples in the basaltic silica range
    # (note that if uncertainty in SiO2 were more significant, we should be resampling this too)
    t = 43 .< ign["SiO2"] .< 51 # Mafic

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_means(ign["Age"][t],ign[elem][t],xmin,xmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t])

    # Plot results
    plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkred,mscolor=:darkred,label="")
    plot!(xlabel="Age (Ma)", ylabel="$elem ($unit)",xlims=(0,4000),framestyle=:box,grid=:off,xflip=true) # Format plot

## --- Same as above, but for Na2O
    xmin = 0
    xmax = 3900
    nbins = 39
    elem = "Na2O"
    unit = "wt. %"

    # Look only at samples in the basaltic silica range
    # (note that if uncertainty in SiO2 were more significant, we should be resampling this too)
    t = 43 .< ign["SiO2"] .< 51 # Mafic

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_means(ign["Age"][t],ign[elem][t],xmin,xmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t])

    # Plot results
    plot(c,m,yerror=(el,eu),seriestype=:scatter,markerstrokecolor=:auto,label="")
    plot!(xlabel="Age (Ma)", ylabel="$elem ($unit)",xlims=(0,4000),framestyle=:box,grid=:off,xflip=true) # Format plot


text:
noisify_fixed(df, param) = DataFrame(
    id = df.id,
    periode = df.periode,
    pecahan = df.pecahan,
    value = [d + (rand([-1,1]) * param) for d in df.value]
)

noisify_varied(df, percent) = DataFrame(
    id = df.id,
    periode = df.periode,
    pecahan = df.pecahan,
    value = [d + (rand([-1,1]) * ceil(percent * abs(d))) for d in df.value]
)


text:
module QCCDevDes_Types
export QCCDevDescription, ZoneInfoDesc, GateZoneDesc, JunctionInfoDesc, JunctionDesc
export AuxZoneDesc, LoadZoneInfoDesc, LoadZoneDesc, OperationTimes, setOperationTimes

using StructTypes

OperationTimes = Dict{Symbol, Int64}()
"""Setter for OperationTimes"""
setOperationTimes(times :: Dict{Symbol, Int64}) = global OperationTimes = times

struct ZoneInfoDesc
    id:: String
    end0:: String
    end1:: String
    capacity:: Int64
end

struct GateZoneDesc
    gateZones:: Array{ZoneInfoDesc}
end

struct JunctionInfoDesc
    id:: String
    type:: String
end
struct JunctionDesc
    junctions :: Array{JunctionInfoDesc}
end

struct AuxZoneDesc
    auxZones:: Array{ZoneInfoDesc}
end

struct LoadZoneInfoDesc
    id:: String
    end0:: String
    end1:: String
end
struct LoadZoneDesc
    loadZones:: Array{LoadZoneInfoDesc}
end

struct QCCDevDescription
    gateZone:: GateZoneDesc
    auxZone:: AuxZoneDesc
    junction:: JunctionDesc
    loadZone:: LoadZoneDesc

end

StructTypes.StructType(::Type{ZoneInfoDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{GateZoneDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{JunctionInfoDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{JunctionDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{AuxZoneDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{LoadZoneInfoDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{LoadZoneDesc})= StructTypes.Struct()
StructTypes.StructType(::Type{QCCDevDescription})= StructTypes.Struct()

end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
using Test
using ObservablePmap
using Distributed
using Observables
using Logging

addprocs(2)
#addprocs(2; exeflags="--project")

@everywhere using ObservablePmap

_summ( els, statuses, messages, ids ) = copy(messages)
@testset "schedule_now=$schedule_now, worker_arg=$worker_arg" for
                    schedule_now in (false,true), worker_arg in (false,true)
    function test_pmap_output(summ, task, v, r)
        all_summaries = []
        on(summ) do x
            push!(all_summaries, x)
        end
        if schedule_now
            i = 0
            while !istaskstarted(task) && i < 10
                #println("Waiting for task...")
                sleep(0.1)
                i += 1
            end
            @test istaskstarted(task)
        else
            @test !istaskstarted(task)
            schedule(task)
        end
        wait(task)

        a,b,c = eachrow(reduce(hcat, all_summaries))
        ua,ub,uc = unique.((a,b,c))
        ra,rb,rc = r

        @test all(!isnothing, match.(ra,ua))
        @test all(!isnothing, match.(rb,ub))
        @test all(!isnothing, match.(rc,uc))
        for (x,ux) in ((a,ua),(b,ub),(c,uc))
            @test all(in((0,1)), diff(indexin(x,ux)))
        end
    end

    v = ["a", "b", "c"]
    args = worker_arg ? (default_worker_pool(), v) : (v,)
    summ, task = opmap(args...; schedule_now=schedule_now,
                    summarize=_summ, on_error=identity) do setmessage, x
        for i in 1:3
            setmessage("$x $i")
        end
    end

    test_pmap_output(summ, task, v,
            [[r"^$", r"^a 1$", r"^a 2$", r"^a 3$"],
             [r"^$", r"^b 1$", r"^b 2$", r"^b 3$"],
             [r"^$", r"^c 1$", r"^c 2$", r"^c 3$"]],
    )

    summ, task = ologpmap(args...; schedule_now=schedule_now,
                    summarize=_summ, on_error=identity) do x
        for i in 1:3
            @info "$x $i"
        end
    end
    test_pmap_output(summ, task, v,
            [[r"^$", r"┌ Info: a 1\n└ @ Main.*$",
                     r"┌ Info: a 2\n└ @ Main.*$",
                     r"┌ Info: a 3\n└ @ Main.*$"],
             [r"^$", r"┌ Info: b 1\n└ @ Main.*$",
                     r"┌ Info: b 2\n└ @ Main.*$",
                     r"┌ Info: b 3\n└ @ Main.*$"],

             [r"^$", r"┌ Info: c 1\n└ @ Main.*$",
                     r"┌ Info: c 2\n└ @ Main.*$",
                     r"┌ Info: c 3\n└ @ Main.*$"]]
    )

end


text:
module Actions
    import ..Wilsonloops:Wilson_loop_set,make_staples,Wilson_loop_set,
            make_cloverloops,Tensor_derivative_set, make_loops
    import ..SUN_generator:Generator
    import ..Rhmc:RHMC

    abstract type GaugeActionParam end

    abstract type FermiActionParam end

    abstract type SmearingParam end

    struct SmearingParam_nosmearing <: SmearingParam
    end

    abstract type SmearingParam_single <: SmearingParam
    end

    abstract type SmearingParam_multi <: SmearingParam
    end

    mutable struct SmearingParam_stout <: SmearingParam_single
        staples_for_stout::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative::Array{Tensor_derivative_set,1}
        staples_for_stout_dag::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative_dag::Array{Tensor_derivative_set,1}
        ρs::Array{Float64,1}
        #ρs::Array{Float64,1}
    end

    mutable struct SmearingParam_stout_multi <: SmearingParam_multi
        staples_for_stout::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative::Array{Tensor_derivative_set,1}
        staples_for_stout_dag::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative_dag::Array{Tensor_derivative_set,1}
        ρs::Array{Array{Float64,1},1}
        #ρs::Array{Float64,1}
    end

    const Nosmearing = SmearingParam_nosmearing
    const Stout = SmearingParam_stout

    mutable struct GaugeActionParam_standard <: GaugeActionParam
        β::Float64
        NTRACE::Int64
    end

    mutable struct GaugeActionParam_autogenerator <: GaugeActionParam
        βs::Array{Float64,1}
        β::Float64
        numactions::Int64
        NTRACE::Int64
        loops::Array{Wilson_loop_set,1}
        staples::Array{Array{Wilson_loop_set,1},1}
        couplinglist::Array{String,1}
    end

    """
    ```Setup_Gauge_action(β;actiontype="standard")````
    - β: Coupling value

    Set up the information about the Gauge action.
    You can set the coupling value β.

    Now only SU(3) case is supported.
    """
    function Setup_Gauge_action(β;actiontype="standard")
        if actiontype == "standard"
            NTRACE =3
            return GaugeActionParam_standard(β,NTRACE)
        end
    end

    Base.@kwdef struct FermiActionParam_Wilson <: FermiActionParam
        hop::Float64  = 0.141139#Hopping parameter
        r::Float64  = 1#Wilson term
        #eps::Float64 = 1e-8
        eps::Float64 = 1e-19
        Dirac_operator::String = "Wilson"
        MaxCGstep::Int64 = 3000
        quench::Bool = false
        smearing::SmearingParam

        function FermiActionParam_Wilson(hop,r,eps,Dirac_operator,MaxCGstep,quench
            ;smearingparameters = "nothing",
            loops_list = nothing,
            coefficients  = nothing,
            numlayers = 1,
            L = nothing)

            smearing = construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)
            #smearing = Nosmearing()
            #if smearingparameters == nothing
            #    smearing = Nosmearing()
            #end
            return new(hop,r,eps,Dirac_operator,MaxCGstep,quench,smearing)
        end

    end

    const Clover_coefficient  = 1.5612

    Base.@kwdef struct FermiActionParam_WilsonClover <: FermiActionParam
        hop::Float64  = 0.141139#Hopping parameter
        r::Float64  = 1#Wilson term
        #eps::Float64 = 1e-8
        eps::Float64 = 1e-19
        Dirac_operator::String = "WilsonClover"
        MaxCGstep::Int64 = 3000
        Clover_coefficient::Float64 = Clover_coefficient
        #CloverFμν::Array{ComplexF64,4}
        internal_flags::Array{Bool,1}
        inn_table::Array{Int64,3}
        _ftmp_vectors::Array{Array{ComplexF64,3},1}
        _is1::Array{Int64,1}
        _is2::Array{Int64,1}
        quench::Bool = false
        SUNgenerator::Union{Nothing,Generator}
        _cloverloops::Array{Wilson_loop_set,2}
        smearing::SmearingParam

        function FermiActionParam_WilsonClover(hop,r,eps,Dirac_operator,MaxCGstep,
            Clover_coefficient,internal_flags,inn_table,_ftmp_vectors,
            _is1,_is2,
            quench, SUNgenerator,_cloverloops;smearingparameters = "nothing",
            loops_list = nothing,
            coefficients  = nothing,
            numlayers = 1,
            L = nothing)

            smearing = construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)

            #if smearingparameters == nothing
            #    smearing = Nosmearing()
            #end
            #smearing = Nosmearing()

            return new(hop,r,eps,Dirac_operator,MaxCGstep,
            Clover_coefficient,internal_flags,inn_table,_ftmp_vectors,
            _is1,_is2,
            quench, SUNgenerator,_cloverloops,smearing)
        end
    end

    function construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)
        if smearingparameters == "nothing"
            smearing = Nosmearing()
        else
            @assert loops_list != nothing "loops should be put if you want to use smearing schemes"
            loops = make_loops(loops_list,L)

            @assert coefficients != nothing "coefficients should be put if you want to use smearing schemes"
            println("stout smearing will be used")
            if numlayers == 1
                smearing = SmearingParam_stout(loops,coefficients)
            else
                numloops = length(loops)
                smearing_single = SmearingParam_stout(loops,rand(numloops))
                smearing = SmearingParam_stout_multi(smearing_single.staples_for_stout,
                    smearing_single.tensor_derivative,
                    smearing_single.staples_for_stout_dag,
                    smearing_single.tensor_derivative_dag,
                    coefficients
                    )

            end
            #println(smearing )
            #exit()
        end
        return smearing
    end

    Base.@kwdef struct FermiActionParam_Staggered <: FermiActionParam
        mass::Float64 = 0.5
        eps::Float64 = 1e-19
        Dirac_operator::String = "Staggered"
        MaxCGstep::Int64 = 3000
        quench::Bool = false
        Nf::Int8 = 4
        rhmc_action::Union{Nothing,RHMC}
        rhmc_MD::Union{Nothing,RHMC}
        smearing::SmearingParam

        function FermiActionParam_Staggered(
            mass,
            eps,
            Dirac_operator,
            MaxCGstep,
            quench,
            Nf;smearingparameters = "nothing",
            loops_list = nothing,
            coefficients  = nothing,
            numlayers = 1,
            L = nothing
            ) where T <: Real

            if Nf == 4 || Nf == 8 # 8 flavors if phi (Mdag M)^{-1} phi
                rhmc_action = nothing
                rhmc_MD = nothing
            else
                #for action: r_action
                #Nf = 8 -> alpha = 1 -> power x^{1/2}
                #Nf = 2 -> alpha = 1/4 -> power x^1/8
                #Nf = 1 -> alpha = 1/8  -> power x^1/16
                order = Nf //16

                rhmc_action = RHMC(order,n=15)

                #for MD: r_MD
                #Nf = 8 -> alpha = 1 -> power x^{1}
                #Nf = 2 -> alpha = 1/4 -> power x^1/4
                #Nf = 1 -> alpha = 1/8  -> power x^1/8
                order = Nf // 8
                #rhmcorder = 8 ÷ Nf
                rhmc_MD = RHMC(order,n=10)
            end

            smearing = construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)

            return new(
            mass,
            eps,
            Dirac_operator,
            MaxCGstep,
            quench,
            Nf,
            rhmc_action,
            rhmc_MD,
            smearing
            )
        end
    end

    function GaugeActionParam_autogenerator(βs,loops,NC,couplinglist)
        @assert length(βs) == length(loops) "The number of loops should be the number of βs!"
        numactions = length(loops)
        β = βs[1]
        staples = Array{Array{Wilson_loop_set,1},1}(undef,numactions)
        for (i,loop) in enumerate(loops)
            staples[i] = make_staples(loop)
            #=
            println("$i-th action")
            display(loop)
            println("$i-th action's staple")
            for mu=1:4
                println("$mu -direction")
                display(staples[i][mu])
            end
            =#
        end

        return GaugeActionParam_autogenerator(βs,β,numactions,NC,loops,staples,couplinglist)
    end

    Base.@kwdef struct FermiActionParam_Domainwall <: FermiActionParam
        N5::Int64
        r::Float64 = 1
        M::Float64 = -1
        m::Float64
        ωs::Array{Float64,1}
        b::Float64
        c::Float64
        eps::Float64 = 1e-19
        Dirac_operator::String = "Domainwall"
        MaxCGstep::Int64 = 3000
        quench::Bool = false
        smearing::SmearingParam

        function FermiActionParam_Domainwall(N5,r,M,m,ωs,b,c,eps,Dirac_operator,MaxCGstep,quench
            ;smearingparameters = "nothing",
            loops_list = nothing,
            coefficients  = nothing,
            numlayers = 1,
            L = nothing)

            smearing = construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)
            #smearing = Nosmearing()
            #if smearingparameters == nothing
            #    smearing = Nosmearing()
            #end
            return new(N5,r,M,m,ωs,b,c,eps,Dirac_operator,MaxCGstep,quench,smearing)
        end

        function FermiActionParam_Domainwall(N5,r,M,m,eps,Dirac_operator,MaxCGstep,quench
            ;smearingparameters = "nothing",
            loops_list = nothing,
            coefficients  = nothing,
            numlayers = 1,
            L = nothing)

            b = 1
            c = 1
            ωs = ones(Float64,Ns)

            smearing = construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)
            #smearing = Nosmearing()
            #if smearingparameters == nothing
            #    smearing = Nosmearing()
            #end
            return new(N5,r,M,m,ωs,b,c,eps,Dirac_operator,MaxCGstep,quench,smearing)
        end

    end

    function show_parameters_action(fparam::FermiActionParam_Staggered)
        println("#--------------------------------------------------")
        println("#Fermion Action parameters")
        println("#type: ",typeof(fparam))
        println("mass = ",fparam.mass)
        println("eps = ",fparam.eps)
        println("Dirac_operator = ",fparam.Dirac_operator)
        println("MaxCGstep = ",fparam.MaxCGstep)
        println("quench = ",fparam.quench)
        println("Nf = ",fparam.Nf)
        println("#--------------------------------------------------")
    end

    function show_parameters_action(fparam::FermiActionParam_Wilson)
        println("#--------------------------------------------------")
        println("#Fermion Action parameters")
        println("#type: ",typeof(fparam))
        println("hop = ",fparam.hop)
        println("r = ",fparam.r)
        println("eps = ",fparam.eps)
        println("Dirac_operator = ",fparam.Dirac_operator)
        println("MaxCGstep = ",fparam.MaxCGstep)
        println("quench = ",fparam.quench)
        println("#--------------------------------------------------")
    end

    function show_parameters_action(fparam::FermiActionParam_WilsonClover)
        println("#--------------------------------------------------")
        println("#Fermion Action parameters")
        println("#type: ",typeof(fparam))
        println("hop = ",fparam.hop)
        println("r = ",fparam.r)
        println("eps = ",fparam.eps)
        println("Dirac_operator = ",fparam.Dirac_operator)
        println("MaxCGstep = ",fparam.MaxCGstep)
        println("Clover_coefficient = ",fparam.Clover_coefficient)
        println("quench = ",fparam.quench)
        println("#--------------------------------------------------")
    end

    function show_parameters_action(gparam::GaugeActionParam_standard)
        println("#--------------------------------------------------")
        println("#Gauge Action parameters")
        println("#type: ",typeof(gparam))
        println("β = ",gparam.β)
        println("NC = ",gparam.NTRACE)
        println("#--------------------------------------------------")
    end

    function show_parameters_action(gparam::GaugeActionParam_autogenerator)
        println("#--------------------------------------------------")
        println("#Gauge Action parameters")
        println("#type: ",typeof(gparam))
        println("#Num. of action terms: ",gparam.numactions)
        println("#coupling βs = ",gparam.βs)
        #=
        println("#actions that we consider: ")
        for (i,loop) in enumerate(gparam.loops)
            println("#------------------------------")
            println("#$i-th action term consists of")
            display(loop)
            println("#------------------------------")
        end
        =#
        println("NC = ",gparam.NTRACE)
        println("#--------------------------------------------------")
    end

    function FermiActionParam_WilsonClover(hop,r,eps,MaxCGstep,NV,Clover_coefficient,NC;quench=false)
        #CloverFμν = zeros(ComplexF64,3,3,NV,6)
        inn_table= zeros(Int64,NV,4,2)
        internal_flags = zeros(Bool,2)
        _ftmp_vectors = Array{Array{ComplexF64,3},1}(undef,6)
        _is1 = zeros(Int64,NV)
        _is2 = zeros(Int64,NV)
        for i=1:6
            _ftmp_vectors[i] = zeros(ComplexF64,3,NV,4)
        end
        #return FermiActionParam_WilsonClover(hop,r,eps,"WilsonClover",MaxCGstep,Clover_coefficient,CloverFμν,inn_table,internal_flags,_ftmp_vectors,_is1,_is2,quench)
        if NC ≥ 4
            SUNgenerator = Generator(NC)
        else
            SUNgenerator = nothing
        end

        _cloverloops = Array{Wilson_loop_set,2}(undef,3,4)
        for μ=1:3
            for ν=μ+1:4
                _cloverloops[μ,ν] = make_cloverloops(μ,ν)
            end
        end

        return FermiActionParam_WilsonClover(hop,r,eps,"WilsonClover",MaxCGstep,Clover_coefficient,inn_table,internal_flags,_ftmp_vectors,_is1,_is2,quench,
                        SUNgenerator,_cloverloops)
    end

    function FermiActionParam_Wilson(hop,r,eps,MaxCGstep;quench=false)
        return FermiActionParam_Wilson(hop,r,eps,"Wilson",MaxCGstep,quench)
    end

    function FermiActionParam_Wilson(hop,r,eps,Dirac,MaxCGstep;quench=false)
        return FermiActionParam_Wilson(hop,r,eps,Dirac,MaxCGstep,quench)
    end

    """
    ```Setup_Fermi_action(Dirac_operator= "Wilson")```

    Set up the information about the Fermion action.

    Now only WilsonFermion case is supported.

    # For example
    ```julia
        fparam = Setup_Fermi_action()
    ```

    The default values are

    ```julia
    hop::Float64  = 0.141139
    r::Float64  = 1
    eps::Float64 = 1e-19
    Dirac_operator::String = "Wilson"
    MaxCGstep::Int64 = 3000
    ```

    - hop : hopping parameter
    - r : Wilson term
    - eps : convergence criteria in the CG method
    - MaxCGstep : maximum number of the CG steps

    If you want to change the parameters for the Wilson Fermions,
    please do as follows.

    ```julia
        fparam = FermiActionParam_Wilson(hop,r,eps,MaxCGstep)
    ```

    """
    function Setup_Fermi_action(Dirac_operator= "Wilson")
        if Dirac_operator == "Wilson"
            return FermiActionParam_Wilson()
        end
    end

    function SmearingParam_stout(loops_smearing,ρs)
        num = length(loops_smearing)
        @assert num == length(ρs) "number of ρ elements in stout smearing scheme should be $num. Now $(length(ρs))"
        staplesforsmear_set = Array{Wilson_loop_set,1}[]
        staplesforsmear_dag_set = Array{Wilson_loop_set,1}[]
        println("staple for stout smearing")

        tensor_derivative = Array{Tensor_derivative_set,1}(undef,num)
        tensor_derivative_dag = Array{Tensor_derivative_set,1}(undef,num)

        for i=1:num
            loop_smearing = loops_smearing[i]

            staplesforsmear = Array{Wilson_loop_set,1}(undef,4)
            staplesforsmear_dag = Array{Wilson_loop_set,1}(undef,4)

            staple = make_staples(loop_smearing)

            for μ=1:4
                staplesforsmear_dag[μ] = staple[μ]##make_plaq_staple(μ)
                staplesforsmear[μ] = staplesforsmear_dag[μ]'
                #println("$μ -direction")
                #display(staplesforsmear[μ])
                #staplesforsmear[μ] = make_plaq_staple(μ)
                #staplesforsmear_dag[μ] = staplesforsmear[μ]'
                #println("dagger: $μ -direction")
                #display(staplesforsmear_dag[μ])
            end
            tensor_derivative[i] = Tensor_derivative_set(staplesforsmear)
            tensor_derivative_dag[i] = Tensor_derivative_set(staplesforsmear_dag)

            push!(staplesforsmear_set,staplesforsmear )
            push!(staplesforsmear_dag_set,staplesforsmear_dag)

        end

        return SmearingParam_stout(
            staplesforsmear_set,
            tensor_derivative,
            staplesforsmear_dag_set,
            tensor_derivative_dag,
            ρs
            )
    end

end


text:
using Query
using DataFrames
using Statistics
using Test

@testset "a.b Syntax (dplyr API)" begin

    df = DataFrame(name=repeat(["John", "Sally", "Kirk"],inner=[1],outer=[2]),
                   age=vcat([10., 20., 30.],[10., 20., 30.].+3),
                   children=repeat([3,2,2],inner=[1],outer=[2]),state=[:a,:a,:a,:b,:b,:b])

    x = @from i in df begin
        @group i by i.state into g
        @select {group=key(g),mage=mean(g.age), oldest=maximum(g.age), youngest=minimum(g.age)}
        @collect DataFrame
    end

    @test x isa DataFrame
    @test size(x) == (2,4)
    @test x[1,:mage] == 20
    @test x[2,:mage] == 23
    @test x[1,:oldest] == 30
    @test x[2,:oldest] == 33
    @test x[1,:youngest] == 10
    @test x[2,:youngest] == 13
end


text:
using Documenter, ParametricOptInterface

makedocs(
    modules = [ParametricOptInterface],
    doctest  = false,
    clean = true,
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(assets = ["assets/favicon.ico"],
                    mathengine = Documenter.MathJax(),
                    prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "ParametricOptInterface.jl",
    authors = "Tomás Gutierrez, and contributors",
    pages = [
        "Home" => "index.md",
        "manual.md",
        "example.md"
        # "reference.md"
    ]
)

deploydocs(
    repo = "github.com/jump-dev/ParametricOptInterface.jl.git",
)


text:
mutable struct OptimizerList <: Optimizer
    opts::Array{Optimizer}
end

"Promote Optimizer to OptimizerList"
function update(opt::Optimizer,
                g_t::Array{Array{Float64, N}, 1}) where {N}
    # length of list
    n = length(g_t)

    opt = OptimizerList([deepcopy(opt) for k in 1:n])

    update(opt, g_t)
end

function update(opt::OptimizerList,
                g_t::Array{Array{Float64, N}, 1})  where {N}
    # length of list
    n = length(g_t)

    δ = Array{Array{Float64}}(n)

    for i in 1:n
        δ[i] = update(opt.opts[i], g_t[i])
    end

    return δ
end


text:
module TestMonoidalUndirectedWiringDiagrams
using Test

using Catlab.Theories, Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra.CSets: subpart

# Categorical interface
#######################

const HDOb = HypergraphDiagramOb{Symbol,Symbol}
const HDHom = HypergraphDiagramHom{Symbol,Symbol}

# Category
#---------

# Domains and codomains
A, B = HDOb([:A]), HDOb([:B])
f = HDHom(singleton_diagram(A⊗B, name=:f), dom=BitArray([1,0]))
g = HDHom(singleton_diagram(B⊗A, name=:g), codom=BitArray([0,1]))
@test subpart(f.diagram, :name) == [:f]
@test subpart(g.diagram, :name) == [:g]
@test dom(f) == A
@test codom(f) == B
@test dom(g) == B
@test codom(g) == A

# Composition and identities
fg = compose(f,g)
@test dom(fg) == A
@test codom(fg) == A
@test nboxes(fg.diagram) == 2
@test njunctions(fg.diagram) == 3
@test compose(id(dom(f)), f) == f
@test compose(f, id(codom(f))) == f

# Symmetric monoidal category
#----------------------------

# Monoidal products
X, Y = HDOb([:A,:B]), HDOb([:C,:D])
I = munit(HDOb)
@test otimes(A,B) == X
@test otimes(X,I) == X
@test otimes(I,X) == X
@test dom(f⊗g) == A⊗B
@test codom(f⊗g) == B⊗A

# Braiding
@test compose(braid(X,Y),braid(Y,X)) == id(otimes(X,Y))

# Hypergraph category
#--------------------

@test compose(create(X), mcopy(X)) == dunit(X)
@test compose(mmerge(X), delete(X)) == dcounit(X)
# XXX: Not equal due to ordering of outer ports, but isomorphic.
#@test (id(B)⊗dunit(A)) ⋅ (id(B)⊗f⊗id(A)) ⋅ (dcounit(B)⊗id(A)) == dagger(f)

end


text:
include(joinpath("..", "required.jl"))

## Conditional arguments
outcome = "bart"
subset_qc = true
# save_figures = true

# Subset to QC data (optional)
if (@isdefined subset_qc) && subset_qc == true
    @rput subset_qc
end

## Train & apply models (RF or BART)
model_script = joinpath("src", "others", "x_direct-models_$(outcome).R")
@rput model_script
@time begin
    R"""
    source(model_script)
    """
end # about 1 minute
@rget results

## Fix missing values
predictions = results[:predictions]

## Get comparison layers
# Load species names
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name

# Load distributions
raw_distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "raster", "distributions_raw.tif"), i
    ) for i in eachindex(spenames)
]
sdm_distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "raster", "distributions_bart.tif"), i
    ) for i in eachindex(spenames)
]

# Subset to Quebec
if (@isdefined subset_qc) && subset_qc == true
    coords_qc = (left=-80.0, right=-55.0, bottom=45.0, top=63.0)
    raw_distributions = [d[coords_qc] for d in raw_distributions]
    sdm_distributions = [d[coords_qc] for d in sdm_distributions]
end

## Arrange predictions as layers
richness_pred = SimpleSDMResponse(predictions, :richness, raw_distributions[1])
lcbd_pred = SimpleSDMResponse(predictions, :lcbd, raw_distributions[1])
replace!(richness_pred, missing => nothing)
replace!(lcbd_pred, missing => nothing)
pred = (richness=richness_pred, lcbd=lcbd_pred)

## Plot predictions
richness_plot = plot_layer(
    richness_pred;
    c=:cividis,
    title="Direct richness predictions ($(uppercase(outcome)))",
    colorbar_title="Predicted number of species",
    dpi=200,
)

## Plot differences

# Get comparisons richness & lcbd
Yraw = Ymatrix(raw_distributions)
Ysdm = Ymatrix(sdm_distributions)
raw, sdm = [
    (
        richness=richness(Y, raw_distributions[1]),
        lcbd=lcbd(Y, raw_distributions[1]),
    ) for Y in (Yraw, Ysdm)
]

# Custom calculation function
function difference(prediction, comparison; absolute=false)
    # Calculate difference
    diff_layer = prediction - comparison
    # Get absolute values (optional)
    if absolute
        replace!(x -> !isnothing(x) ? abs(x) : x, diff_layer.grid)
    end
    return diff_layer
end

# Calculate differences
diff_pred_raw = difference(pred.richness, raw.richness)
diff_pred_sdm = difference(pred.richness, sdm.richness)
diff_sdm_raw = difference(sdm.richness, raw.richness)
diffs = [diff_pred_raw, diff_pred_sdm, diff_sdm_raw]
map(x -> extrema(filter(!isnothing, x.grid)), diffs)
maxlim = 30

# Custom plot function
function difference_plot(layer, lim; title="")
    limrange = (-lim, lim)
    diff_map = plot_layer(
        layer;
        c=cgrad(:PuOr; rev=true),
        clim=limrange,
        title="Richness difference",
        colorbar_title="Difference",
    )
    diff_hist = histogram(
        [
            filter(x -> !isnothing(x) && x > 0, layer.grid),
            filter(x -> !isnothing(x) && x <= 0, layer.grid),
        ];
        bins=:rice,
        c=cgrad(:PuOr; rev=true),
        legend=:none,
        ylim=limrange, # xlabel = "Difference",
        title="Distribution of difference values",
        orientation=:horizontal,
    )
    diff_title = plot(; annotation=(0.5, 0.5, "$(title)"), framestyle=:none)
    l = @layout [t{0.01h}; a{0.6w} b{0.38w}]
    diff_plot = plot(diff_title, diff_map, diff_hist; size=(800, 400), dpi=200, layout=l)
    return diff_plot
end

# Plot differences
plot_pred_sdm = difference_plot(
    diff_pred_sdm, maxlim; title="Direct richness $(uppercase(outcome)) vs SDM richness"
)
plot_pred_raw = difference_plot(
    diff_pred_raw, maxlim; title="Direct richness $(uppercase(outcome)) vs raw richness"
)
plot_sdm_raw = difference_plot(diff_sdm_raw, maxlim; title="SDM richness vs raw richness")

## Extras

# Absolute richness difference
richness_absdiff = plot_layer(
    difference(pred.richness, sdm.richness; absolute=true);
    c=:inferno, # clim = (-Inf, Inf),
    title="Direct richness predictions difference ($(uppercase(outcome)))",
    colorbar_title="Difference from SDM-predicted richness (absolute)",
    dpi=200,
)

# LCBD predictions
lcbd_plot = plot_layer(
    lcbd_pred;
    c=:viridis,
    title="Direct LCBD predictions ($(uppercase(outcome)))",
    colorbar_title="LCBD scores",
    clim=(0, 1),
    dpi=200,
)

# LCBD difference
lcbd_diff = plot_layer(
    difference(pred.lcbd, sdm.lcbd; absolute=true);
    c=:inferno, # clim = (-Inf, Inf),
    title="Direct richness predictions difference ($(uppercase(outcome)))",
    colorbar_title="Difference in predicted LCBD (absolute)",
    dpi=200,
)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
function one_update_step_hyperparameters!(grammar, sequences)
    forests = map(x->run_chartparser(x, grammar), sequences)
    expected_rule_usages_for_each_sequence =
        map(x->calculate_outside_probs_and_expected_rule_usages!(x, grammar), forests)

    cats = union([keys(d) for d in expected_rule_usages_for_each_sequence]...)
    for cat in cats
        d = sum(d[cat] for d in expected_rule_usages_for_each_sequence if haskey(d, cat))
        set_params!(grammar.rule_cond.dists[dependent_components(cat)], d)
    end
    old_likelihood = prod(score(forest) for forest in forests)
end

function learn_hyperparameters!(grammar, sequences)
    old_likelihood = zero(LogProb)
    old_old_likelihood = zero(LogProb)
    i = 1
    while true
        println("iteration step $i")
        old_likelihood = one_update_step_hyperparameters!(grammar, sequences)
        println("old likelihood $old_likelihood")
        if abs(old_likelihood.value - old_old_likelihood.value) < 0.001
            break
        end
        i += 1
    end
end


text:
function insert_element(values)
  push!(values, -10)
end

list = [1,2,3]
insert_element(list)
println(list)

f(x,y) = x^3 + x*y - y
g(a, b = 5) = a + b
k(x; a1=1, a2=2) = x * (a1+a2)

println(f(3,2))

println("g(3) + f(1,2) = $( g(3) + f(1,2) )")
println("k(2) = $(k(2))")

cubes = map(x->x^3, [1,2,3,5])
println("$cubes")

println("$([1:3])")

a = split("A,B,C,D",",")
println("typeof(a): $(typeof(a))")
println("a: $(a)")

str::String = "This is a string"
println(str)


text:
using Catalyst, Test
using LinearAlgebra

include("test_networks.jl")

rn = @reaction_network begin
    (1, 2), A + B <--> C
    (3, 2), D <--> E
    (0.1, 0.2), E <--> F
    6, F --> G
    7, H --> G
    5, 0 --> K
    3, A + B --> Z + C
end

S = netstoichmat(rn)
C = conservationlaws(S)

@test size(C,1) == 3

b = [ 0, 0, 0, 1, 1, 1, 1, 1, 0, 0 ]
@test any(b == C[i,:] for i in 1:size(C,1))

# For the A + B <--> C subsystem one of these must occur
# as a conservation law
D = [ 1 -1 0 0 0 0 0 0 0 0;
      -1 1 0 0 0 0 0 0 0 0
      0  1 1 0 0 0 0 0 0 0 ]
@test any(D[j,:] == C[i,:] for i in 1:size(C,1), j in 1:size(D,1))

Ss_standard = map(netstoichmat, reaction_networks_standard)
Cs_standard = map(conservationlaws, Ss_standard)
@test all(size(C, 1) == 0 for C in Cs_standard)

Ss_hill = map(netstoichmat, reaction_networks_hill)
Cs_hill = map(conservationlaws, Ss_hill)
@test all(size(C, 1) == 0 for C in Cs_hill)

function consequiv(A, B)
    rank([A; B]) == rank(A) == rank(B)
end

Ss_constraint = map(netstoichmat, reaction_networks_constraint)
Cs_constraint = map(conservationlaws, Ss_constraint)

@test all(consequiv.(Matrix{Int}.(Cs_constraint), reaction_network_constraints))