text:
export visualize

"""
visualize(states::AbstractArray; statenames = string.(1:length(states)), quantiles = (0.1, 0.99), aspect = (1,1,1), resolution = (1920, 1080), statistics = false, title = "Field = ")
# Description
Visualize 3D states
# Arguments
- `states`: Array{Array{Float64,3},1}. An array of arrays containing different fields
# Keyword Arguments
- `statenames`: Array{String,1}. An array of stringnames
- `aspect`: Tuple{Int64,Int64,Float64}. Determines aspect ratio of box for volumes
- `resolution`: Resolution of preliminary makie window
- `statistics`: boolean. toggle for displaying statistics
# Return
- `scene`: Scene. A preliminary scene object for manipulation
"""
function visualize(
    states::AbstractArray;
    statenames = string.(1:length(states)),
    units = ["" for i in eachindex(states)],
    aspect = (1, 1, 1),
    resolution = (1920, 1080),
    statistics = false,
    title = "Field = ",
    bins = 300,
)
    # Create scene
    scene, layout = layoutscene(resolution = resolution)
    lscene = layout[2:4, 2:4] = LScene(scene)
    width = round(Int, resolution[1] / 4) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    if statistics
        llscene =
            layout[4, 1] = Axis(
                scene,
                xlabel = @lift(statenames[$statenode] * units[$statenode]),
                xlabelcolor = :black,
                ylabel = "pdf",
                ylabelcolor = :black,
                xlabelsize = 40,
                ylabelsize = 40,
                xticklabelsize = 25,
                yticklabelsize = 25,
                xtickcolor = :black,
                ytickcolor = :black,
                xticklabelcolor = :black,
                yticklabelcolor = :black,
            )
        layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)
    end

    # x,y,z are for determining the aspect ratio of the box
    if (typeof(aspect) <: Tuple) & (length(aspect) == 3)
        x, y, z = aspect
    else
        x, y, z = size(states[1])
    end

    # Clim sliders
    upperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    clims = @lift((
        quantile($state[:], $lowerclim_node),
        quantile($state[:], $upperclim_node),
    ))
    cmap_rgb = @lift(to_colormap($colornode))
    titlename = @lift(title * $statename) # use padding and appropriate centering

    # Statistics
    if statistics
        histogram_node = @lift(histogram($state, bins = bins))
        xs = @lift($histogram_node[1])
        ys = @lift($histogram_node[2])
        pdf = GLMakie.AbstractPlotting.barplot!(
            llscene,
            xs,
            ys,
            color = :red,
            strokecolor = :red,
            strokewidth = 1,
        )
        @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
        @lift(GLMakie.AbstractPlotting.ylims!(
            llscene,
            extrema($histogram_node[2]),
        ))
        vlines!(
            llscene,
            @lift($clims[1]),
            color = :black,
            linewidth = width / 100,
        )
        vlines!(
            llscene,
            @lift($clims[2]),
            color = :black,
            linewidth = width / 100,
        )
    end

    # Volume Plot
    volume!(
        lscene,
        0..x,
        0..y,
        0..z,
        state,
        camera = cam3d!,
        colormap = cmap_rgb,
        colorrange = clims,
    )
    # Camera
    cam = cameracontrols(scene.children[1])
    eyeposition = Float32[2, 2, 1.3]
    lookat = Float32[0.82, 0.82, 0.1]
    # Title
    supertitle =
        layout[1, 2:4] = Label(scene, titlename, textsize = 50, color = :black)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end
    lowerclim_string = @lift(
        "lower clim quantile = " *
        @sprintf("%0.2f", $lowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[1])
    )
    upperclim_string = @lift(
        "upper clim quantile = " *
        @sprintf("%0.2f", $upperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[2])
    )
    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )
    layout[1, 1] = Label(scene, "Menu", width = width, textsize = 50)

    # Modify Axis
    axis = scene.children[1][OldAxis]
    # axis[:names][:axisnames] = ("↓ Zonal [m] ", "Meriodonal [m]↓ ", "Depth [m]↓ ")
    axis[:names][:axisnames] = ("↓", "↓ ", "↓ ")
    axis[:names][:align] =
        ((:left, :center), (:right, :center), (:right, :center))
    # need to adjust size of ticks first and then size of axis names
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)
    # axis[:ticks][:ranges_labels].val # current axis labels
    xticks = collect(range(-0, aspect[1], length = 2))
    yticks = collect(range(-0, aspect[2], length = 6))
    zticks = collect(range(-0, aspect[3], length = 2))
    ticks = (xticks, yticks, zticks)
    axis[:ticks][:ranges] = ticks
    xtickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    xtickslabels[end] = "1e6"
    ytickslabels = ["", "south", "", "", "north", ""]
    ztickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    labels = (xtickslabels, ytickslabels, ztickslabels)
    axis[:ticks][:labels] = labels

    display(scene)
    # Change the default camera position after the fact
    # note that these change dynamically as the plot is manipulated
    return scene
end

# 2D visualization
function visualize(
    states::Array{Array{S, 2}, 1};
    statenames = string.(1:length(states)),
    units = ["" for i in eachindex(states)],
    aspect = (1, 1, 1),
    resolution = (2412, 1158),
    title = "Zonal and Temporal Average of ",
    xlims = (0, 1),
    ylims = (0, 1),
    bins = 300,
) where {S}
    # Create scene
    scene, layout = layoutscene(resolution = resolution)
    lscene =
        layout[2:4, 2:4] = Axis(
            scene,
            xlabel = "South to North [m]",
            xlabelcolor = :black,
            ylabel = "Depth [m]",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = 25,
            yticklabelsize = 25,
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
            titlesize = 50,
        )
    width = round(Int, resolution[1] / 4) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    interpolationlabels = ["contour", "heatmap"]
    interpolationchoices = [true, false]
    interpolationnode = Node(interpolationchoices[1])

    # Statistics
    llscene =
        layout[4, 1] = Axis(
            scene,
            xlabel = @lift(statenames[$statenode] * " " * units[$statenode]),
            xlabelcolor = :black,
            ylabel = "pdf",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = 25,
            yticklabelsize = 25,
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
        )
    layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)

    # Clim sliders
    upperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    #ylims = @lift(range($lowerval, $upperval, length = $))
    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    unit = @lift(units[$statenode])
    oclims = @lift((
        quantile($state[:], $lowerclim_node),
        quantile($state[:], $upperclim_node),
    ))
    cmap_rgb = colornode
    clims = @lift(
        $oclims[1] != $oclims[2] ? (minimum($oclims), maximum($oclims)) :
        (minimum($oclims) - 1, maximum($oclims) + 1)
    )
    xlims = Array(range(xlims[1], xlims[2], length = 4)) #collect(range(xlims[1], xlims[2], length = size(states[1])[1]))
    ylims = Array(range(ylims[1], ylims[2], length = 4)) #@lift(collect(range($lowerval], $upperval, length = size($state)[2])))
    # newrange = @lift(range($lowerval, $upperval, length = 4))
    # lscene.yticks = @lift(Array($newrange))
    titlename = @lift(title * $statename) # use padding and appropriate centering
    layout[1, 2:4] = Label(scene, titlename, textsize = 50)
    # heatmap
    heatmap1 = heatmap!(
        lscene,
        xlims,
        ylims,
        state,
        interpolate = interpolationnode,
        colormap = cmap_rgb,
        colorrange = clims,
    )

    # statistics
    histogram_node = @lift(histogram($state, bins = bins))
    xs = @lift($histogram_node[1])
    ys = @lift($histogram_node[2])
    pdf = GLMakie.AbstractPlotting.barplot!(
        llscene,
        xs,
        ys,
        color = :red,
        strokecolor = :red,
        strokewidth = 1,
    )
    @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
    @lift(GLMakie.AbstractPlotting.ylims!(llscene, extrema($histogram_node[2])))
    vlines!(llscene, @lift($clims[1]), color = :black, linewidth = width / 100)
    vlines!(llscene, @lift($clims[2]), color = :black, linewidth = width / 100)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end

    interpolationmenu =
        Menu(scene, options = zip(interpolationlabels, interpolationchoices))
    on(interpolationmenu.selection) do s
        interpolationnode[] = s
        heatmap1 = heatmap!(
            lscene,
            xlims,
            ylims,
            state,
            interpolate = s,
            colormap = cmap_rgb,
            colorrange = clims,
        )
    end

    newlabel = @lift($statename * " " * $unit)
    cbar = Colorbar(scene, heatmap1, label = newlabel)
    cbar.width = Relative(1 / 3)
    cbar.height = Relative(5 / 6)
    cbar.halign = :center
    # cbar.flipaxisposition = true
    # cbar.labelpadding = -350
    cbar.labelsize = 50

    lowerclim_string = @lift(
        "clim quantile = " *
        @sprintf("%0.2f", $lowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[1])
    )
    upperclim_string = @lift(
        "clim quantile = " *
        @sprintf("%0.2f", $upperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[2])
    )

    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(
            scene,
            "plotting options",
            width = width,
            textsize = 30,
            padding = (0, 0, 10, 0),
        ),
        interpolationmenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )

    layout[2:4, 5] = vgrid!(
        Label(
            scene,
            "Color Bar",
            width = width / 2,
            textsize = 50,
            padding = (25, 0, 0, 00),
        ),
        cbar,
    )
    layout[1, 1] = Label(scene, "Menu", width = width, textsize = 50)
    display(scene)
    return scene
end


text:
# Example Configuration for nearly neutral simuation
# Command line to run:  julia run.jl examples/ia_example3
# Creates output file examples/nn_example3.csv

function dfe_mixd( x::Int64 )
  dfe_mixed( x::Int64; adv_probability=0.01,
    alpha_disadv=dfe_disadv_alpha, beta_disadv=dfe_disadv_beta,
    alpha_adv=dfe_adv_alpha, beta_adv=dfe_adv_beta )
end
const nn_simtype = 1    # nearly neutral infinite alleles
mu_list_flag =true
type_str="I1"
@everywhere const N_list     = [100,200]   # popsize N
const N_mu_list= [1.0,2.0]   # Mutation rate as a multiple of 1.0/N
const mu_list = [0.1, 0.02]
const L = 1           # number of loci # not used in infinite alleles
const ngens = 1000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of N
const use_poplist=false  # set to false for "production" runs with large ngens
dfe=dfe_mixd
dfe_adv_prob=0.01
dfe_adv_alpha=1.0
dfe_disadv_alpha=0.2
dfe_adv_beta=0.005
dfe_disadv_beta=0.5
dfe_str = "dfe_mixd advprob=0.01 alpha_adv=1.0 alpha_disadv=0.2 beta_adv=0.005 beta_disadv=0.5"

ext: jl

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
struct RaggedArray{TO, TV}
    offsets::TO
    values::TV
end
whichdevice(ra::RaggedArray) = whichdevice(ra.offsets)

function RaggedArray(vv::Vector{Vector{T}}) where T
    offsets = [1; 1 .+ accumulate(+, collect(length(v) for v in vv))]
    values = Vector{T}(undef, offsets[end]-1)
    for (i, v) in enumerate(vv)
        r = offsets[i]:offsets[i+1]-1
        values[r] .= v
    end
    RaggedArray(offsets, values)
end
@define_cu(RaggedArray, :offsets, :values)

function Base.getindex(ra::RaggedArray, i)
    @assert 1 <= i < length(ra.offsets)
    r = ra.offsets[i]:ra.offsets[i+1]-1
    @assert 1 <= r[1] && r[end] <= length(ra.values)
    return @view ra.values[r]
end
function Base.getindex(ra::RaggedArray, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    return ra.values[r[i]]
end
function Base.setindex!(ra::RaggedArray, v, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    ra.values[r[i]] = v
end

function find_varind(black, white, ::Type{TI}=Int) where TI
    nel = length(black)
    nel == length(white) || throw("Black and white vectors should be of the same length")
    varind = zeros(TI, nel)
    k = 1
    for i in 1:nel
        if !black[i] && !white[i]
            varind[i] = k
            k += 1
        end
    end
    return varind
end

function find_black_and_white(dh)
    black = falses(getncells(dh.grid))
    white = falses(getncells(dh.grid))
    if haskey(dh.grid.cellsets, "black")
        for c in grid.cellsets["black"]
            black[c] = true
        end
    end
    if haskey(dh.grid.cellsets, "white")
        for c in grid.cellsets["white"]
            white[c] = true
        end
    end

    return black, white
end

YoungsModulus(p) = getE(p)
PoissonRatio(p) = getν(p)

function compliance(Ke, u, dofs)
    comp = zero(eltype(u))
    for i in 1:length(dofs)
        for j in 1:length(dofs)
            comp += u[dofs[i]]*Ke[i,j]*u[dofs[j]]
        end
    end
    comp
end

function meandiag(K::AbstractMatrix)
    z = zero(eltype(K))
    for i in 1:size(K, 1)
        z += abs(K[i, i])
    end
    return z / size(K, 1)
end

density(var, xmin) = var*(1-xmin) + xmin

macro debug(expr)
    return quote
        if DEBUG[]
            $(esc(expr))
        end
    end
end


text:
# Replicate the Biogeme Swissmetro example
# https://biogeme.epfl.ch/examples/swissmetro/01logit.html

using Distributed
procs = addprocs(4)
Distributed.@everywhere begin
    using Pkg; Pkg.activate(joinpath(@__DIR__, "..", ".."))
end

@everywhere using Dagger, DiscreteChoiceModels

@testset "Biogeme swissmetro distributed" begin
    using Dagger
    using DiscreteChoiceModels
    using Test
    using StatsBase
    using CSV

    @test length(workers()) == 4

    data = DTable(
        CSV.File(joinpath(dirname(Base.source_path()), "../data/biogeme_swissmetro.dat"), delim='\t'),
        2683
    )

    @test length(data.chunks) == 4

    # filter to wanted cases
    data = filter(data) do row
        ((row.PURPOSE == 1) | (row.PURPOSE == 3)) && row.CHOICE != 0
    end

    @test length(data) == 6768

    data = map(data) do r
        merge((
            avtr=(r.TRAIN_AV == 1) && (r.SP != 0),
            avsm=r.SM_AV == 1,
            avcar=(r.CAR_AV == 1) && (r.SP != 0)
        ), pairs(r))
    end

    model = multinomial_logit(
            @utility(begin
                1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
                2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
                3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

                αswissmetro = 0, fixed  # fix swissmetro ASC to zero
            end),
            :CHOICE,
            data,
            availability=[
                1 => :avtr,
                2 => :avsm,
                3 => :avcar,
            ]
        )

    coefs = Dict(zip(coefnames(model), round.(coef(model), sigdigits=3)))
    ses = Dict(zip(coefnames(model), round.(stderror(model), sigdigits=3)))

    @test coefs[:αcar] ≈ -0.155 # biogeme says 0.155 which would round to 0.16, but it's actually 0.1546
    @test coefs[:αswissmetro] == 0  # fixed
    @test coefs[:αtrain] ≈ -0.701
    @test coefs[:βcost] ≈ -1.08
    @test coefs[:βtravel_time] ≈ -1.28

    @test ses[:αcar] ≈ 0.0432
    @test isnan(ses[:αswissmetro])  # fixed
    @test ses[:αtrain] ≈ 0.0549
    @test ses[:βcost] ≈ 0.0518
    @test ses[:βtravel_time] ≈ 0.0569

    idx = Dict(zip(coefnames(model), 1:length(coefs)))
    vcovmat = round.(vcov(model), sigdigits=3)
    @test vcovmat[idx[:αtrain], idx[:αcar]] ≈ 0.00138
    @test vcovmat[idx[:βcost], idx[:αcar]] ≈ 0.000485
    @test vcovmat[idx[:βcost], idx[:αtrain]] ≈ 8.22e-6
    @test vcovmat[idx[:βtravel_time], idx[:αcar]] ≈ -0.00144
    @test vcovmat[idx[:βtravel_time], idx[:αtrain]] ≈ -0.00225
    @test vcovmat[idx[:βtravel_time], idx[:βcost]] ≈ 0.000550

    @test all(isnan.(vcovmat[idx[:αswissmetro], :]))
    @test all(isnan.(vcovmat[:, idx[:αswissmetro]]))

    @test round(loglikelihood(model), digits=3) ≈ -5331.252
    @test round(model.init_ll, digits=3) ≈ -6964.663

    sleep(2)  # let dagger shut down
end

# clean up
rmprocs(procs)


text:
function ueberpruefeMittelpunkte(PX1,PX2,theta,ParaX1,ParaX2,a,b,X1,X2,N,JJ,VX1,VX2,HX1,HX2,Abstand,MitteSuSx,MitteSuSy,nx,ny,Pixel)
counterb = 2;
eps = 1e-7;
XWerte = zeros(size(Abstand,1),size(Abstand,2));
YWerte = zeros(size(Abstand,1),size(Abstand,2));
XWertevec = zeros(size(Abstand,1)*size(Abstand,2),1);
YWertevec = zeros(size(Abstand,1)*size(Abstand,2),1);

t1(Z,S,k) = (X[k]-PX1[Z,S])/cosd(theta[S]);
t2(Z,S,k) = (Y[k]-PX2[Z,S])/sind(theta[S]);

 SPX1,SPX2,Plot = berechneallgemeineSchnittpunkte(PX1,PX2,theta,ParaX1,ParaX2,a,b,X1,X2,N,JJ,VX1,VX2,HX1,HX2,counterb);

Z1 = 1;
S1 = 1;
Mittex = zeros(size(Abstand,1),size(Abstand,2));
Mittey = zeros(size(Abstand,1),size(Abstand,2));

if MitteSuSx(Z1,S1) == SPX1[Z1,S1]+1/2*(SPX1[Z1+1,S1]-SPX1[Z1,S1]) && MitteSuSy(Z1,S1) == SPX2[Z1,S1]+1/2*(SPX2[Z1+1,S1]-SPX2[Z1,S1])
    Mittex(Z,S) = SPX1[Z,S]+1/2*(SPX1[Z+1,S]-SPX1[Z,S]);
    Mittey(Z,S) = SPX2[Z,S]+1/2*(SPX2[Z+1,S]-SPX2[Z,S]);

veccounter = 1;
        for m = 1:JJ# Anzahl Winkel
            for Spalte = m*N-(N-1):m*N
                    counter = 1;
                    for Zeile = 1:size(Abstand,1)
                        if almostequal(Abstand[Zeile,Spalte], 0)
                        else
                            if (Mittex(Zeile,Spalte) >= a-eps) && (Mittex(Zeile,Spalte) <= b+eps) && (Mittey(Zeile,Spalte) >= a-eps) && (Mittey(Zeile,Spalte) <= b+eps)
                                XWerte[counter,Spalte] = Mittex(Zeile,Spalte);
                                YWerte[counter,Spalte] = Mittey(Zeile,Spalte);
                                XWertevec[veccounter] = Mittex(Zeile,Spalte);
                                YWertevec[veccounter] = Mittey(Zeile,Spalte);
                                veccounter = veccounter+1;
                                counter = counter+1;
                            end
                        end
                    end
            end
        end
        scatter!(XWertevec,YWertevec,markersize=3,markercolor=:cyan);
        scatter!(Pixel[:,2],Pixel[:,3],markersize=3,markercolor=:yellow);
else
    println("Überprüft erneut eure Formeln.")
    XWerte = [];
    YWerte = [];
end

print("MPX1: ");
show(stdout,"text/plain",XWerte);
print("\n\n MPX2: ");
show(stdout,"text/plain",YWerte);
display(Plot);

return XWerte, YWerte;

end


text:
struct Component{P,S}
    prior::P
    suffstats::S
end

Component(prior::Distribution) = Component(prior, empty_suffstats(prior))

Component(params::NTuple{N,<:Real}) where N = Component(NormalInverseChisq(params...))

Base.show(io::IO, c::Component{P,S}) where {P,S} = print(io, "$P$(params(c)) w/ n=$(nobs(c))")

posterior_predictive(gc::Component) =
    posterior_predictive(posterior_canon(gc.prior, gc.suffstats))

posterior_predictive(d::NormalInverseChisq) =
    LocationScale(d.μ, sqrt((1+d.κ)*d.σ2/d.κ), TDist(d.ν))

function posterior_predictive(d::NormalInverseWishart)
    df = d.nu - length(d.mu) + 1
    Λ = convert(AbstractMatrix, d.Lamchol) * (d.kappa+1)/(d.kappa*df)
    Λchol = cholesky(Λ)
    Distributions.GenericMvTDist(df, d.mu, PDMat(Λ, Λchol))
end

Distributions.logpdf(c::Component, x) = logpdf(posterior_predictive(c), x)
Distributions.pdf(c::Component, x) = logpdf(posterior_predictive(c), x)

Distributions.params(c::Component) = params(posterior_canon(c.prior, c.suffstats))

empty_suffstats(d::NormalInverseChisq) = NormalStats(0, 0, 0, 0)
function empty_suffstats(d::NormalInverseWishart)
    d = d.dim
    MvNormalStats(zeros(d), zeros(d), zeros(d,d), 0.)
end

# struct NormalStats <: SufficientStats
#     s::Float64    # (weighted) sum of x
#     m::Float64    # (weighted) mean of x
#     s2::Float64   # (weighted) sum of (x - μ)^2
#     tw::Float64    # total sample weight
# end

# based on teh wikipedia algorithm for online variance/mean
function add(ss::NormalStats, x)
    tw = ss.tw+1
    Δ = x - ss.m
    m = ss.m + Δ/tw
    s2 = tw>1 ? ss.s2 + Δ*(x-m) : 0
    return NormalStats(ss.s+x, m, s2, tw)
end

function sub(ss::NormalStats, x)
    tw = ss.tw-1
    if tw ≤ 0
        return NormalStats(0,0,0,0)
    else
        Δ = x - ss.m
        m = ss.m - Δ/tw
        s2 = ss.s2 - Δ*(x-m)
        return NormalStats(ss.s-x, m, s2, tw)
    end
end

function add(ss::MvNormalStats, x::AbstractVector)
    tw = ss.tw + 1
    Δ = x - ss.m
    m = ss.m + Δ/tw
    s2 = tw>1 ? ss.s2 + Δ*(x.-m)' : zeros(size(ss.s2))
    return MvNormalStats(ss.s+x, m, Hermitian(s2), tw)
end

function sub(ss::MvNormalStats, x::AbstractVector)
    tw = ss.tw-1
    if tw ≤ 0
        n = length(ss.m)
        return MvNormalStats(zeros(n), zeros(n), zeros(n,n), 0)
    else
        Δ = x - ss.m
        m = ss.m - Δ./tw
        s2 = ss.s2 - Δ*(x-m)'
        return MvNormalStats(ss.s-x, m, s2, tw)
    end
end

add(c::Component, x) = Component(c.prior, add(c.suffstats, x))
sub(c::Component, x) = Component(c.prior, sub(c.suffstats, x))

nobs(c::Component) = Int(c.suffstats.tw)
Base.isempty(c::Component) = nobs(c) == 0

"""
    marginal_lhood(c::Component)

The marginal likelihood of data fit so far by `o`, which is the integral of the
likelihood given the mean and variance under the prior on those parameters.
"""
marginal_lhood(c::Component) =
    marginal_lhood(c.prior, c.suffstats)

marginal_lhood(prior, suffstats) = exp(marginal_log_lhood(prior, suffstats))

function marginal_lhood(prior::NormalInverseChisq{Float64}, suffstats::NormalStats)
    μ0, σ20, κ0, ν0 = params(prior)
    μn, σ2n, κn, νn = params(posterior_canon(prior, suffstats))
    n = νn - ν0
    gamma(νn*0.5)/gamma(ν0*0.5) *
        sqrt(κ0/κn) *
        (ν0*σ20)^(ν0*0.5) * (νn*σ2n)^(-νn*0.5) /
        π^(n*0.5)
end

"""
    marginal_log_lhood(c::Component)

The marginal likelihood of data fit so far by `o`, which is the integral of the
likelihood given the mean and variance under the prior on those parameters.
"""
marginal_log_lhood(c::Component) =
    marginal_log_lhood(c.prior, c.suffstats)

function marginal_log_lhood(prior::NormalInverseChisq{Float64}, suffstats::NormalStats)
    μ0, σ20, κ0, ν0 = params(prior)
    μn, σ2n, κn, νn = params(posterior_canon(prior, suffstats))
    n = νn - ν0
    lgamma(νn*0.5) - lgamma(ν0*0.5) +
        0.5*(log(κ0)-log(κn)) +
        (0.5*ν0)*(log(ν0)+log(σ20)) - (0.5*νn)*(log(νn)+log(σ2n)) -
        (0.5*n)*log(π)
end

function marginal_log_lhood(prior::NormalInverseWishart, suffstats::MvNormalStats)
    μ0, Λchol0, κ0, ν0 = params(prior)
    μn, Λcholn, κn, νn = params(posterior_canon(prior, suffstats))
    n = νn - ν0
    d = length(μ0)

    (-0.5*n*d) * logπ +
        logmvgamma(d, νn*0.5) - logmvgamma(d, ν0*0.5) +
        logdet(Λchol0)*ν0*0.5 - logdet(Λcholn)*νn*0.5 +
        d*0.5*(log(κ0) - log(κn))
end


text:
strip_nothing_(::Type{Union{Nothing, T}}) where T = T
strip_nothing_(::Type{T}) where T = T

progress_noop_(n) = nothing

mutable struct SplitCapOutput{S}
    work_buffer::Vector{UInt8}
    complete_buffers::Channel{Vector{UInt8}}
    stream::S
end

# Since Julia (as of 1.5) doesn't support task migration,
# continue with new task after each buffer write, to rebalance across threads
function write_one_and_continue_(output::SplitCapOutput, free_buffers::Channel, pending::Threads.Atomic{Int})
    i = iterate(output.complete_buffers)
    if i === nothing
        Threads.atomic_sub!(pending, 1)
        return nothing
    end
    b, _ = i
    write(output.stream, b)
    empty!(b)
    put!(free_buffers, b)
    Threads.@spawn write_one_and_continue_($output, $free_buffers, $pending)
end

function splitcap(
    ::Type{KeyType},
    ::Type{StreamType},
    reader::PcapReader,
    record2key,
    key2stream,
    progress_callback = progress_noop_
) where {KeyType, StreamType}
    buffer_size = 1024 * 1024 * 2
    max_pending_buffers = 4
    outputs = Dict{KeyType, SplitCapOutput{StreamType}}()
    free_buffers = Channel{Vector{UInt8}}(Inf)
    n = 0
    pending = Threads.Atomic{Int}(0)
    try
        while !eof(reader)
            record = read(reader)
            dst = record2key(record)
            if dst isa KeyType
                output = get!(outputs, dst) do
                    stream = key2stream(dst)
                    buffer = sizehint!(UInt8[], buffer_size + 1500)
                    append!(buffer, reader.raw_header)
                    output = SplitCapOutput{StreamType}(
                        buffer,
                        Channel{Vector{UInt8}}(max_pending_buffers),
                        stream)
                    Threads.atomic_add!(pending, 1)
                    Threads.@spawn write_one_and_continue_($output, $free_buffers, $pending)
                    output
                end
                append!(output.work_buffer, record.raw)
                if length(output.work_buffer) >= buffer_size
                    put!(output.complete_buffers, output.work_buffer)
                    if isready(free_buffers)
                        output.work_buffer = take!(free_buffers)
                    else
                        output.work_buffer = sizehint!(UInt8[], buffer_size + 1500)
                    end
                end
            end
            n += 1
            progress_callback(n)
            GC.safepoint()
        end
        for output in values(outputs)
            if !isempty(output.work_buffer)
                put!(output.complete_buffers, output.work_buffer)
            end
            close(output.complete_buffers)
        end
        while pending[] != 0
            sleep(0.1)
        end
    finally
        for output in values(outputs)
            close(output.stream)
        end
    end
    nothing
end

function splitcap(
    reader::PcapReader,
    record2key,
    key2stream,
    progress_callback = progress_noop_
)
    KeyType = strip_nothing_(Core.Compiler.return_type(record2key, (PcapRecord,)))
    StreamType = Core.Compiler.return_type(key2stream, (KeyType,))
    splitcap(KeyType, StreamType, reader, record2key, key2stream, progress_callback)
end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################
"""
    constraint_mc_residual

Equality constraint that describes the residual definition, which depends on the
criterion assigned to each individual measurement in data["meas"]["m"]["crit"].
"""
function constraint_mc_residual(pm::_PMD.AbstractUnbalancedPowerModel, i::Int; nw::Int=_IM.nw_id_default)

    cmp_id = get_cmp_id(pm, nw, i)
    res = _PMD.var(pm, nw, :res, i)
    var = _PMD.var(pm, nw, _PMD.ref(pm, nw, :meas, i, "var"), cmp_id)
    dst = _PMD.ref(pm, nw, :meas, i, "dst")
    rsc = _PMD.ref(pm, nw, :se_settings)["rescaler"]
    crit = _PMD.ref(pm, nw, :meas, i, "crit")
    conns = get_active_connections(pm, nw, _PMD.ref(pm, nw, :meas, i, "cmp"), cmp_id)

    for (idx, c) in enumerate(conns)
        if (occursin("ls", crit) || occursin("lav", crit)) && isa(dst[idx], _DST.Normal)
            μ, σ = occursin("w", crit) ? (_DST.mean(dst[idx]), _DST.std(dst[idx])) : (_DST.mean(dst[idx]), 1.0)
        end
        if isa(dst[idx], Float64)
            JuMP.@constraint(pm.model, var[c] == dst[idx])
            JuMP.@constraint(pm.model, res[idx] == 0.0)
        elseif crit ∈ ["wls", "ls"] && isa(dst[idx], _DST.Normal)
            JuMP.@constraint(pm.model,
                res[idx] * rsc^2 * σ^2 == (var[c] - μ)^2
            )
        elseif crit == "rwls" && isa(dst[idx], _DST.Normal)
            JuMP.@constraint(pm.model,
                res[idx] * rsc^2 * σ^2 >= (var[c] - μ)^2
            )
        elseif crit ∈ ["wlav", "lav"] && isa(dst[idx], _DST.Normal)
            JuMP.@NLconstraint(pm.model,
                res[idx] * rsc * σ == abs(var[c] - μ)
            )
        elseif crit == "rwlav" && isa(dst[idx], _DST.Normal)
            JuMP.@constraint(pm.model,
                res[idx] * rsc * σ >= (var[c] - μ)
            )
            JuMP.@constraint(pm.model,
                res[idx] * rsc * σ >= - (var[c] - μ)
            )
        elseif crit == "mle"
            #TODO: enforce min and max in the meas dictionary and just with haskey make it optional for extendedbeta
            pkg_id = any([ dst[idx] isa d for d in [ExtendedBeta{Float64}, _Poly.Polynomial]]) ? _PMDSE : _DST
            lb = ( !isa(dst[idx], _DST.MixtureModel) && !isinf(pkg_id.minimum(dst[idx])) ) ? pkg_id.minimum(dst[idx]) : -10
            ub = ( !isa(dst[idx], _DST.MixtureModel) && !isinf(pkg_id.maximum(dst[idx])) ) ? pkg_id.maximum(dst[idx]) : 10
            if any([ dst[idx] isa d for d in [_DST.MixtureModel, _Poly.Polynomial]]) lb = _PMD.ref(pm, nw, :meas, i, "min") end
            if any([ dst[idx] isa d for d in [_DST.MixtureModel, _Poly.Polynomial]]) ub = _PMD.ref(pm, nw, :meas, i, "max") end

            shf = abs(Optim.optimize(x -> -pkg_id.logpdf(dst[idx],x),lb,ub).minimum)
            f = Symbol("df_",i,"_",c)

            fun(x) = rsc * ( - shf + pkg_id.logpdf(dst[idx],x) )
            grd(x) = pkg_id.gradlogpdf(dst[idx],x)
            hes(x) = heslogpdf(dst[idx],x)
            JuMP.register(pm.model, f, 1, fun, grd, hes)
            JuMP.add_NL_constraint(pm.model, :($(res[idx]) == - $(f)($(var[c]))))
        else
            error("SE criterion of measurement $(i) not recognized")
        end
    end
end


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking
#gr(size=(600,600));

# ### snippet 4.1

# No attempt has been made to condense this too fewer lines of code

noofsteps = 20;
noofwalks = 15;
pos = Array{Float64, 2}(rand(Uniform(-1, 1), noofsteps, noofwalks));
pos[1, :] = zeros(noofwalks);
csum = cumsum(pos, dims=1);

#-

f = Plots.font("DejaVu Sans", 6)
mx = minimum(csum) * 0.9

# Plot and annotate the random walks

p1 = plot(csum, leg=false, title="Random walks ($(noofwalks))")
plot!(p1, csum[:, Int(floor(noofwalks/2))], leg=false, title="Random walks ($(noofwalks))", color=:black)
plot!(p1, [5], seriestype="vline")
annotate!(5, mx, text("step 4", f, :left))
plot!(p1, [9], seriestype="vline")
annotate!(9, mx, text("step 8", f, :left))
plot!(p1, [17], seriestype="vline")
annotate!(17, mx, text("step 16", f, :left))

# Generate 3 plots of densities at 3 different step numbers (4, 8 and 16)

p2 = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3);
plt = 1
for step in [4, 8, 16]
  indx = step + 1 # We aadded the first line of zeros
  global plt
  fit = fit_mle(Normal, csum[indx, :])
  x = (fit.μ-4*fit.σ):0.01:(fit.μ+4*fit.σ)
  p2[plt] = density(csum[indx, :], legend=false, title="$(step) steps")
  plot!( p2[plt], x, pdf.(Normal( fit.μ , fit.σ ) , x ), fill=(0, .5,:orange))
  plt += 1
end
p3 = plot(p2..., layout=(1, 3))
plot(p1, p3, layout=(2,1))
#src savefig("fig4_2.pdf")

# ### snippet 4.2

prod(1 .+ rand(Uniform(0, 0.1), 10))

# ### snippet 4.3

growth = [prod(1 .+ rand(Uniform(0, 0.1), 10)) for i in 1:10000];
fit = fit_mle(Normal, growth)
plot(Normal(fit.μ , fit.σ ), fill=(0, .5,:orange), lab="Normal distribution")
density!(growth, lab="'sample' distribution")
savefig("s4_3.pdf") #src

# ### snippet 4.4

big = [prod(1 .+ rand(Uniform(0, 0.5), 12)) for i in 1:10000];
small = [prod(1 .+ rand(Uniform(0, 0.01), 12)) for i in 1:10000];
fitb = fit_mle(Normal, big)
fits = fit_mle(Normal, small)
p1 = plot(Normal(fitb.μ , fitb.σ ), lab="Big normal distribution", fill=(0, .5,:orange))
p2 = plot(Normal(fits.μ , fits.σ ), lab="Small normal distribution", fill=(0, .5,:orange))
density!(p1, big, lab="'big' distribution")
density!(p2, small, lab="'small' distribution")
plot(p1, p2, layout=(1, 2))
#src savefig("s4_4.pdf")

# ### snippet 4.5

log_big = [log(prod(1 .+ rand(Uniform(0, 0.5), 12))) for i in 1:10000];
fit = fit_mle(Normal, log_big)
plot(Normal(fit.μ , fit.σ ), fill=(0, .5,:orange), lab="Normal distribution")
density!(log_big, lab="'sample' distribution")
#src savefig("s4_5.pdf")

# ### snippet 4.6

# Grid of 1001 steps

p_grid = range(0, step=0.001, stop=1);

# all priors = 1.0

prior = ones(length(p_grid));

# Binomial pdf

likelihood = [pdf(Binomial(9, p), 6) for p in p_grid];

# As Uniform priar has been used, unstandardized posterior is equal to likelihood

posterior = likelihood .* prior;

# Scale posterior such that they become probabilities

posterior = posterior / sum(posterior);

# Sample using the computed posterior values as weights
# In this example we keep the number of samples equal to the length of p_grid,
# but that is not required.

samples = sample(p_grid, Weights(posterior), length(p_grid));
p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:length(p_grid), samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")

# analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")

# quadratic approximation

plot!( p[2], x, pdf.(Normal( 0.67 , 0.16 ) , x ), lab="Normal approximation", fill=(0, .5,:orange))
plot(p..., layout=(1, 2))

# End of `clip-01-06.jl`


text:
using Mimi

include("utils/load_parameters.jl")
include("components/RCPSSPScenario.jl")
include("components/CO2emissions.jl")
include("components/CO2cycle.jl")
include("components/CO2forcing.jl")
include("components/CH4emissions.jl")
include("components/CH4cycle.jl")
include("components/CH4forcing.jl")
include("components/N2Oemissions.jl")
include("components/N2Ocycle.jl")
include("components/N2Oforcing.jl")
include("components/LGemissions.jl")
include("components/LGcycle.jl")
include("components/LGforcing.jl")
include("components/SulphateForcing.jl")
include("components/TotalForcing.jl")
include("components/ClimateTemperature.jl")
include("components/PermafrostSiBCASA.jl")
include("components/PermafrostJULES.jl")
include("components/PermafrostTotal.jl")

include("models/climate_model_def.jl")


text:
# An SVG backend for compose.

export SVG

# Format a floating point number into a decimal string of reasonable precision.
function svg_fmt_float(x::Float64)
    # All svg (in our use) coordinates are in millimeters. This number gives the
    # largest deviation from the true position allowed in millimeters.
    const eps = 0.01
    a = @sprintf("%0.8f", round(x / eps) * eps)
    n = length(a)

    while a[n] == '0'
        n -= 1
    end

    if a[n] == '.'
        n -= 1
    end

    a[1:n]
end

# Format a color for SVG.
svg_fmt_color(c::ColorValue) = @sprintf("#%s", hex(c))
svg_fmt_color(c::Nothing) = "none"

# The SVG image we generate are going to use millimeters everywhere.
type SVGMeasure <: NativeMeasure
    value::Float64
end

*(u::Float64, v::SVGMeasure) = SVGMeasure(u * v.value)
+(u::SVGMeasure, v::SVGMeasure) = SVGMeasure(u.value + v.value)
-(u::SVGMeasure, v::SVGMeasure) = SVGMeasure(u.value - v.value)
convert(::Type{Float64}, u::SVGMeasure) = u.value
convert(::Type{SVGMeasure}, u::Float64) = SVGMeasure(u)
function convert(::Type{SimpleMeasure{MillimeterUnit}}, u::SVGMeasure)
    SimpleMeasure{MillimeterUnit}(u.value)
end

type SVG <: Backend
    # Image size in millimeters.
    width::SVGMeasure
    height::SVGMeasure

    # Output stream.
    f::IO

    # Should f be closed when the backend is finished?
    close_stream::Bool

    # Current level of indentation.
    indentation::Int

    # Javascript fragments are placed in a function with a unique name. This is
    # a map of unique function names to javascript code.
    scripts::Dict{String, String}

    # Clip-paths that need to be defined at the end of the document.
    clippaths::Vector{Vector{Point}}

    # Miscelaneous embedded objects included immediately before the </svg> tag,
    # such as extra javascript or css.
    embobj::Set{String}

    # Keep track of which properties that are push are empty to we can avoid
    # printiing them.
    empty_properties::Vector{Bool}

    # Keep track of which properties have links, so we know when to insert
    # closing </a> tags.
    linked_properties::Vector{Bool}

    # Keep track of mask objects.
    mask_properties::Vector{Bool}

    # Emit the graphic on finish when writing to a buffer.
    emit_on_finish::Bool

    function SVG(f::IO,
                 width::MeasureOrNumber,
                 height::MeasureOrNumber,
                 emit_on_finish::Bool=true)
        img = new()
        img.width  = native_measure(width,  img)
        img.height = native_measure(height, img)
        img.f = f
        img.close_stream = false
        img.indentation = 0
        img.clippaths = Array(Vector{Point}, 0)
        img.scripts = Dict{String, String}()
        img.embobj = Set{String}()
        img.empty_properties = Array(Bool, 0)
        img.linked_properties = Array(Bool, 0)
        img.mask_properties = Array(Bool, 0)
        img.emit_on_finish = emit_on_finish

        write(img.f, @sprintf("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" version=\"1.1\" width=\"%smm\" height=\"%smm\" viewBox=\"0 0 %s %s\" style=\"stroke:black;fill:black\" stroke-width=\"0.5\">\n",
              svg_fmt_float(img.width.value), svg_fmt_float(img.height.value),
              svg_fmt_float(img.width.value), svg_fmt_float(img.height.value)))
        img
    end

    # Write to a file.
    function SVG(filename::String,
                 width::MeasureOrNumber,
                 height::MeasureOrNumber)
        f = open(filename, "w")
        img = SVG(f, width, height)
        img.close_stream = true
        img
    end

    # Write to buffer.
    function SVG(width::MeasureOrNumber,
                 height::MeasureOrNumber;
                 emit_on_finish::Bool=true)
        img = SVG(IOString(), width, height, emit_on_finish=emit_on_finish)
        img.close_stream = false
        img
    end
end

function finish(img::SVG)
    for obj in img.embobj
        write(img.f, obj)
        write(img.f, "\n")
    end

    if length(img.scripts) > 0
        write(img.f, "<script type=\"application/ecmascript\"><![CDATA[\n")
        for (fn_name, js) in img.scripts
            @printf(img.f, "function %s(evt) {\n%s\n}\n\n", fn_name, js)
        end
        write(img.f, "]]></script>\n")
    end

    if length(img.clippaths) > 0
        write(img.f, "<defs>\n")
        for (i, clippath) in enumerate(img.clippaths)
            write(img.f, "<clipPath id=\"clippath$(i)\">\n  <path d=\"")
            print_svg_path(img.f, clippath)
            write(img.f, "\" />\n</clipPath\n>")
        end
        write(img.f, "</defs>\n")
    end

    write(img.f, "</svg>\n")
    if img.close_stream
        close(img.f)
    end

    # If we are writing to a buffer. Collect the string and emit it.
    if img.emit_on_finish && typeof(img.f) == IOString
        seek(img.f, 0)
        emit(Emitable("image/svg+xml", readall(img.f)))
        close(img.f)
    end
end

function root_box(img::SVG)
    NativeBoundingBox(
        SVGMeasure(0.0),
        SVGMeasure(0.0),
        img.width,
        img.height)
end

native_zero(::SVG) = SVGMeasure(0.0)

function indent(img::SVG)
    for i in 1:img.indentation
        write(img.f, "  ")
    end
end

function native_measure(u::SimpleMeasure{PixelUnit},
                        img::SVG)
    native_measure(convert(SimpleMeasure{MillimeterUnit}, u), img)
end

function native_measure(u::SimpleMeasure{MillimeterUnit},
                        img::SVG)
    SVGMeasure(u.value)
end

# Draw

function print_svg_path(out, points::Vector{Point})
    @printf(out, "M%s,%s L",
            svg_fmt_float(points[1].x.value),
            svg_fmt_float(points[1].y.value))
    for point in points[2:]
        @printf(out, " %s %s",
                svg_fmt_float(point.x.value),
                svg_fmt_float(point.y.value))
    end
end

function draw(img::SVG, form::Lines)
    n = length(form.points)
    if n <= 1; return; end
    indent(img)
    paths = make_paths(form.points)
    for path in paths
        write(img.f, "<path d=\"")
        print_svg_path(img.f, path)
        write(img.f, "\" />\n")
    end
end

function draw(img::SVG, form::Curve)
    indent(img)
    @printf(img.f, "<path d=\"M%s,%s C%s,%s %s,%s %s,%s\" />\n",
        svg_fmt_float(form.anchor0.x.value),
        svg_fmt_float(form.anchor0.y.value),
        svg_fmt_float(form.ctrl0.x.value),
        svg_fmt_float(form.ctrl0.y.value),
        svg_fmt_float(form.ctrl1.x.value),
        svg_fmt_float(form.ctrl1.y.value),
        svg_fmt_float(form.anchor1.x.value),
        svg_fmt_float(form.anchor1.y.value))
end

function draw(img::SVG, form::Polygon)
    n = length(form.points)
    if n <= 1; return; end

    indent(img)
    write(img.f, "<path d=\"")
    print_svg_path(img.f, form.points)
    write(img.f, " z\" />\n")
end

minmax(a, b) = a < b ? (a,b) : (b,a)

function draw(img::SVG, form::Ellipse)
    cx = form.center.x.value
    cy = form.center.y.value
    rx = sqrt((form.x_point.x.value - cx)^2 +
              (form.x_point.y.value - cy)^2)
    ry = sqrt((form.y_point.x.value - cx)^2 +
              (form.y_point.y.value - cy)^2)
    theta = radians2degrees(atan2(form.x_point.y.value - cy,
                                  form.x_point.x.value - cx))

    indent(img)
    eps = 1e-6

    if abs(rx - ry) < eps
        @printf(img.f, "<circle cx=\"%s\" cy=\"%s\" r=\"%s\" />\n",
                svg_fmt_float(cx), svg_fmt_float(cy), svg_fmt_float(rx))
    else
        @printf(img.f, "<ellipse cx=\"%s\" cy=\"%s\" rx=\"%s\" ry=\"%s\"",
                svg_fmt_float(cx), svg_fmt_float(cy),
                svg_fmt_float(rx), svg_fmt_float(ry))

        if abs(theta) >= eps
            @printf(img.f, " transform=\"rotate(%s %s %s)\"",
                    svg_fmt_float(theta),
                    svg_fmt_float(cx),
                    svg_fmt_float(cy))
        end

        write(img.f, " />\n")
    end
end

function draw(img::SVG, form::Text)
    indent(img)
    @printf(img.f, "<text x=\"%s\" y=\"%s\"",
            svg_fmt_float(form.pos.x.value),
            svg_fmt_float(form.pos.y.value))

    if is(form.halign, hcenter)
        print(img.f, " text-anchor=\"middle\"")
    elseif is(form.halign, hright)
        print(img.f, " text-anchor=\"end\"")
    end

    if is(form.valign, vcenter)
        print(img.f, " style=\"dominant-baseline:central\"")
    elseif is(form.valign, vtop)
        print(img.f, " style=\"dominant-baseline:text-before-edge\"")
    end

    if !isidentity(form.t)
        @printf(img.f, " transform=\"rotate(%s, %s, %s)\"",
                svg_fmt_float(radians2degrees(atan2(form.t.M[2,1], form.t.M[1,1]))),
                svg_fmt_float(form.pos.x.value),
                svg_fmt_float(form.pos.y.value))
    end

    # TODO: escape special characters
    @printf(img.f, ">%s</text>\n",
            pango_to_svg(form.value))
end

# Applying properties

function push_property(img::SVG, property::Property)
    if property === empty_property
        push!(img.empty_properties, true)
        push!(img.linked_properties, false)
        push!(img.mask_properties, false)
    else
        push!(img.empty_properties, false)
        indent(img)

        # There can only be one property of each type. E.g, defining 'fill'
        # twice in the same element is not valid svg.
        properties = Dict{Type, PropertyPrimitive}()
        rawproperties = {}
        p = property
        while !is(p, empty_property)
            if typeof(p.primitive) == SVGAttribute
                push!(rawproperties, p.primitive)
            else
                properties[typeof(p.primitive)] = p.primitive
            end
            p = p.next
        end

        # Special-case for the mask property.
        if haskey(properties, SVGDefineMask)
            write(img.f, @sprintf("<mask id=\"%s\">", properties[SVGDefineMask].id))
            push!(img.mask_properties, true)
        else
            push!(img.mask_properties, false)
        end

        # Special-case for links.
        if haskey(properties, SVGLink)
            write(img.f, @sprintf("<a xlink:href=\"%s\">",
                  properties[SVGLink].target === nothing ? '#' : properties[SVGLink].target))
            push!(img.linked_properties, true)
        else
            push!(img.linked_properties, false)
        end

        write(img.f, "<g")
        for p in chain(values(properties), rawproperties)
            apply_property(img, p)
        end
        write(img.f, ">\n")

        img.indentation += 1
    end
end

function pop_property(img::SVG)
    is_linked = pop!(img.linked_properties)
    is_mask = pop!(img.mask_properties)
    if !pop!(img.empty_properties)
        img.indentation -= 1
        indent(img)
        write(img.f, "</g>")
        if is_linked
            write(img.f, "</a>")
        end
        if is_mask
            write(img.f, "</mask>")
        end
        write(img.f, "\n")
    end
end

# Nop catchall
function apply_property(img::SVG, p::PropertyPrimitive)
end

function apply_property(img::SVG, p::Stroke)
    @printf(img.f, " stroke=\"%s\"", svg_fmt_color(p.value))
end

function apply_property(img::SVG, p::Fill)
    @printf(img.f, " fill=\"%s\"", svg_fmt_color(p.value))
end

function apply_property(img::SVG, p::LineWidth)
    @printf(img.f, " stroke-width=\"%s\"", svg_fmt_float(p.value.value))
end

function apply_property(img::SVG, p::Visible)
    @printf(img.f, " visibility=\"%s\"",
            p.value ? "visible" : "hidden")
end

function apply_property(img::SVG, p::Opacity)
    @printf(img.f, " opacity=\"%s\"", fmt_float(p.value))
end

function apply_property(img::SVG, p::SVGID)
    @printf(img.f, " id=\"%s\"", escape_string(p.value))
end

function apply_property(img::SVG, p::SVGClass)
    @printf(img.f, " class=\"%s\"", escape_string(p.value))
end

function apply_property(img::SVG, p::Font)
    @printf(img.f, " font-family=\"%s\"", escape_string(p.family))
end

function apply_property(img::SVG, p::FontSize)
    @printf(img.f, " font-size=\"%s\"", svg_fmt_float(p.value.value))
end

function apply_property(img::SVG, p::Clip)
    push!(img.clippaths, p.points)
    i = length(img.clippaths)
    write(img.f, " clip-path=\"url(#clippath$(i))\"")
end

function apply_property(img::SVG, p::SVGEmbed)
    add!(img.embobj, p.markup)
end

function apply_property(img::SVG, p::SVGMask)
    @printf(img.f, " mask=\"url(#%s)\"", p.id)
end

function apply_property(img::SVG, p::SVGAttribute)
    @printf(img.f, " %s=\"%s\"", p.attribute, p.value)
end

for event in events
    prop_name = lowercase(string(event))
    @eval begin
        function apply_property(img::SVG, p::($event))
            fn_name = add_script(img, p.value, object_id(p.value))
            @printf(img.f, " %s=\"%s(evt)\"", ($prop_name), fn_name)
        end
    end
end

# Add some javascript. Return the unique name generated for the wrapper function
# containing the given code.
function add_script(img::SVG, js::String, obj_id::Integer)
    fn_name = @sprintf("js_chunk_%x", obj_id)
    if !haskey(img.scripts, fn_name)
        img.scripts[fn_name] = replace(js, r"^"m, "  ")
    end
    fn_name
end


text:
@testset "Invalid Scenario Tree Checks" begin
    st = PH.ScenarioTree()
    PH.add_node(st, root(st))
    PH.add_node(st, root(st))

    @test_throws(ErrorException,
                 PH.solve_extensive(st, two_stage_model, ()->Ipopt.Optimizer()),
                 )
    @test_throws(ErrorException,
                 PH.solve(st, two_stage_model, PH.ScalarPenaltyParameter(2.0))
                 )
end

@testset "Probability Checks" begin
    st = PH.ScenarioTree()
    PH.add_leaf(st, root(st), 0.25)
    PH.add_leaf(st, root(st), 0.5)

    @test_throws(ErrorException,
                 PH.solve_extensive(st, two_stage_model, ()->Ipopt.Optimizer())
                 )
    @test_throws(ErrorException,
                 PH.solve(st, two_stage_model, PH.ScalarPenaltyParameter(2.0))
                 )

    st = PH.ScenarioTree()
    PH.add_leaf(st, root(st), 0.25)
    PH.add_leaf(st, root(st), 0.5)
    PH.add_leaf(st, root(st), 0.5)

    @test_throws(ErrorException,
                 PH.solve_extensive(st, two_stage_model, ()->Ipopt.Optimizer())
                 )
    @test_throws(ErrorException,
                 PH.solve(st, two_stage_model, PH.ScalarPenaltyParameter(2.0))
                 )
end


text:
using Test

@testset "AutoHOOT.jl" begin
    @testset "$filename" for filename in ["autodiff_test.jl", "graphops_test.jl"]
        println("Running $filename")
        include(filename)
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""HTTPGetAction describes an action based on HTTP Get requests.

    IoK8sApiCoreV1HTTPGetAction(;
        host=nothing,
        httpHeaders=nothing,
        path=nothing,
        port=nothing,
        scheme=nothing,
    )

    - host::String : Host name to connect to, defaults to the pod IP. You probably want to set \"Host\" in httpHeaders instead.
    - httpHeaders::Vector{IoK8sApiCoreV1HTTPHeader} : Custom headers to set in the request. HTTP allows repeated headers.
    - path::String : Path to access on the HTTP server.
    - port::IoK8sApimachineryPkgUtilIntstrIntOrString : Name or number of the port to access on the container. Number must be in the range 1 to 65535. Name must be an IANA_SVC_NAME.
    - scheme::String : Scheme to use for connecting to the host. Defaults to HTTP.
"""
mutable struct IoK8sApiCoreV1HTTPGetAction <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    httpHeaders::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1HTTPHeader} } # spec name: httpHeaders
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    port::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: port
    scheme::Any # spec type: Union{ Nothing, String } # spec name: scheme

    function IoK8sApiCoreV1HTTPGetAction(;host=nothing, httpHeaders=nothing, path=nothing, port=nothing, scheme=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("httpHeaders"), httpHeaders)
        setfield!(o, Symbol("httpHeaders"), httpHeaders)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("scheme"), scheme)
        setfield!(o, Symbol("scheme"), scheme)
        o
    end
end # type IoK8sApiCoreV1HTTPGetAction

const _property_map_IoK8sApiCoreV1HTTPGetAction = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("httpHeaders")=>Symbol("httpHeaders"), Symbol("path")=>Symbol("path"), Symbol("port")=>Symbol("port"), Symbol("scheme")=>Symbol("scheme"))
const _property_types_IoK8sApiCoreV1HTTPGetAction = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("httpHeaders")=>"Vector{IoK8sApiCoreV1HTTPHeader}", Symbol("path")=>"String", Symbol("port")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("scheme")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1HTTPGetAction }) = collect(keys(_property_map_IoK8sApiCoreV1HTTPGetAction))
Swagger.property_type(::Type{ IoK8sApiCoreV1HTTPGetAction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1HTTPGetAction[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1HTTPGetAction }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1HTTPGetAction[property_name]

function check_required(o::IoK8sApiCoreV1HTTPGetAction)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1HTTPGetAction }, name::Symbol, val)
end


text:
"""
	grdblend(cmd0::String="", arg1=nothing, arg2=nothing, kwargs...)

Reads a listing of grid files and blend parameters, or up to 2 GTMgrid types, and creates
a grid by blending the other grids using cosine-taper weights.

Full option list at [`grdblend`]($(GMTdoc)grdblend.html)

Parameters
----------

- $(GMT.opt_I)
    ($(GMTdoc)grdblend.html#i)
- $(GMT.opt_R)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdblend(....) form.
    ($(GMTdoc)grdblend.html#g)

- **C** | **clobber** :: [Type => Str | []]      ``Arg = f|l|o|u[±]``

    Clobber mode: Instead of blending, simply pick the value of one of the grids that covers a node.
    ($(GMTdoc)grdproject.html#c)
- **N** | **nodata** :: [Type => Str | Number]

    No data. Set nodes with no input grid to this value [Default is NaN].
    ($(GMTdoc)grdblend.html#n)
- **Q** | **headless** :: [Type => Bool]

    Create plain header-less grid file (for use with external tools). Requires that the output
    grid file is a native format (i.e., not netCDF). DO NOT USE WITH **G**.
    ($(GMTdoc)grdblend.html#q)
- **W** | **no_blend** :: [Type => Str | []]

    Do not blend, just output the weights used for each node [Default makes the blend].
    Append ``z`` to write the weight*z sum instead.
    ($(GMTdoc)grdblend.html#w)
- **Z** | **scale** :: [Type => Number]

    Scale output values by scale before writing to file.
    ($(GMTdoc)grdblend.html#z)
- $(GMT.opt_V)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_r)
"""
function grdblend(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:G :I :R :V_params :f :n :r])
	cmd  = parse_these_opts(cmd, d, [[:C :clobber], [:N :nodata], [:Q :headless], [:W :no_blend], [:Z :scale]])

	cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2)
	return common_grd(d, "grdblend " * cmd, arg1, arg2)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdblend(arg1, arg2=nothing, cmd0::String=""; kw...) = grdblend(cmd0, arg1, arg2; kw...)


text:
using AliasedSeabed
using Test
using SimradEK60

filename = joinpath(dirname(@__FILE__), "data/JR280_-D20121202-T170437.raw")

ps = SimradEK60.load(filename)
ps38 = [p for p in ps if p.frequency == 38000]
_Sv = Sv(ps38) # Volume backscatter
ntheta = alongshipangle(ps38) # Split beam angle
nphi = athwartshipangle(ps38)

_m1 = m1(ntheta)
_m2 = m2(nphi)
a = blackwell_asbmask(_Sv, ntheta, nphi)

@test sum(_m1) == 66165
@test sum(_m2) == 39929
@test sum(a) == 178399


text:
# using Pkg; Pkg.add("SineFit")

begin
    using Printf
    using DSP, SineFit
    using Plots, StatsPlots
end

struct SinregResult
    period::Float64
    params::SineFit.WaveFitParams
    plot
end

"""Sinusoidal regression function.

Fits a sinusoid to data provided as as `x` and `y` and returns a named tuple
containing both the fitted curve and the parameters of that curve.
Used to produce initial values for modeling data windows and later?|
to find peaks in modeled julian day values to align the result to
a cumulative age timeline.

[x] Vector of `x` values of input data

[y] Vector of `y` values of input data

[fixed_period] Optional variable for fixing the period of the sinusoid
in the depth domain. Defaults to `NA`, period is not fixed. Supply a
single value to fix the period.

[plot] Should the fitting result be plotted? `TRUE/FALSE`

RETURNS: A list containing a vector of parameters of the fitted sinusoid
and the fitted values belonging to each `x` value.
Fitting parameters:

`I` | sineparams.vertoffst = the mean annual value of the sinusoid (height)
`A` | sineparams.amplitude = the amplitude of the sinusoid
`Dper` | sineparams.frequency = the period of the sinusoid in `x` domain
`peak` | sineparams.phaseshft = the location of the peak in the sinusoid
`R2adj` = the adjusted `R^2` value of the fit
`p` = the p-value of the fit

R Examples
```
# Create dummy data
x <- seq(1000, 11000, 1000)
y <- sin((2 * pi * (seq(1, 11, 1) - 8 + 7 / 4)) / 7)
sinlist <- sinreg(x, y, plot = FALSE) # Run the function
@export
```
"""
function sinreg(
        x::Union{AbstractArray{Number, 1}, AbstractVector},
        y::Union{AbstractArray{Number, 1}, AbstractVector};
        fixed_period=nothing, plot::Bool=false, verbose::Bool=false
    )

    # Determine sine period
    if fixed_period == nothing
        # Get periodogram from FFT
        ssp = periodogram(y)
        # Get frequency of maximum density
        per = 1 / ssp.freq[ssp.power .== maximum(ssp.power)]
        println(per)
        dper = per * diff([extrema(x)...])
    else
        dper = fixed_period
    end

    # SineFit
    sineparams = calculate_wave_shape(x, y)
    sinpred = SineFit.sin(x, sineparams)
    fstring = @sprintf(
                "y = %.2f sin(%.2f x + %.2f) + %.2f",
                sineparams.amplitude, sineparams.frequency,
                sineparams.phaseshft, sineparams.vertoffst
                )

    # Print results
    if verbose
        print(fstring)
    end

    # Draw plot
    if plot
        # Plot results
        scatter(x, y, label="Data")
        plot = plot!(x, sinpred, lw=2, label="SinFit")
        annotate!(plot, [(-12, -0.8, (fstring, font("serif", 14, :red, :left,))),])
    else
        plot = nothing
    end

    # Gather results
    results = SinregResult(dper, sineparams, plot)

    # return results
    return results
end

x =  collect(-10π:0.01:10π)
y =  Base.sin.(x .* 2.1).+ rand(length(x))

result = sinreg(x, y, plot=true)

result.plot


text:
using Distributed # addprocs, @everywhere
using BenchmarkTools # @btime
using UnicodePlots # barplot on terminal
using Statistics # mean
using Base.Threads
using PiSMC

# checks if we're running Julia in multithread mode
PiSMC.Multithread.check_julia_num_threads()

qty_workers = 8
addprocs(qty_workers)

@everywhere begin
    using Pkg
    Pkg.activate(joinpath(@__DIR__, ".."))

    using PiSMC

    # Warmup
    PiSMC.run_naive(2)
    PiSMC.run_singlethread(2)
    PiSMC.Multithread.run_multithread(2)
    PiSMC.run_distributed(2)
end

@info("Started $(length(workers())) workers")
@info("Threads available: $(nthreads())")

num_sims = 100_000_000

let
    naive_result = PiSMC.run_naive(num_sims)
    singlethread_result = PiSMC.run_singlethread(num_sims)
    multithread_result = PiSMC.Multithread.run_multithread(num_sims)
    distributed_result = PiSMC.run_distributed(num_sims)

    @info("Naive result: $naive_result")
    @info("Singlethread result: $singlethread_result")
    @info("Multithread result: $multithread_result")
    @info("Distributed result: $distributed_result")
end

@info("Benchmarking naive...")
naive_bench = @benchmark PiSMC.run_naive(num_sims)
@info("Benchmarking singlethread...")
singlethread_bench = @benchmark PiSMC.run_singlethread(num_sims)
@info("Benchmarking multithread...")
multithread_bench = @benchmark PiSMC.Multithread.run_multithread(num_sims)
@info("Benchmarking distributed...")
distributed_bench = @benchmark PiSMC.run_distributed(num_sims)

plt = barplot( [ "naive", "singlethread", "multithread ($(nthreads()) threads)", "Distributed ($qty_workers workers)"],
              [ round(mean(v.times)/1E9, digits=4) for v in [ naive_bench, singlethread_bench, multithread_bench, distributed_bench ]],
        title="Benchmark (seconds)")

show(plt)


text:
"""
    VaR(lat::Lattice,alpha::Float32)
implements the nested version of the Value-at-Risk at level ``\\alpha`` for the lattice ``lat``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function VaR(lat::Lattice,alpha::Float32)
    # look for smallest x such that  P(-states <= x) >alpha , i.e.
    T0 = deepcopy(lat);

    for i=(length(T0.state)-1):-1:1
		for j=1:length(T0.state[i])
			T0.state[i][j] = RiskMeasures.VaR(T0.state[i+1],T0.probability[i+1][j,:,1],alpha)
		end
	end
    return(T0.state[1],T0)
end

"""
    CTE(lat::Lattice,alpha::Float32)
implements the nested version of the Conditional Tail Expectation level ``\\alpha`` for the lattice ``lat``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function CTE(lat::Lattice,alpha::Float32)
    T0 = deepcopy(lat);

	for i=(length(T0.state)-1):-1:1
		for j=1:length(T0.state[i])
			T0.state[i][j] = RiskMeasures.CTE(T0.state[i+1],T0.probability[i+1][j,:,1],alpha)
		end
	end

    return(T0.state[1],T0)
end

"""
    EVaR2(lat::Lattice,beta::Float32)
implements the nested version of the Entropic Value-at-Risk at level ``\\beta`` for the lattice ``lat``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version. Here the optimization is done with a simple
goldenSearch algorithm.
"""
function EVaR2(lat::Lattice,beta::Float32)
	T0 = deepcopy(lat);

	for i=(length(T0.state)-1):-1:1
		for j=1:length(T0.state[i])
			T0.state[i][j] = RiskMeasures.EVaR2(T0.state[i+1],T0.probability[i+1][j,:,1],beta)[1]
		end
	end

    return(T0.state[1],T0)
end

"""
    EVaR(lat::Lattice,beta::Float32)
implements the nested version of the Entropic Value-at-Risk at level ``\\beta`` for the lattice ``lat``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version. The optimization routine uses JuMP and Ipopt.
"""
function EVaR(lat::Lattice,beta::Float32)
	T0 = deepcopy(lat);

    for i=(length(T0.state)-1):-1:1
		for j=1:length(T0.state[i])
			T0.state[i][j] = RiskMeasures.EVaR(T0.state[i+1],T0.probability[i+1][j,:,1],beta)[1]
		end
	end
    return(T0.state[1],T0)
end

"""
    AVaR(lat::Lattice,alpha::Float32)
implements the nested version of the Average Value-at-Risk at level ``\\alpha`` for the lattice ``lat``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function AVaR(lat::Lattice, alpha::Float32)
    T0 = deepcopy(lat);

    for i=(length(T0.state)-1):-1:1
		for j=1:length(T0.state[i])
			T0.state[i][j] = RiskMeasures.AVaR(T0.state[i+1],T0.probability[i+1][j,:,1],alpha)[1]
		end
	end
    return(T0.state[1],T0)
end


text:
using TensorBoardLogger, Logging
using TensorBoardLogger: preprocess, summary_impl
using Test

test_log_dir = "test_logs/"

@testset "Initialization" begin
    isdir(test_log_dir) && rm(test_log_dir, force=true, recursive=true)

    # Check tb_append
    TBLogger(test_log_dir*"run")
    @test isdir(test_log_dir*"run")
    TBLogger(test_log_dir*"run") #tb_append by default
    @test isdir(test_log_dir*"run_1")
    TBLogger(test_log_dir*"run", tb_increment)
    @test isdir(test_log_dir*"run_2")

    # check tb_overwrite
    close(open(test_log_dir*"run_2/testfile", "w"))
    TBLogger(test_log_dir*"run_2", tb_overwrite)
    @test !isfile(test_log_dir*"run_2/testfile")
end

@testset "create log directory" begin
    import TensorBoardLogger.init_logdir
    rm(test_log_dir, recursive=true, force=true)

    # Check tb_append
    init_logdir(test_log_dir*"run")
    @test isdir(test_log_dir*"run")
    init_logdir(test_log_dir*"run") #tb_append by default
    @test isdir(test_log_dir*"run_1")
    init_logdir(test_log_dir*"run", tb_increment)
    @test isdir(test_log_dir*"run_2")

    # check tb_overwrite
    close(open(test_log_dir*"run_2/testfile", "w"))
    init_logdir(test_log_dir*"run_2", tb_overwrite)
    @test !isfile(test_log_dir*"run_2/testfile")
end

@testset "create event file" begin
    import TensorBoardLogger.init_logdir, TensorBoardLogger.create_eventfile
    rm(test_log_dir, recursive=true, force=true)
    logdir = init_logdir(test_log_dir*"run")

    fname, file = create_eventfile(logdir)
    @test isfile(joinpath(logdir, fname))

    fname, file = create_eventfile(logdir, prepend="test/")
    @test isfile(joinpath(logdir, fname))
end

@testset "stepping" begin
    # workaround step being already exported in julia 1.0
    #import TensorBoardLogger.step
    tb_step = TensorBoardLogger.step

    tbl = TBLogger(test_log_dir*"run", tb_overwrite)
    @test tb_step(tbl) == 0

    tbl = TBLogger(test_log_dir*"run", purge_step=12)
    @test tb_step(tbl) == 12

    @test increment_step!(tbl, 1) == 13
    @test tb_step(tbl) == 13

    @test set_step!(tbl, 1) == 1
    @test tb_step(tbl) == 1

    @test set_step_increment!(tbl, 10) == 10
    # changing the increment should not affect the step
    @test tb_step(tbl) == 1
end

@testset "resetting" begin
    tbl = TBLogger(test_log_dir*"run", tb_overwrite)
    TensorBoardLogger.add_eventfile(tbl, "pp")
    set_step!(tbl, 100)

    reset!(tbl)
    @test TensorBoardLogger.step(tbl) == 0
    @test length(tbl.all_files) == 1

end


text:
module AutoObjectsReStacker

"""
    restack(x) -> x

An identity function in the sense `restack(x) === x`.  However, it
(recursively) re-construct `x` to beg the compiler to move everything
in the stack.
"""
restack

_getlength(::Type{<:NTuple{N,Any}}) where N = N
_getnames(::Type{<:NamedTuple{names}}) where names = names

# Not using dispatch, to check `issingletontype` first.
@generated function restack(x)
    if Base.issingletontype(x)
        ex = x.instance :: x
    elseif x <: Tuple
        N = _getlength(x)
        ex = :(($(map(i -> :(restack(x[$i])), 1:N)...),))
    elseif x <: NamedTuple
        names = _getnames(x)
        ex = :((; $(map(n -> Expr(:kw, n, :(restack(x.$n))), names)...)))
    elseif parentmodule(x) === Core
        # Workaround the issue with `isstructtype(String)` and
        #  `isstructtype(Symbol)`.
        # https://github.com/JuliaLang/julia/issues/30210
        ex = :x
    elseif isstructtype(x)
        fn = fieldnames(x)
        lastfield = QuoteNode(fn[end])
        new = Expr(:new, x, map(n -> :(restack(getfield(x, $(QuoteNode(n))))), fn)...)

        ex = quote
            if isimmutable(x) && isdefined(x, $lastfield)
                $new
            else
                x
            end
        end
    else
        ex = :x
    end
    return Expr(:block, Expr(:meta, :inline), ex)
end

end  # module


text:
# ---
# title: 287. Find the Duplicate Number
# id: problem287
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers, Binary Search
# link: <https://leetcode.com/problems/find-the-duplicate-number/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` containing `n + 1` integers where each
# integer is in the range `[1, n]` inclusive.
#
# There is only **one duplicate number** in `nums`, return _this  duplicate
# number_.
#
# **Follow-ups:**
#
#   1. How can we prove that at least one duplicate number must exist in `nums`?
#   2. Can you solve the problem **without**  modifying the array `nums`?
#   3. Can you solve the problem using only constant, `O(1)` extra space?
#   4. Can you solve the problem with runtime complexity less than `O(n2)`?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,4,2,2]
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,1,3,4,2]
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,1,2]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 3 * 104`
#   * `nums.length == n + 1`
#   * `1 <= nums[i] <= n`
#   * All the integers in `nums` appear only **once** except for **precisely one integer** which appears **two or more** times.
#
#
## @lc code=start
using LeetCode

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


text:
function available(domain::GroundDomain, state::State)
    iters = ((t for (t, a) in group.actions if available(domain, state, a))
             for (_, group) in domain.actions)
    return Iterators.flatten(iters)
end

function available(domain::GroundDomain, state::State,
                   group::GroundActionGroup, args)
    term = Compound(group.name, args)
    return available(domain, state, group.actions[term])
end

function available(domain::GroundDomain, state::State, term::Term)
    if term.name == get_name(no_op) return true end
    if (term isa Const) term = Compound(term.name, []) end
    action = domain.actions[term.name].actions[term]
    return available(domain, state, action)
end

function available(domain::GroundDomain, state::State, action::GroundAction)
    return satisfy(domain, state, action.preconds)
end


text:
struct CancelToken
    cancelled::Threads.Atomic{Bool}
    restarts::Threads.Atomic{Int}
    cond::Threads.Condition
end

CancelToken() = CancelToken(
    Threads.Atomic{Bool}(false),
    Threads.Atomic{Int}(0),
    Threads.Condition()
)

function Base.close(token::CancelToken)
    lock(token.cond) do
        token.cancelled[] = true
        notify(token.cond)
        notify(Revise.revision_event);
    end
end

Base.isopen(token::CancelToken) = lock(() -> !token.cancelled[], token.cond)

function Base.wait(token::CancelToken)
    try
        lock(() -> wait(token.cond), token.cond)
    catch e
        if e isa InterruptException
            close(token)
        else
            rethrow()
        end
    end
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using Test, StatsBase, SearchModels, KCenters, KNearestCenters, TextSearch, TextClassification, CategoricalArrays
using Downloads, Random, MLDataUtils, JSON3, GZip

Random.seed!(1)
function train_test_split(corpus, labels, at=0.7)
    n = length(corpus)
    iX = shuffle!(collect(1:n))
    n_ = floor(Int, n * at)
    itrain, itest = iX[1:n_], iX[n_+1:end]
    corpus[itrain], labels[itrain], corpus[itest], labels[itest]
end

function folds_split(corpus, labels, folds=3)
    n = length(corpus)
    indexes = shuffle!(collect(1:n))
    folds = kfolds(indexes, folds)

    [(corpus[itrain], labels[itrain], corpus[itest], labels[itest]) for (itrain, itest) in folds]
end

@testset "microtc" begin
    !isfile("emo50k.json.gz") && Downloads.download("https://github.com/sadit/TextClassificationTutorial/raw/main/data/emo50k.json.gz", "emo50k.json.gz")
    labels = []
    corpus = []
    targets = ("♡", "💔")
    GZip.open("emo50k.json.gz") do f
        for line in eachline(f)
            tweet = JSON3.read(line)
            label = tweet["klass"]
            if label in targets
                push!(labels, tweet["klass"])
                push!(corpus, tweet["text"])
            end
        end
    end
    labels = categorical(labels)
    traincorpus, trainlabels, testcorpus, testlabels = train_test_split(corpus, labels, 0.7)
    folds = folds_split(traincorpus, trainlabels)

    function error_function(config::MicroTC_Config)
        S = Float64[]
        for (_traincorpus, _trainlabels, _testcorpus, _testlabels) in folds
            tc = MicroTC(config, _traincorpus, _trainlabels; verbose=true)
            #valX = vectorize_corpus(tc, _testcorpus)
            #ypred = predict.(tc, valX)
            ypred = predict_corpus(tc, _testcorpus)
            push!(S, recall_score(_testlabels, ypred, weight=:macro))
        end

        1.0 - mean(S)
    end

    for t in traincorpus[1:10]
        @show t
    end

    @show countmap(trainlabels)
    @show countmap(testlabels)

    space = MicroTC_ConfigSpace(
       textconfig=TextConfigSpace(
            qlist=[[4], [3], [5]],
            nlist=[[1], [1, 2], []],
            slist=[]
        ),
        # cls=KncConfigSpace(
        #     centerselection=[TextCentroidSelection()],
        #     kernel=[k_(AngleDistance()) for k_ in [DirectKernel]]
        # )
        #cls = KnnClassifierConfigSpace()
        #cls=LiblinearConfigSpace()
    )

    best_list = search_models(space, error_function, 32;
        maxpopulation=8,
        # search hyper-parameters
        bsize=2, mutbsize=8, crossbsize=8,
        tol=0.0, maxiters=30, verbose=true)

    for (i, b) in enumerate(best_list)
        @info "-- perf best_lists[$i]:", i, b[1], b[2]
    end

    cls = MicroTC(best_list[1][1], traincorpus, trainlabels)
    sc = classification_scores(testlabels, predict_corpus(cls, testcorpus))
    @info "*** Performance on test: " sc
    @test sc.accuracy > 0.6
end

# flush(stdout); flush(stderr)
# sort!(P, :score, rev=true)
# plot(P.score, label="", xlabel="rank", ylabel="score", markershape=:auto, title="RS performance on validation set") |> display


text:
using Printf

function exerciselogic(a::Bool, b::Bool)
    st = @sprintf " %5s" a
    st *= @sprintf " %5s" b
    st *= @sprintf " %5s" ~a
    st *= @sprintf " %5s" a | b
    st *= @sprintf " %5s" a & b
    st *= @sprintf " %5s" a $ b
end

println("Julia's logical operations on Bool:")
println("   a     b    not   or    and   xor")
for a in [true, false], b in [true, false]
    println(exerciselogic(a, b))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Printf

"""
    pretty_filesize(s, suffix="B")

Convert a floating point value `s` representing a file size to a more human-friendly
formatted string with one decimal places with a `suffix` defaulting to "B". Depending on
the value of `s` the string will be formatted to show `s` using an SI prefix from bytes,
kiB (1024 bytes), MiB (1024² bytes), and so on up to YiB (1024⁸ bytes).
"""
function pretty_filesize(s, suffix="B")
    # Modified from: https://stackoverflow.com/a/1094933
    for unit in ["", "Ki", "Mi", "Gi", "Ti", "Pi", "Ei", "Zi"]
        abs(s) < 1024 && return @sprintf("%3.1f %s%s", s, unit, suffix)
        s /= 1024
    end
    return @sprintf("%.1f %s%s", s, "Yi", suffix)
end


text:
using DSP #for phase unwrapping

struct AMFMcompN
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  a::Vector{Float64}
  ω::Vector{Float64}
  s::Vector{Float64}
  σ::Vector{Float64}
  θ::Vector{Float64}
  t::Vector{Float64}
  fs::Float64
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
"""
Determines how to react to an exceeded rate limit.

- `ORL_RETURN`: Return a [`Result`](@ref) containing a [`RateLimited`](@ref) exception.
- `ORL_WAIT`: Block and wait for the rate limit to expire.
"""
@enum OnRateLimit ORL_RETURN ORL_WAIT

"""
A signal that a rate limit has been exceeded.
Can contain the amount of time until its expiry.
"""
struct RateLimited <: Exception
    period::Union{Period, Nothing}
end

"""
A generic rate limiter using the `[X-]RateLimit-Remaining` and `[X-]RateLimit-Reset` response headers.
The reset header is assumed to be a Unix timestamp in seconds.
"""
mutable struct RateLimiter
    remaining::Int
    reset::Int

    RateLimiter() = new(1, 0)
end

rate_limit_check(rl::RateLimiter) =
    rl.remaining == 0 && rate_limit_period(rl) > Millisecond(0)

rate_limit_wait(rl::RateLimiter) = sleep(rate_limit_period(rl))

rate_limit_period(rl::RateLimiter) = max(Millisecond(0), unix2datetime(rl.reset) - now(UTC))

function rate_limit_update!(rl::RateLimiter, r::HTTP.Response)
    remaining = tryheader(i -> parse(Int, i), r, "RateLimit-Remaining")
    reset = tryheader(i -> parse(Int, i), r, "RateLimit-Reset")
    if remaining === nothing  || reset === nothing
        @warn "Parsing rate limit headers failed"
    else
        rl.remaining = remaining
        rl.reset = reset
    end
end

function tryheader(f::Function, r::HTTP.Response, header::AStr)
    for h in [header, "X-$header"]
        HTTP.hasheader(r, h) && return f(HTTP.header(r, h))
    end
    return nothing
end


text:
#Run this in Escher to see a rotating cube in the browser.

import ThreeJS

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    eventloop = every(1/60)
    rx = 0.0
    ry = 0.0
    rz = 0.0
    map(eventloop) do _
        rx += 0.5
        ry += 0.5
        rz += 0.5
        ThreeJS.outerdiv() <<
                (ThreeJS.initscene() <<
                    [
                        ThreeJS.mesh(0.0, 0.0, 0.0; rx=rx, ry=ry, rz=rz) <<
                        [
                            ThreeJS.box(5.0, 5.0, 5.0), ThreeJS.material(Dict(:kind=>"phong",:color=>"red"))
                        ],
                        ThreeJS.pointlight(150.0, 150.0, 150.0),
                        ThreeJS.pointlight(-150.0, -150.0, -150.0),
                        ThreeJS.camera(0.0, 0.0, 250.0)
                    ]
                )
        end
end


text:
# Dirichlet
struct Dirichlet{P<:Tuple{Vararg{Real}}} <: Distribution{Tuple{Vararg{Real}},P}
    params :: P
    Dirichlet(; α) = new{typeof(α)}(α)
end

# Dirichlet(α::Real...) = Dirichlet(α)
# Dirichlet(α::Integer...) = Dirichlet(Float64.(α))
# Dirichlet(n::Integer, α::Real)= Dirichlet(ntuple(_->α, n))

function randDirichlet(α::Tuple{Vararg{T}}) where T<:Real
    p = Tuple(randGamma(αi, one(T)) for αi in α)
    p ./ sum(p)
end
randDirichlet(α::Tuple{Vararg{Integer}}) = randDirichlet(Float64.(α))
randDirichlet(n::Integer, α::Real) = randDirichlet(Tuple(α for _ in 1:n))
randDirichlet(n::Integer, α::Integer) = randDirichlet(n, Float64(α))

@cufunc function logpdfDirichlet(x::NTuple{N,Real}, α::NTuple{N,Real}) where N
    # a, b = sum(u -> SVector(u,lgamma(u)), α)
    # s = sum(((u,v) -> (u-one(u))log(v)).(α, x))
    # s - b + lgamma(a)
    T = promote_type(typeof.(x)..., typeof.(α)...)
    a = b = s = zero(T)
    for (xi, αi) in zip(x, α)
        a += T(αi)
        b += lgamma(T(αi))
        s += (T(αi) - one(T))log(T(xi))
    end
    s - b + lgamma(a)
end


text:
#unum-subtraction.jl
#Performs subtraction with unums.  Requires two unums to have the same
#environment signature.

doc"""
  `Unums.frac_sub!(carry, subtrahend::Unum, minuend, guardbit::UInt64)`
  subtracts fraction from the fraction value of unum.
"""
function frac_sub!{ESS,FSS}(carry::UInt64, subtrahend::UnumSmall{ESS,FSS}, minuend::UInt64)
  (carry, subtrahend.fraction) = i64sub(carry, subtrahend.fraction, minuend)
  return carry
end
function frac_sub!{ESS,FSS}(carry::UInt64, subtrahend::UnumLarge{ESS,FSS}, minuend::ArrayNum{FSS})
  i64sub!(carry, subtrahend.fraction, minuend)
end

doc"""
  `Unums.sub(::Unum, ::Unum)` outputs a Unum OR Ubound corresponding to the difference
  of two unums.  This is bound to the (-) operation if options[:usegnum] is not
  set.  Note that in the case of degenerate unums, sub may change the bit values
  of the individual unums, but the values will not be altered.
"""
@universal function sub(a::Unum, b::Unum)
  #some basic checks out of the gate.
  (is_nan(a) || is_nan(b)) && return nan(U)
  is_zero(a) && return additiveinverse!(copy(b))
  is_zero(b) && return copy(a)

  #resolve degenerate conditions in both A and B before calculating the exponents.
  resolve_degenerates!(a)
  resolve_degenerates!(b)

  #go ahead and decode the a and b exponents, these will be used, a lot.
  _aexp = decode_exp(a)
  _bexp = decode_exp(b)

  #check to see if the signs on a and b are mismatched.
  if ((a.flags $ b.flags) & UNUM_SIGN_MASK) != z16
    (_aexp >= _bexp) ? unum_sum(a, b, _aexp, _bexp) : additiveinverse!(unum_sum(b, a, _bexp, _aexp))
  else
    is_inward(b, a) ? unum_diff(a, b, _aexp, _bexp) : additiveinverse!(unum_diff(b, a, _bexp, _aexp))
  end
end

#import the Base add operation and bind it to the add and add! functions
import Base.-
@bind_operation(-, sub)

doc"""
  `Unums.unum_diff(::Unum, ::Unum, _aexp, _bexp)` outputs a Unum OR Ubound
  corresponding to the difference of two unums.  This function as a prerequisite
  must have the b be inward of a
"""
@universal function unum_diff(a::Unum, b::Unum, _aexp::Int64, _bexp::Int64)
  #basic secondary checks which eject early results.
  is_inf(a) && return is_inf(b) ? nan!(zero(U)) : copy(a)
  is_inf_ulp(a) && return inf_ulp_sub(a, b)
  #there is a corner case that b winds up being infinity (and a does not; same
  #with mmr.)

  if (is_exact(a) && is_exact(b))
    diff_exact(a, b, _aexp, _bexp)
  else
    diff_inexact(a, b, _aexp, _bexp)
  end
end

@universal function diff_exact(a::Unum, b::Unum, _aexp::Int64, _bexp::Int64)
  #track once whether or not a and b are subnormal
  _a_subnormal = is_exp_zero(a)
  _b_subnormal = is_exp_zero(b)

  #modify the exponent such that they are.
  _aexp += _a_subnormal * 1
  _bexp += _b_subnormal * 1

  #calculate the shift between _aexp and _bexp.
  _shift = to16(_aexp - _bexp)

  if _shift > max_fsize(FSS)
    #then go down one previous exact unum and decrement.
    return inner_ulp!(make_exact!(copy(a)))
  end

  #copy the b unum as the temporary result, nuke its ubit just to be sure
  #(we might be calling this function from diff_inexact)
  result = make_exact!(copy(b))
  #set the sign to the sign of the dominant figure.
  coerce_sign!(result, a)

  if _shift == z16
    carry::UInt64 = ((!_a_subnormal) & _b_subnormal) * o64
    guardbit::Bool = false
  else

    carry = (!_a_subnormal) * o64
    guardbit = get_bit(result.fraction, (max_fsize(FSS) + o16) - _shift)
    rsh_and_set_ubit!(result, _shift, true)
    (_b_subnormal) || frac_set_bit!(result, _shift)
  end

  #subtract fractionals parts together, and reset the carry.
  carry = frac_sub!(carry, result, a.fraction)

  if (carry == z64)

    #nb we only need to check the top bit.
    is_frac_zero(result) && return zero!(result, (guardbit != z64) * UNUM_UBIT_MASK)

    _mexp = min_exponent(ESS)

    if (_aexp != _mexp)
      #next count how many zeros are in the front of the fraction.
      places_to_shift = clz(result.fraction) + o16

      #check to make sure that places_to_shift doesn't get too big.
      (_aexp - places_to_shift < _mexp) && (places_to_shift = to16(_aexp - _mexp))

      frac_lsh!(result, places_to_shift)

      if (guardbit)
        result.fraction -= 1
      end
      (result.esize, result.exponent) = encode_exp(_aexp - places_to_shift)
    end
  else
    #check to see if we need to the guard bit to set ubit.
    (guardbit != z64) && make_ulp!(result)
    #set the exponent
    (result.esize, result.exponent) = encode_exp(_aexp)
  end

  is_exact(result) && exact_trim!(result)
  trim_and_set_ubit!(result)

  return result
end

@universal function diff_inexact(a::Unum, b::Unum, _aexp::Int64, _bexp::Int64)
  #do a second is_inward check.  If this is_inward check fails, then the result
  #can have an opposite sign, because it's an ulp that goes wierdly.

  a_out, a_in = is_exact(a) ? (a, a) : (outer_exact(a), inner_exact(a))
  b_out, b_in = is_exact(b) ? (b, b) : (outer_exact(b), inner_exact(b))

  #println("a:")
  #print("out:"); describe(a_out)
  #print("in:");  describe(a_in)
  #println("b:")
  #print("out:"); describe(b_out)
  #print("in:");  describe(b_in)

  if is_inf(a_out)
    first_value = a_out
  elseif is_zero(a_out)
    first_value = b_in
  elseif is_zero(b_in)
    first_value = a_out
  else
    first_value = is_inward(b_in, a_out) ?
      diff_exact(a_out, b_in, decode_exp(a_out), decode_exp(b_in)) :
      diff_exact(b_in, a_out, decode_exp(b_in), decode_exp(a_out))
  end

  if is_inf(b_out)
    second_value = b_out
  elseif is_zero(a_in)
    second_value = b_out
  elseif is_zero(b_out)
    second_value = a_in
  else
    second_value = is_inward(b_out, a_in) ?
      diff_exact(a_in, b_out, decode_exp(a_in), decode_exp(b_out)) :
      diff_exact(b_out, a_in, decode_exp(b_out), decode_exp(a_in))
  end

  outer_hull_lub = lub(first_value) > lub(second_value) ? lub(first_value) : lub(second_value)
  outer_hull_glb = glb(first_value) < glb(second_value) ? glb(first_value) : glb(second_value)

  outer_hull_top = is_zero(outer_hull_lub) ? neg_sss(U) : lower_ulp(outer_hull_lub)
  outer_hull_bot = is_zero(outer_hull_glb) ? pos_sss(U) : upper_ulp(outer_hull_glb)

  #next find the outer hull of these values, and return it.
  resolve_as_utype!(outer_hull_bot, outer_hull_top)
end

@universal function inf_ulp_sub(a::Unum, b::Unum)
  is_inf_ulp(b) && return B(cheapest_inf_ulp(U, UNUM_SIGN_MASK), cheapest_inf_ulp(U)) #all real numbers
  #calculate inner_value by doing an exact subtraction of b from big_exact.
  #ensure that it's an ulp.
  diff_bound = is_ulp(b) ? outer_exact(b) : b
  inner_value = diff_exact(a, diff_bound, max_exponent(ESS), decode_exp(diff_bound))
  #coerce sign and ulp to cover zero and exact cases.
  outer_ulp!(make_exact!(coerce_sign!(inner_value, a)))

  #check for the dominant sign.  If a was positive make it (inner_value -> mmr)
  #if a was negative, make it (-mmr -> inner_value )
  if ((@signof a) == z16)
    resolve_as_utype!(inner_value, pos_mmr(U))
  else
    resolve_as_utype!(neg_mmr(U), inner_value)
  end
end

@universal function subtract_ubit!(x::Unum, s::UInt16)
  borrowed = frac_sub_ubit!(x, s)
  if borrowed
    exponent = decode_exp(x)
    if exponent <= min_exponent(ESS)
      x.exponent = z64
    else
      (x.esize, x.exponent) = encode_exp(exponent - 1)
      frac_rsh!(x, 1)
    end
  end
  return make_ulp!(x)
end

doc"""
  `Unums.diff_overlap(::Unum, ::Unum, ::Int64, ::Int64)`
  calculates the ubound that results when two unums have overlapping ulps.
"""
@universal function diff_overlap(a::Unum, b::Unum, _aexp, _bexp)
  #it's possible there is a smarter algorithm than this.  For now, though, it'll do.

  #first, decide which ulp unum is dominant.
  (d, e, _dexp, _eexp) = (a.fsize < b.fsize) ? (a, b, _aexp, _bexp) : (b, a, _bexp, _aexp)

  d_shell = coerce_sign!(outer_exact(d), a)
  e_shell = coerce_sign!(outer_exact(e), a)
  top = inner_ulp!(make_exact!(diff_exact(d_shell, e, decode_exp(d_shell), _eexp)))
  bot = inner_ulp!(additiveinverse!(make_exact!(diff_exact(e_shell, d, decode_exp(e_shell), _dexp))))
  return is_positive(a) ? B(bot, top) : B(top, bot)
end

#unary subtraction creates a new unum and flips it.
@universal function -(x::Unum)
  additiveinverse!(copy(x))
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module CellValuesBench

using Gridap.FieldValues
using Gridap.CellValues

include("../test/CellValuesTests/Mocks.jl")

function doloop(x)
  for xi in x
  end
end

l = 1000000

println("+++ CellValuesBench ( length = $l ) +++")

sv = 1.0
scv = ConstantCellValue(sv,l)
sv2 = 1.1
scv2 = ConstantCellValue(sv2,l)
sa = [sv, sv, sv]
sca = ConstantCellValue(sa,l)
sa2 = [sv sv; sv sv; sv sv]
sca2 = ConstantCellValue(sa2,l)

print("ConstantCellValue ->"); @time doloop(scv)
print("ConstantCellValue ->"); @time doloop(scv)

print("ConstantCellValue ->"); @time doloop(sca2)
print("ConstantCellValue ->"); @time doloop(sca2)

sv = 1.0
scv = TestCellValue(sv,l)
sv2 = 1.1
scv2 = TestCellValue(sv2,l)
sa = [sv, sv, sv]
sca = TestCellArray(sa,l)
sa2 = [sv sv; sv sv; sv sv]
sca2 = TestCellArray(sa2,l)

print("TestCellValue ->"); @time doloop(scv)
print("TestCellValue ->"); @time doloop(scv)

print("TestCellArray ->"); @time doloop(sca2)
print("TestCellArray ->"); @time doloop(sca2)

for op in (:+,:-)
  @eval begin
    scv3 = $op(scv)
    print("CellValueUnary($(string($op))) ->"); @time doloop(scv3)
    print("CellValueUnary($(string($op))) ->"); @time doloop(scv3)
  end
end

for op in (:+,:-)
  @eval begin
    scv3 = $op(sca)
    print("CellArrayUnary($(string($op))) ->"); @time doloop(scv3)
    print("CellArrayUnary($(string($op))) ->"); @time doloop(scv3)
  end
end

sca3 = cellsum(sca2,dim=2)
print("CellArrayFromCellSum ->"); @time doloop(sca3)
print("CellArrayFromCellSum ->"); @time doloop(sca3)

sca3 = cellsum(sca,dim=1)
print("CellValueFromCellSum ->"); @time doloop(sca3)
print("CellValueFromCellSum ->"); @time doloop(sca3)

sca3 = cellnewaxis(sca2,dim=2)
print("CellArrayFromNewAxis ->"); @time doloop(sca3)
print("CellArrayFromNewAxis ->"); @time doloop(sca3)

for op in (:+,:-,:*,:/,:(outer),:(inner))
  @eval begin
    scv3 = $op(scv,scv2)
    print("CellValueBinary($(string($op))) ->"); @time doloop(scv3)
    print("CellValueBinary($(string($op))) ->"); @time doloop(scv3)
  end
end

for op in (:+,:-,:*,:/,:(outer),:(inner))
  @eval begin
    sca3 = $op(sca,sca2)
    print("CellArrayBinary($(string($op))) ->"); @time doloop(sca3)
    print("CellArrayBinary($(string($op))) ->"); @time doloop(sca3)
  end
end

for op in (:+,:-,:*,:/,:(outer),:(inner))
  @eval begin
    sca3 = $op(scv,sca2)
    print("CellValueArrayBinary($(string($op))) ->"); @time doloop(sca3)
    print("CellValueArrayBinary($(string($op))) ->"); @time doloop(sca3)
  end
end

end # module CellValuesBench


text:
module InfoTheory

using Distances
using NearestNeighbors

export differential_entropy
export IntConfusionMatrix, count
#export NCDMultiset, ncd1, delta_diameter
export compress, clen
export ncd
export ncd1_sequence, ncdm, ncd1plot
export ncdmatrix

include("knn_entropy.jl")
include("confusion_matrix.jl")

include("utils.jl")

# NCD-related functions
include("utils/isinstalled.jl")
include("compress.jl")
include("ncd.jl")
include("ncd_matrix.jl")
#include("ncdmultiset.jl")

end # module


text:
using DistributedFactorGraphs
using IncrementalInference

using Test

##

mutable struct MutableLinearConditional{N, T <: SamplableBelief} <: AbstractRelativeRoots
    Z::T
    # timestamp::DateTime
end

function MutableLinearConditional{N}() where N
    newval = MvNormal(zeros(N), diagm(ones(N)))
    MutableLinearConditional{N,typeof(newval)}(newval)
end
MutableLinearConditional(n::Int=1) = MutableLinearConditional{n}()
MutableLinearConditional(nm::Distributions.ContinuousUnivariateDistribution) = MutableLinearConditional{1, typeof(nm)}(nm)
MutableLinearConditional(nm::MvNormal) = MutableLinearConditional{length(nm.μ), typeof(nm)}(nm)
MutableLinearConditional(nm::BallTreeDensity) = MutableLinearConditional{Ndim(nm), typeof(nm)}(nm)

getDimension(::Type{MutableLinearConditional{N,<:SamplableBelief}}) where {N} = N
getManifolds(::Type{MutableLinearConditional{N,<:SamplableBelief}}) where {N} = tuple([:Euclid for i in 1:N]...)

IIF.getSample(cf::CalcFactor{<:MutableLinearConditional}, N::Int=1) = (reshape(rand(cf.factor.Z,N),:,N), )
function (s::CalcFactor{<:MutableLinearConditional})(meas,
                                                     X1,
                                                     X2  )
#
    return meas .- (X2 .- X1)
end

@testset "testing dead reckoning tether" begin

# test error message and then define method for MutableLinearConditional
@test_throws ErrorException getFactorMean(MutableLinearConditional(Normal(0.0,0.1)))

IIF.getFactorMean(fct::MutableLinearConditional) = getFactorMean(fct.Z)
##

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, ContinuousScalar)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], Prior( Normal(0.0,0.1) ))

# Drive around in line
for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, ContinuousScalar)
    pp = LinearRelative(Normal(1.0,0.1))
    addFactor!(fg, [psym;nsym], pp )
end

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, ContinuousScalar, tags=[:LANDMARK])
p2br = LinearRelative(Normal(1.0,0.1))
addFactor!(fg, [:x0; :l1], p2br )

## async solving with dead-reckon branch

addVariable!(fg, :deadreckon_x0, ContinuousScalar, solvable=0)

drec = MutableLinearConditional(Normal(0.0,0.1))

addFactor!(fg, [:x0; :deadreckon_x0], drec, solvable=0)

#
@test length(map( x->x.label, getVariables(fg, solvable=1))) == 8
@test length(map( x->x.label, getVariables(fg, solvable=0))) == 9
#
# # make sure
@test length(getEliminationOrder(fg, solvable=1)) == 8
# check default
@test length(getEliminationOrder(fg)) == 8

# default check
vo = getEliminationOrder(fg)
@test length(vo) == 8

tree = buildTreeFromOrdering!(fg,vo)

tree2, smt, hists = solveTree!(fg);

@test !isInitialized(fg, :deadreckon_x0)

val = accumulateFactorMeans(fg, [:x0deadreckon_x0f1])

@test norm(val - calcVariablePPE(fg, :x0).suggested) < 1e-8

#TODO improve test
rebaseFactorVariable!(fg, :x0deadreckon_x0f1, [:x1; :deadreckon_x0])
@test issetequal(ls2(fg, :x0), [:x1, :l1])
@test issetequal(ls2(fg, :x1), [:x0, :deadreckon_x0, :x2])

end


text:
@bitmask_flag FlagBits::UInt32 begin
    BIT_A = 1
    BIT_B = 2
    BIT_C = 3
end

@bitmask_flag FlagBits2::UInt32 begin
    BIT_A_2 = 1
    BIT_B_2 = 2
    BIT_C_2 = 3
end

@testset "Vulkan enum wrappers" begin
    @test BIT_A & BIT_B == FlagBits(0)
    @test BIT_A | BIT_B == BIT_C
    @test xor(BIT_A, BIT_B) == BIT_C
    @test BIT_A < BIT_B && BIT_B < BIT_C && BIT_A < BIT_C
    @test BIT_A & 1 == 1 & BIT_A == BIT_A
    @test BIT_A | 2 == 2 | BIT_A == BIT_C
    @test xor(BIT_A, 2) == xor(2, BIT_A) == BIT_C
    @test 1 < BIT_B && BIT_B < 3
    @test BIT_A <= BIT_B <= BIT_C
    @test BIT_A == BIT_A
    @test BIT_B in BIT_C
    @test BIT_A ∉ BIT_B
    @test 2 in BIT_C
    @test BIT_B in 3
    @test_throws ErrorException BIT_A & BIT_A_2
    @test_throws ErrorException BIT_A == BIT_A_2
    @test_throws ErrorException BIT_A in BIT_A_2
end