text:
using Pythia

pythia = newPythia()

icxx"""$pythia->readString("Beams:idA = 2212");"""
icxx"""$pythia->readString("Beams:idB = 2212");"""
icxx"""$pythia->readString("Beams:eCM = 8000.");"""
icxx"""$pythia->readString("HardQCD:gg2qqbar = on");"""
icxx"""$pythia->init();"""

ρ⁺ = Array(Float64,0)
ρ⁻ = Array(Float64,0)
ρ⁰ = Array(Float64,0)
for e in take(events(pythia),5)
    for p in particles(e)
        id = icxx"$p.id();"
        m = icxx"$p.m();"
        if id == 213
            push!(ρ⁺,m)
        elseif id == -213
            push!(ρ⁻,m)
        elseif id == 113
            push!(ρ⁰,m)
        end
    end
end

particles = [ρ⁺; ρ⁻; ρ⁰]
colors = [
    [:ρ⁺ for _ in 1:length(ρ⁺)];
    [:ρ⁻ for _ in 1:length(ρ⁻)];
    [:ρ⁰ for _ in 1:length(ρ⁰)]]
plot(y = particles, color = colors, Geom.point)

function track_decay(it,id)
    p = it[id]
    daughter1 = icxx"$p.daughter1();"
    daughter2 = icxx"$p.daughter2();"
    if daughter1 == daughter2 == 0
        println("$id ($(name(it[id]))) does not decay")
        return
    end
    for d in daughter1:daughter2
        println("$id ($(name(it[id]))) decays to $d ($(name(it[d])))")
        track_decay(it,d)
    end
end

for e in take(events(pythia),5)
    ps = particles(e)
    for (i,p) in enumerate(ps)
        # Look for Delta0 or 3112
        id = icxx"$p.id();"
        (2114 == id || 3112 == id) && track_decay(ps,i-1)
    end
    println()
end


text:
struct Identity <: LocalOperator
end

kernelvals(biop::Identity, x) = nothing
integrand(op::Identity, kernel, x, g, f) = dot(f[1], g[1])
scalartype(op::Identity) = Union{}

struct NCross <: LocalOperator
end

kernelvals(op::NCross, mp) = nothing
integrand(op::NCross, kernel, x, g, f) = dot(g[1], normal(x) × f[1])
scalartype(op::NCross) = Union{}

function _alloc_workspace(qd, g, f, tels, bels)
    q = qd[1]
    τ = tels[1]
    w, p = q[1], neighborhood(τ,q[2])
    a = (w, p, g(p), f(p))
    A = Vector{typeof(a)}(undef,length(qd))
end

const LinearRefSpaceTriangle = Union{RTRefSpace, NDRefSpace}
function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels)
    u, w = trgauss(6)
    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)

    return qd, A
end

function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels)
    u, w = trgauss(6)
    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

const LinearRefSpaceTetr = Union{NDLCCRefSpace, NDLCDRefSpace}
function quaddata(op::LocalOperator, g::LinearRefSpaceTetr, f::LinearRefSpaceTetr, tels, bels)
    o, x, y, z = CompScienceMeshes.euclidianbasis(3)
    reftet = simplex(x,y,z,o)
    qps = quadpoints(reftet, 3)
    qd = [(w, parametric(p)) for (p,w) in qps]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector)

    u, w = legendre(6, 0.0, 1.0)
    qd = [(w[i],u[i]) for i in eachindex(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector)

    u, w = trgauss(6)
    qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector)

     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
     reftet = simplex(x,y,z,o)
     qps = quadpoints(reftet, 6)
     qd = [(w, parametric(p)) for (p,w) in qps]
     A = _alloc_workspace(qd, g, f, tels, bels)
     return qd, A
end

function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A))
    for i in eachindex(qd)
        q = qd[i]
        w, p = q[1], neighborhood(τ,q[2])
        A[i] = (w, p, ψ(p), ϕ(p))
    end
    return A
end

ext: jl


text:
using Documenter, PredictableComponents

makedocs(
    modules = [CoupledFields,PredictableComponents],
    clean = false,
    format = :html,
    sitename = "Predictable\nComponents.jl",
    pages = [
        "Home" => "index.md",
        "Guide" => Any["man/Example 1.md","man/Example 2.md"],
        "Library" => "lib/library.md"
    ],
    # Use clean URLs, unless built as a "local" build
    html_prettyurls = !("local" in ARGS),
#    html_canonical = "https://juliadocs.github.io/Documenter.jl/stable/",
)

deploydocs(
    repo = "github.com/Mattriks/PredictableComponents.jl.git",
    target = "build",
    julia = "0.6",
    osname = "linux",
    deps = nothing,
    make = nothing,
)


text:
using ParameterEstimocean.InverseProblems: vectorize, forward_run!, transpose_model_output
using Oceananigans.Architectures: arch_array
using CairoMakie
using LaTeXStrings

include("visualize_profile_predictions_utils.jl")

# Temporary hack -- couldn't get Makie to render LaTeX properly
guide = ["⁰","¹","²","³","⁴","⁵","⁶","⁷","⁸","⁹"]

"""
    visualize!(ip::InverseProblem, parameters;
                    field_names = [:u, :v, :b, :e],
                    directory = pwd(),
                    filename = "realizations.png"
                    )

    For visualizing 1-dimensional time series predictions.
"""
function visualize!(ip::InverseProblem, parameters;
                    field_names = [:u, :v, :b, :e],
                    directory = pwd(),
                    filename = "realizations.png"
                    )

    isdir(directory) || mkdir(directory)

    model = ip.simulation.model

    n_fields = length(field_names)

    observations = vectorize(ip.observations)

    forward_run!(ip, parameters)

    # Vector of SyntheticObservations objects, one for each observation
    predictions = transpose_model_output(ip.time_series_collector, observations)

    fig = Figure(resolution = (200*(length(field_names)+1), 200*(length(ip.observations)+1)), font = "CMU Serif")
    colors = [:black, :red, :blue]

    function empty_plot!(fig_position)
        ax = fig_position = Axis(fig_position)
        hidedecorations!(ax)
        hidespines!(ax, :t, :b, :l, :r)
    end

    for (oi, observation) in enumerate(observations)

        i = oi + 1
        prediction = predictions[oi]

        targets = observation.times
        snapshots = round.(Int, range(1, length(targets), length=3))

        Qᵇ = arch_array(CPU(), model.tracers.b.boundary_conditions.top.condition)[1,oi]
        Qᵘ = arch_array(CPU(), model.velocities.u.boundary_conditions.top.condition)[1,oi]
        fv = arch_array(CPU(), model.coriolis)[1,oi].f

        empty_plot!(fig[i,1])
        text!(fig[i,1], "Qᵇ = $(tostring(Qᵇ)) m⁻¹s⁻³\nQᵘ = $(tostring(Qᵘ)) m⁻¹s⁻²\nf = $(tostring(fv)) s⁻¹",
                    position = (0, 0),
                    align = (:center, :center),
                    textsize = 15,
                    justification = :left)

        for (j, field_name) in enumerate(field_names)

            middle = j > 1 && j < n_fields
            remove_spines = j == 1 ? (:t, :r) : j == n_fields ? (:t, :l) : (:t, :l, :r)
            axis_args = j == n_fields ? (ylabelposition=:right, yaxisposition=:right) : NamedTuple()

            if j == 1 || j == n_fields
                axis_args = merge(axis_args, (ylabel="z (m)",))
            end

            j += 1 # reserve the first column for row labels

            info = field_guide[field_name]

            grid = observation.grid

            z = field_name ∈ [:u, :v] ? grid.zᵃᵃᶠ[1:grid.Nz] : grid.zᵃᵃᶜ[1:grid.Nz]

            to_plot = field_name ∈ keys(prediction.field_time_serieses)

            if to_plot

                # field data for each time step
                obsn = getproperty(observation.field_time_serieses, field_name)
                pred = getproperty(prediction.field_time_serieses, field_name)
                obsn = view(arch_array(CPU(), obsn.data), 1,1,1:grid.Nz,:)

                z_offset = pred.data.offsets[3]
                z_start = 1 - z_offset
                z_end = grid.Nz - z_offset

                pred = view(arch_array(CPU(), parent(parent(pred.data))), 1,oi,z_start:z_end,:)

                om = order_of_magnitude(maximum(abs.(obsn)))
                scaling = 10^(-om)
                om_string = string([getindex(guide, parse(Int64, c)+1) for c in string(abs(om))]...)
                om_string = om < 0 ? "⁻" * om_string : om_string
                xlabel = string("$(info.name) (10$(om_string) $(info.units))")

                ax = Axis(fig[i,j]; xlabelpadding=0,
                                    xtickalign=1,
                                    ytickalign=1,
                                    xlabel,
                                    xticks=LinearTicks(3),
                                    axis_args...)

                hidespines!(ax, remove_spines...)

                middle && hideydecorations!(ax, grid=false)

                lins = []
                for (color_index, t) in enumerate(snapshots)
                    l = lines!([obsn[:,t] .* scaling...], z; color = (colors[color_index], 0.4))
                    push!(lins, l)
                    l = lines!([pred[:,t] .* scaling...], z; color = (colors[color_index], 1.0), linestyle = :dash)
                    push!(lins, l)
                end

                times = @. round((targets[snapshots] - targets[snapshots[1]]) / 86400, sigdigits=2)

                legendlabel(time) = ["Observation, t = $time days", "Model, t = $time days"]
                Legend(fig[1,2:3], lins, vcat([legendlabel(time) for time in times]...), nbanks=2)
                lins = []
            else

                empty_plot!(fig[i,j])
            end
        end
    end

    save(joinpath(directory, filename), fig, px_per_unit = 2.0)
    return nothing
end

# visualize!(ip::InverseProblem, parameters; kwargs...) = visualize!(model_time_series(ip, parameters); kwargs...)

# function visualize_and_save!(calibration, validation, parameters, directory; fields=[:u, :v, :b, :e])
#         isdir(directory) || makedir(directory)

#         path = joinpath(directory, "results.txt")
#         o = open_output_file(path)
#         write(o, "Training relative weights: $(calibration.relative_weights) \n")
#         write(o, "Validation relative weights: $(validation.relative_weights) \n")
#         write(o, "Training default parameters: $(validation.default_parameters) \n")
#         write(o, "Validation default parameters: $(validation.default_parameters) \n")

#         write(o, "------------ \n \n")
#         default_parameters = calibration.default_parameters
#         train_loss_default = calibration(default_parameters)
#         valid_loss_default = validation(default_parameters)
#         write(o, "Default parameters: $(default_parameters) \nLoss on training: $(train_loss_default) \nLoss on validation: $(valid_loss_default) \n------------ \n \n")

#         train_loss = calibration(parameters)
#         valid_loss = validation(parameters)
#         write(o, "Parameters: $(parameters) \nLoss on training: $(train_loss) \nLoss on validation: $(valid_loss) \n------------ \n \n")

#         write(o, "Training loss reduction: $(train_loss/train_loss_default) \n")
#         write(o, "Validation loss reduction: $(valid_loss/valid_loss_default) \n")
#         close(o)

#         parameters = calibration.parameters.ParametersToOptimize(parameters)

#         for inverse_problem in [calibration, validation]

#             all_data = inverse_problem.observations
#             simulation = inverse_problem.simulation
#             set!(simulation.model, parameters)

#             for data_length in Set(length.(getproperty.(all_data, :t)))

#                 observations = [d for d in all_data if length(d.t) == data_length]
#                 days = observations[1].t[end]/86400

#                 new_ip = InverseProblem()

#                 visualize!(simulation, observations, parameters;
#                             fields = fields,
#                             filename = joinpath(directory, "$(days)_day_simulations.png"))
#             end
#         end

#     return nothing
# end


text:
function grapfOutput!(data::alpheusData)
    x = data.preprocessor.x
    H = data.output.H
    U = data.output.U
    J = data.input.J
    upstream = data.options.upstream
    dowstream = data.options.downstream
    solver = data.options.solver
    z = data.output.z

    plot(x, H ,title="Altura",xlabel="X (m)", ylabel="H (m)", label="Altura (m)",legend=:topleft)
    scatter!(x[z.==1.0], H[z.==1.0], color = "green", label = "Turbina")
    savefig("altura_$(solver)_$(upstream)_$(dowstream)_$(J).png")

    plot(x, U ,title="Velocidade",xlabel="X (m) ", ylabel="U (m/s)", label="Velocidade (m/s)")
    scatter!(x[z.==1.0], U[z.==1.0], color = "green", label = "Turbina")
    savefig("velocidade_$(solver)_$(upstream)_$(dowstream)_$(J).png")

    nothing
end

function textOutput!(data::alpheusData)
    x = data.preprocessor.x
    H = data.output.H
    U = data.output.U
    W = data.output.W
    J = data.input.J
    upstream = data.options.upstream
    dowstream = data.options.downstream
    solver = data.options.solver
    z = data.output.z

    open("output_$(J)","w") do file
        write(file, "$(J) $(data.preprocessor.x'z) $(H'z) $(U'z) $(W'z)")
    end

    nothing
end

function grapfOutputEnum!(data::alpheusData)
    x = data.preprocessor.x
    W = data.output.W_enum
    J = data.input.J
    upstream = data.options.upstream
    dowstream = data.options.downstream
    solver = data.options.solver

    W = W*g*density*data.input.Q

    plot(x, W ,title="Potência",xlabel="X (m)", ylabel="Y (W)", label="Potência(W)",legend=:topright)
    savefig("energiaEnum_$(solver)_$(upstream)_$(dowstream)_$(J).png")

    nothing
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
function is_diag_dom(Q::SparseMatrixCSC{Float64,Int64})
    for i = 1:size(Q,2)
      if 3 * Q[i,i] < sum(Q[:,i]) + sum(Q[i,:])
        return false
      end
    end

    return true
end

function dist_diag_dom(Q::SparseMatrixCSC{Float64,Int64})
    delta_dist = 0.0
    for i = 1:size(Q,2)
      delta_dist = max(delta_dist,sum(Q[:,i]) + sum(Q[i,:]) - 3 * Q[i,i])
    end

    return delta_dist
end

#function bisection_search_delta!(iter::Class_iterate, kkt_solver::abstract_KKT_system_solver, min_delta::Float64, max_delta::Float64,  pars::Class_parameters, timer::class_advanced_timer)
#max_exponent = ceil(log(max_delta/min_delta))

function int_bisection(f::Function, a::Int64, b::Int64)
    for i = 1:10000
        if a == b - 1
          return b
        end
        midpoint = ceil((a + b)/2.0)
        if f(midpoint) == true
          b = midpoint
        else
          a = midpoint
        end
    end
end

function ipopt_strategy!(iter::Class_iterate, kkt_solver::abstract_KKT_system_solver, pars::Class_parameters, timer::class_advanced_timer)
    step_info = Blank_ls_info()

    MAX_IT = 500
    DELTA_ZERO = pars.delta.zero
    #DELTA_MIN = iter.point.mu * LinearAlgebra.norm(iter.point.y,Inf) / 100.0 #
    DELTA_MIN = pars.delta.min
    DELTA_MAX = pars.delta.max

    #Q = get_lag_hess(iter)
    #@show dist_diag_dom(Q)

    num_fac = 0

    inertia = 0
    status = :none
    tau = 1.5 * diag_min(kkt_solver)
    #Q = Symmetric(kkt_solver.M,:L)
    #eigvals,vec, = eigs(Q,nev=1,which=:SR,maxiter=10,tol=1e4)
    #@show eigvals
    #v = randn(length(iter.point.x))
    #v = v /LinearAlgebra.norm(v,2)
    #@show dot(v,Q * v)
    delta = DELTA_ZERO
    # see if we can succeed with delta = DELTA_ZERO
    if tau > 0.0
      tau = 0.0
      inertia = factor!( kkt_solver, delta, timer )
      num_fac += 1
      if inertia == 1
        return :success, num_fac, delta
      end
    end
    if pars.output_level >= 4
        println(pd("delta"), pd("inertia"))
    end
    for i = 1:MAX_IT
        if pars.output_level >= 4
          println("delta=",rd(delta), "inertia=", pd(inertia))
        end

        if i == 1
          if get_delta(iter) != 0.0
            delta = max(DELTA_MIN - tau, get_delta(iter) * pars.delta.dec)
          else
            delta = choose_delta_start(iter, kkt_solver, pars) - tau
          end
        else
            delta = delta * pars.delta.inc
        end

        inertia = factor!( kkt_solver, delta, timer )
        num_fac += 1

        n = length(iter.point.x)

        if inertia == 1
          return :success, num_fac, delta
        elseif is_diag_dom(kkt_solver.Q[1:n,1:n])
            println("WARNING: Inertia calculation incorrect")
            @warn("Inertia calculation incorrect")
        end

        if delta > DELTA_MAX
          dx = LinearAlgebra.norm(kkt_solver.dir.x,Inf)
          dy = LinearAlgebra.norm(kkt_solver.dir.y,Inf)
          ds = LinearAlgebra.norm(kkt_solver.dir.s,Inf)

          my_warn("ipopt_strategy failed with delta_max=$DELTA_MAX, delta=$delta, i=$i")
          my_warn("num_fac=$num_fac, inertia=$inertia, status=$status, dir_x=$dx, dir_y=$dx, dir_s=$ds")

          #error("ipopt_strategy failed with too big a delta")
          return :failure, iter, delta
        end
    end

    error("max it")
end

function choose_delta_start(iter::Class_iterate, kkt_solver::abstract_KKT_system_solver, pars::Class_parameters)
    return pars.delta.start

    # do something clever to reduce number of factorizations

end


text:
# Discounted HEDEGEd KOM

"""
    dHedgePPM{SymbolType}( c::Int, [ β::Float64 = 1.0, γ::Float64 = 1.0 ] )

Creates a Discounted HEDGE Predictor for `SymbolType` with context length `c`, learning
parameter ``β`` and discounting parameter ``γ``.

## Examples
```julia-repl
julia> p = dHedgePPM{Char}( 4 )
DiscretePredictors.dHedgePPM{Char}([*] (0)
, Char[], 4, 1.0, 1.0, [1.0, 1.0, 1.0, 1.0, 1.0], Dict{Char,Float64}[Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}()])

julia> p = dHedgePPM{Char}( 4, 0.8 )
DiscretePredictors.dHedgePPM{Char}([*] (0)
, Char[], 4, 0.8, 1.0, [1.0, 1.0, 1.0, 1.0, 1.0], Dict{Char,Float64}[Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}()])

julia> p = dHedgePPM{Char}( 4, 0.9, 0.8 )
DiscretePredictors.dHedgePPM{Char}([*] (0)
, Char[], 4, 0.9, 0.8, [1.0, 1.0, 1.0, 1.0, 1.0], Dict{Char,Float64}[Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}()])
```

Reference:
Raj, Vishnu, and Sheetal Kalyani. "An aggregating strategy for shifting experts in discrete sequence prediction." arXiv preprint arXiv:1708.01744 (2017).
"""
mutable struct dHedgePPM{T} <: BasePredictor{T}
    model::Trie{T,Int}
    context::Vector{T}
    cxt_length::Int
    β::Float64
    γ::Float64
    weights::Vector{Float64}
    last_prediction::Vector{Dict{T,Float64}}

    dHedgePPM{T}( _c::Int, _β::Float64 = 1.0, _γ::Float64 = 1.0 ) where {T} = new( Trie{T,Int}(), Vector{T}(),
                                                _c, _β, _γ,
                                                ones(Float64,_c+1), fill(Dict{T,Float64}(),_c+1) )
end

function add!( p::dHedgePPM{T}, sym::T ) where {T}
    p.model.value += 1;

    # Update weights
    for expert_idx = 1:p.cxt_length+1
        expert_best_symbol  = get_best_symbol( p.last_prediction[expert_idx] )
        expert_loss         = (expert_best_symbol==sym) ? 0 : 1
        p.weights[expert_idx] = (p.weights[expert_idx]^p.γ) * (p.β^expert_loss)
    end

    # Create node path
    push!( p.context, sym );
    buffer  = p.context[1:end];
    while !isempty(buffer)
        p.model[buffer] = (haskey(p.model,buffer)) ? p.model[buffer]+1 : 1
        popfirst!( buffer )
    end

    if length(p.context) > p.cxt_length
        popfirst!( p.context )
    end

    nothing
end

function predict( p::dHedgePPM{T} ) where {T}
    # Clear away last predictions
    p.last_prediction   = fill(Dict{T,Float64}(),p.cxt_length+1)

    # Calculate normalized weights
    norm_weights    = p.weights ./ sum(p.weights)
    # println( "    norm_weights: ", norm_weights )

    # Get prediction from 0-order expert
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    # Save this for later
    p.last_prediction[1]    = deepcopy( symbols )

    # Apply 0-order expert weight
    for sym ∈ keys(symbols)
        symbols[sym] *= norm_weights[1]
    end
    # println( "After Order-0: ", symbols );

    # Do prediction from all other higher orders
    no_of_active_experts    = min( p.model.value, p.cxt_length )
    for expert_idx  = 1:no_of_active_experts
        # Predict
        this_expert_prediction = predict_from_subcontext( p, p.context[1:expert_idx] )
        # Normalize probability out from experts - jugaad
        denom   = 0;
        for (sym,prob) ∈ this_expert_prediction denom += prob; end
        for sym ∈ keys(this_expert_prediction) this_expert_prediction[sym] /= denom;    end

        # Save
        p.last_prediction[expert_idx+1] = deepcopy( this_expert_prediction )
        # Apply weight to expert and add to final prediction
        for sym ∈ keys( symbols )
            symbols[sym] += (this_expert_prediction[sym]*norm_weights[expert_idx+1])
        end
        # println( "After order-$expert_idx: ", symbols )
    end

    return symbols;
end

function info_string( p::dHedgePPM{T} ) where {T}
    return @sprintf( "dHedgePPM(%d,\$\\beta\$ = %3.2f,\$\\gamma\$ = %3.2f)",
                            p.cxt_length, p.β, p.γ )
end

function unique_string( p::dHedgePPM{T} ) where {T}
    return @sprintf( "dHedgePPM_%02d_%03d_%03d", p.cxt_length,
                        trunc(Int,100*p.β), trunc(Int,100*p.γ)  )
end

function predict_from_subcontext( p::dHedgePPM{T}, sub_cxt::Vector{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    buffer  = Vector{T}()
    for i = length(sub_cxt):-1:1
        pushfirst!( buffer, sub_cxt[i] )
        # println( "    Buffer: ", buffer, " -> ", p.model[buffer] );
        # Apply escape probability
        esc_prob = 1/p.model[buffer]
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        list_of_children    = children( p.model, buffer )
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            symbols[k] += list_of_children[k].value/p.model[buffer]
        end
    end

    return symbols
end


text:
Base.eachcol(tdf::TimeDataFrame) = eachcol(dataframe(tdf))
Base.eachrow(tdf::TimeDataFrame) = eachrow(dataframe(tdf))
ncol(tdf::TimeDataFrame) = DataFrames.ncol(dataframe(tdf))
nrow(tdf::TimeDataFrame) = DataFrames.nrow(dataframe(tdf))


text:
# Copyright 2016, Wen-chih Chen, Yueh-shan Chung
# JuMP4DEA
# A Modeling Language Package for Solving Data Envelopment Analysis Problems

module JuMP4DEA

# An algebraic modeling language for Julia
import JuMP
# Solver-independent functions and low-level interface for Mathematical Programming
import MathProgBase

# function
export solveDEA

# solving dea problems by a smart algorithm
# it is computationally efficient and can solve large scale DEA problems via any solvers (even with LP size limitation)
include("solveDEA.jl")

end


text:
using Documenter, FlightPlot

makedocs(
    sitename = "FlightPlot.jl",
    modules = [FlightPlot],
    pages = Any[
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/rjdverbeek-tud/FlightPlot.jl.git",
)


text:
push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, DocumenterGHActionContainer

DocMeta.setdocmeta!(DocumenterGHActionContainer, :DocTestSetup, :(using DocumenterGHActionContainer), recursive = true)

makedocs(sitename = "DocumenterGHActionContainer",
         modules = [DocumenterGHActionContainer],
         pages = [
             "Home" => "index.md",
             "API" => "api.md"
         ]
)

deploydocs(
    repo   = "github.com/Nosferican/DocumenterGHActionContainer.jl.git",
    push_preview = true
)


text:
# Execution control

export
    launch, CuDim

typealias CuDim Union(Int, (Int, Int), (Int, Int, Int))
dim3(g::Int) = (g, 1, 1)
dim3(g::(Int, Int)) = (g[1], g[2], 1)
dim3(g::(Int, Int, Int)) = g

function launch(f::CuFunction, grid::CuDim, block::CuDim, args::Tuple;
                shmem_bytes::Int=4, stream::CuStream=default_stream())
	griddim = dim3(grid)
	blockdim = dim3(block)

    @assert all(dim->(dim > 0), grid)
    @assert all(dim->(dim > 0), block)

    @assert all([isbits(arg) || isa(arg, DevicePtr) for arg in args])
    kernel_args = [[arg] for arg in args]

    @cucall(:cuLaunchKernel, (
        Ptr{Void},  			# function
        Cuint, Cuint, Cuint,  	# grid dimensions (x, y, z)
        Cuint, Cuint, Cuint,  	# block dimensions (x, y, z)
        Cuint,  				# shared memory bytes,
        Ptr{Void}, 				# stream
        Ptr{Ptr{Void}}, 		# kernel parameters
        Ptr{Ptr{Void}}), 		# extra parameters
        f.handle, griddim[1], griddim[2], griddim[3], blockdim[1], blockdim[2],
        blockdim[3], shmem_bytes, stream.handle, kernel_args, C_NULL)
end


text:
export acqNumFGFrames, acqNumBGFrames, acqOffsetFieldShift, acqFramePeriod,
       acqNumPeriods, acqNumPatches, acqNumPeriodsPerPatch, acqFov,
       acqGradientDiag,
       rxNumFrequencies, rxFrequencies, rxTimePoints,
       measFGFrameIdx, measBGFrameIdx, measBGFrameBlockLengths

rxNumFrequencies(f::MPIFile) = floor(Int,rxNumSamplingPoints(f) ./ 2 .+ 1)

function rxFrequencies(f::MPIFile)
  numFreq = rxNumFrequencies(f)
  a = collect(0:(numFreq-1))./(numFreq-1).*rxBandwidth(f)
  return a
end

function rxTimePoints(f::MPIFile)
  numTP = rxNumSamplingPoints(f)
  a = collect(0:(numTP-1))./(numTP).*dfCycle(f)
  return a
end

function acqGradientDiag(f::MPIFile)
  g = acqGradient(f)
  g_ = reshape(g,9,size(g,3),size(g,4))
  return g_[[1,5,9],:,:]
end

function acqFov(f::MPIFile)
  if size(dfStrength(f)[1,:,:],1) == 3
    return  2*dfStrength(f)[1,:,:] ./ abs.( acqGradientDiag(f)[:,1,:] )
  else
    return  2*dfStrength(f)[1,:,:] ./ abs.( acqGradientDiag(f)[1,1,1] )
  end
end

acqFramePeriod(b::MPIFile) = dfCycle(b) * acqNumAverages(b) * acqNumPeriodsPerFrame(b)

# numPeriods is the total number of DF periods in a measurement.
acqNumPeriods(f::MPIFile) = acqNumFrames(f)*acqNumPeriodsPerFrame(f)

function acqOffsetFieldShift(f::MPIFile)
    return acqOffsetField(f) ./ reshape( acqGradient(f),9,1,:)[[1,5,9],:,:]
end

acqNumFGFrames(f::MPIFile) = acqNumFrames(f) - acqNumBGFrames(f)
acqNumBGFrames(f::MPIFile) = sum(measIsBGFrame(f))

function measBGFrameIdx(f::MPIFile)
  idx = zeros(Int64, acqNumBGFrames(f))
  j = 1
  mask = measIsBGFrame(f)
  for i=1:acqNumFrames(f)
    if mask[i]
      idx[j] = i
      j += 1
    end
  end
  return idx
end

function measFGFrameIdx(f::MPIFile)
  mask = measIsBGFrame(f)
  if !any(mask)
    #shortcut
    return 1:acqNumFrames(f)
  end
  idx = zeros(Int64, acqNumFGFrames(f))
  j = 1
  for i=1:acqNumFrames(f)
    if !mask[i]
      idx[j] = i
      j += 1
    end
  end
  return idx
end

function measBGFrameBlockLengths(mask)
  len = Vector{Int}(undef,0)

  groupIdxStart = -1
  for i=1:(length(mask)+1)
    if i <= length(mask) && mask[i] && groupIdxStart == -1
      groupIdxStart = i
    end
    if groupIdxStart != -1 && ((i == length(mask)+1) || !mask[i])
      push!(len, i-groupIdxStart)
      groupIdxStart = - 1
    end
  end
  return len
end

function acqNumPatches(f::MPIFile)
  # not valid for varying gradients / multi gradient
  shifts = acqOffsetFieldShift(f)
  return size(unique(shifts,dims=3),3)
end

function acqNumPeriodsPerPatch(f::MPIFile)
  return div(acqNumPeriodsPerFrame(f), acqNumPatches(f))
end

export unflattenOffsetFieldShift

unflattenOffsetFieldShift(f::MPIFile) = unflattenOffsetFieldShift(acqOffsetFieldShift(f))
function unflattenOffsetFieldShift(shifts::Array)
  # not valid for varying gradients / multi gradient
  uniqueShifts = unique(shifts, dims=2)
  numPeriodsPerFrame = size(shifts,2)
  numUniquePatch = size(uniqueShifts,2)

  allPeriods = 1:numPeriodsPerFrame

  flatIndices = Vector{Vector{Int64}}()

  for i=1:numUniquePatch
    temp = allPeriods[vec(sum(shifts .== uniqueShifts[:,i],dims=1)).==3]
    push!(flatIndices, temp)
  end
  return flatIndices
end

# We assume that systemMatrixWithBG has already reordered the BG data
# to the end
systemMatrix(f::MPIFile) = systemMatrixWithBG(f)[1:acqNumFGFrames(f),:,:,:]

function measDataTD(f, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),
                  receivers=1:rxNumChannels(f))

  data1 = measData(f,frames,periods,receivers)

  if measIsFastFrameAxis(f)
    data2 = permutedims(data1, invperm([4,1,2,3]))
  else
    data2 = data1
  end

  if measIsFourierTransformed(f)
    dataTD = irfft(data2,2*size(data2,1)-1, 1)
  else
    dataTD = data2
  end
  return dataTD
end


text:
module timeevolution_master

export master, master_nh, master_h, master_dynamic, master_nh_dynamic

import ..integrate, ..recast!

using ...bases, ...states, ...operators
using ...operators_dense, ...operators_sparse

const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void}

"""
    timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_h as derivative function.

Further information can be found at [`master`](@ref).
"""
function master_h(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                rates::DecayRates=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Void}=nothing,
                kwargs...)
    check_master(rho0, H, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_nh(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_nh as derivative function.

In this case the given Hamiltonian is assumed to be the non-hermitian version:
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
Further information can be found at [`master`](@ref).
"""
function master_nh(tspan, rho0::DenseOperator, Hnh::Operator, J::Vector;
                rates::DecayRates=nothing,
                Hnhdagger::Operator=dagger(Hnh),
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Void}=nothing,
                kwargs...)
    check_master(rho0, Hnh, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master(tspan, rho0, H, J; <keyword arguments>)

Time-evolution according to a master equation.

There are two implementations for integrating the master equation:

* [`master_h`](@ref): Usual formulation of the master equation.
* [`master_nh`](@ref): Variant with non-hermitian Hamiltonian.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
        operator type.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                rates::DecayRates=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Void}=nothing,
                kwargs...)
    isreducible = check_master(rho0, H, J, Jdagger, rates)
    if !isreducible
        tmp = copy(rho0)
        dmaster_h_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
        return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...)
    else
        Hnh = copy(H)
        if typeof(rates) == Matrix{Float64}
            for i=1:length(J), j=1:length(J)
                Hnh -= 0.5im*rates[i,j]*Jdagger[i]*J[j]
            end
        elseif typeof(rates) == Vector{Float64}
            for i=1:length(J)
                Hnh -= 0.5im*rates[i]*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= 0.5im*Jdagger[i]*J[i]
            end
        end
        Hnhdagger = dagger(Hnh)
        tmp = copy(rho0)
        dmaster_nh_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
        return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...)
    end
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic non-hermitian Hamiltonian and J.

In this case the given Hamiltonian is assumed to be the non-hermitian version.
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger)`
or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
at [`master_dynamic`](@ref).
"""
function master_nh_dynamic(tspan, rho0::DenseOperator, f::Function;
                rates::DecayRates=nothing,
                fout::Union{Function,Void}=nothing,
                kwargs...)
    tmp = copy(rho0)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic Hamiltonian and J.

There are two implementations for integrating the master equation with dynamic
operators:

* [`master_dynamic`](@ref): Usual formulation of the master equation.
* [`master_nh_dynamic`](@ref): Variant with non-hermitian Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `f`: Function `f(t, rho) -> (H, J, Jdagger)` or `f(t, rho) -> (H, J, Jdagger, rates)`
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, rho0::DenseOperator, f::Function;
                rates::DecayRates=nothing,
                fout::Union{Function,Void}=nothing,
                kwargs...)
    tmp = copy(rho0)
    dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

# Automatically convert Ket states to density operators
master(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master(tspan, dm(psi0), H, J; kwargs...)
master_h(tspan, psi0::Ket, H::Operator, J::Vector; kwargs...) = master_h(tspan, dm(psi0), H, J; kwargs...)
master_nh(tspan, psi0::Ket, Hnh::Operator, J::Vector; kwargs...) = master_nh(tspan, dm(psi0), Hnh, J; kwargs...)
master_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_dynamic(tspan, dm(psi0), f; kwargs...)
master_nh_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_nh_dynamic(tspan, dm(psi0), f; kwargs...)

# Recasting needed for the ODE solver is just providing the underlying data
function recast!(x::Vector{Complex128}, rho::DenseOperator)
    rho.data = reshape(x, size(rho.data))
end
recast!(rho::DenseOperator, x::Vector{Complex128}) = nothing

function integrate_master(tspan, df::Function, rho0::DenseOperator,
                        fout::Union{Void, Function}; kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    x0 = reshape(rho0.data, length(rho0))
    state = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data)
    integrate(tspan_, df, x0, state, dstate, fout; kwargs...)
end

# Time derivative functions
#   * dmaster_h
#   * dmaster_nh
#   * dmaster_h_dynamic -> callback(t, rho) -> dmaster_h
#   * dmaster_nh_dynamic -> callback(t, rho) -> dmaster_nh
# dmaster_h and dmaster_nh provide specialized implementations depending on
# the type of the given decay rate object which can either be nothing, a vector
# or a matrix.

function dmaster_h(rho::DenseOperator, H::Operator,
                    rates::Void, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, H, rho, 0, drho)
    operators.gemm!(1im, rho, H, 1, drho)
    for i=1:length(J)
        operators.gemm!(1, J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[i], 1, drho)

        operators.gemm!(-0.5, Jdagger[i], tmp, 1, drho)

        operators.gemm!(1., rho, Jdagger[i], 0, tmp)
        operators.gemm!(-0.5, tmp, J[i], 1, drho)
    end
    return drho
end

function dmaster_h(rho::DenseOperator, H::Operator,
                    rates::Vector{Float64}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, H, rho, 0, drho)
    operators.gemm!(1im, rho, H, 1, drho)
    for i=1:length(J)
        operators.gemm!(rates[i], J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[i], 1, drho)

        operators.gemm!(-0.5, Jdagger[i], tmp, 1, drho)

        operators.gemm!(rates[i], rho, Jdagger[i], 0, tmp)
        operators.gemm!(-0.5, tmp, J[i], 1, drho)
    end
    return drho
end

function dmaster_h(rho::DenseOperator, H::Operator,
                    rates::Matrix{Float64}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, H, rho, 0, drho)
    operators.gemm!(1im, rho, H, 1, drho)
    for j=1:length(J), i=1:length(J)
        operators.gemm!(rates[i,j], J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[j], 1, drho)

        operators.gemm!(-0.5, Jdagger[j], tmp, 1, drho)

        operators.gemm!(rates[i,j], rho, Jdagger[j], 0, tmp)
        operators.gemm!(-0.5, tmp, J[i], 1, drho)
    end
    return drho
end

function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
                    rates::Void, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, Hnh, rho, 0, drho)
    operators.gemm!(1im, rho, Hnh_dagger, 1, drho)
    for i=1:length(J)
        operators.gemm!(1, J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[i], 1, drho)
    end
    return drho
end

function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
                    rates::Vector{Float64}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, Hnh, rho, 0, drho)
    operators.gemm!(1im, rho, Hnh_dagger, 1, drho)
    for i=1:length(J)
        operators.gemm!(rates[i], J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[i], 1, drho)
    end
    return drho
end

function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
                    rates::Matrix{Float64}, J::Vector, Jdagger::Vector,
                    drho::DenseOperator, tmp::DenseOperator)
    operators.gemm!(-1im, Hnh, rho, 0, drho)
    operators.gemm!(1im, rho, Hnh_dagger, 1, drho)
    for j=1:length(J), i=1:length(J)
        operators.gemm!(rates[i,j], J[i], rho, 0, tmp)
        operators.gemm!(1, tmp, Jdagger[j], 1, drho)
    end
    return drho
end

function dmaster_h_dynamic(t::Float64, rho::DenseOperator, f::Function,
                    rates::DecayRates,
                    drho::DenseOperator, tmp::DenseOperator)
    result = f(t, rho)
    @assert 3 <= length(result) <= 4
    if length(result) == 3
        H, J, Jdagger = result
        rates_ = rates
    else
        H, J, Jdagger, rates_ = result
    end
    check_master(rho, H, J, Jdagger, rates_)
    dmaster_h(rho, H, rates_, J, Jdagger, drho, tmp)
end

function dmaster_nh_dynamic(t::Float64, rho::DenseOperator, f::Function,
                    rates::DecayRates,
                    drho::DenseOperator, tmp::DenseOperator)
    result = f(t, rho)
    @assert 4 <= length(result) <= 5
    if length(result) == 4
        Hnh, Hnh_dagger, J, Jdagger = result
        rates_ = rates
    else
        Hnh, Hnh_dagger, J, Jdagger, rates_ = result
    end
    check_master(rho, Hnh, J, Jdagger, rates_)
    dmaster_nh(rho, Hnh, Hnh_dagger, rates_, J, Jdagger, drho, tmp)
end

function check_master(rho0::DenseOperator, H::Operator, J::Vector, Jdagger::Vector, rates::DecayRates)
    isreducible = true # test if all operators are sparse or dense
    check_samebases(rho0, H)
    if !(isa(H, DenseOperator) || isa(H, SparseOperator))
        isreducible = false
    end
    for j=J
        @assert isa(j, Operator)
        if !(isa(j, DenseOperator) || isa(j, SparseOperator))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    for j=Jdagger
        @assert isa(j, Operator)
        if !(isa(j, DenseOperator) || isa(j, SparseOperator))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    @assert length(J)==length(Jdagger)
    if typeof(rates) == Matrix{Float64}
        @assert size(rates, 1) == size(rates, 2) == length(J)
    elseif typeof(rates) == Vector{Float64}
        @assert length(rates) == length(J)
    end
    isreducible
end

end #module


text:
abstract type MeanFunction end

"""
    ZeroMean{T<:Real} <: MeanFunction

Returns `zero(T)` everywhere.
"""
struct ZeroMean{T<:Real} <: MeanFunction end

"""
This is an AbstractGPs-internal workaround for AD issues; ideally we would just extend Base.map
"""
_map_meanfunction(::ZeroMean{T}, x::AbstractVector) where {T} = zeros(T, length(x))

function ChainRulesCore.rrule(::typeof(_map_meanfunction), m::ZeroMean, x::AbstractVector)
    map_ZeroMean_pullback(Δ) = (NoTangent(), NoTangent(), ZeroTangent())
    return _map_meanfunction(m, x), map_ZeroMean_pullback
end

ZeroMean() = ZeroMean{Float64}()

"""
    ConstMean{T<:Real} <: MeanFunction

Returns `c` everywhere.
"""
struct ConstMean{T<:Real} <: MeanFunction
    c::T
end

_map_meanfunction(m::ConstMean, x::AbstractVector) = fill(m.c, length(x))

"""
    CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy. Must be able to be mapped over an
`AbstractVector` of inputs.
"""
struct CustomMean{Tf} <: MeanFunction
    f::Tf
end

_map_meanfunction(f::CustomMean, x::AbstractVector) = map(f.f, x)


text:
using InterpolatedRejectionSampling
using InterpolatedRejectionSampling: get_knots, integrate, normalize_interp,
    is_normalized, get_interp, maxmapreduce, Cells, sample, Support, iterate,
    propose_sample, get_extrema, Envelope, rsample, irsample, irsample!

using Test
using Interpolations

using Random: seed!
seed!(1234)

slice = [missing,missing,0.3]
knots = (π.*[0.0, sort(rand(20))..., 1.0],
         π.*[0.0, sort(rand(19))..., 1.0],
         π.*[0.0, sort(rand(18))..., 1.0])

coefs = [sin(x)*sin(y)*sin(z) for x=knots[1],y=knots[2],z=knots[3]]
interp = LinearInterpolation(knots,coefs)

ranges = (range(0, stop=π, length=20),
          range(0, stop=π, length=19),
          range(0, stop=π, length=18))
coefs_ranges = [sin(x)*sin(y)*sin(z) for x=ranges[1],y=ranges[2],z=ranges[3]]
interp_ranges = LinearInterpolation(ranges,coefs_ranges)

@testset "Utilities" begin
    @test get_knots(interp) === knots
    @test isapprox(integrate(interp), 8.0; atol = 0.5)

    @test get_knots(interp_ranges) === ranges
    @test isapprox(integrate(interp_ranges), 8.0; atol = 0.5)

    ninterp = normalize_interp(interp)
    @test isa(ninterp, AbstractExtrapolation)
    @test isapprox(integrate(ninterp), one(Float64))
    @test is_normalized(ninterp)

    ninterp = normalize_interp(interp_ranges)
    @test isa(ninterp, AbstractExtrapolation)
    @test isapprox(integrate(ninterp), one(Float64))
    @test is_normalized(ninterp)

    pnt = (0.2,1.0,0.9)
    @test get_interp(interp, pnt) === interp(pnt...)
    @test get_interp(interp_ranges, pnt) === interp_ranges(pnt...)
    @test maxmapreduce(x -> x^2, [0.0, 0.5, -2.0, 1.5]) == 4.0
end

@testset "Cells" begin
    cells = Cells(interp)
    @test isa(cells, Cells)

    s = sample(cells)
    @test isa(s, CartesianIndex{3})

    cells = Cells(interp,slice)
    @test isa(cells, Cells)

    s = sample(cells)
    @test isa(s, CartesianIndex{3})

    cells = Cells(interp_ranges)
    @test isa(cells, Cells)
end

@testset "Support" begin
    cells = Cells(interp)
    s = sample(cells)
    supp = Support(cells, s)
    @test isa(supp, Support{Float64,3})

    @test isa(iterate(supp), Tuple{NTuple{2,Float64}, Int})
    @test isa(iterate(supp,2), Tuple{NTuple{2,Float64}, Int})
    @test isa(iterate(supp,3), Tuple{NTuple{2,Float64}, Int})
    @test isa(iterate(supp,4), Nothing)

    samp = propose_sample(supp)
    @test isa(samp, NTuple{3,Float64})

    spnts = get_extrema(supp)
    @test length(spnts) == 8
    @test eltype(spnts) == NTuple{3,Float64}
end

@testset "Envelope/rsample" begin
    cells = Cells(interp)
    s = sample(cells)
    envelope = Envelope(cells, s)
    @test isa(envelope, Envelope)
    @test isa(rsample(envelope), NTuple{3,Float64})
    @test isa(rsample(interp), NTuple{3,Float64})
    @test isa(rsample(interp, slice), NTuple{3,Float64})
end

@testset "irsample/irsample!" begin
    @test isa(irsample(knots, coefs), NTuple{3,Float64})

    samp = irsample(knots, coefs, 2)
    @test isa(samp, Matrix{Float64})
    @test size(samp) == (3,2)

    slices = Matrix{Union{Missing,Float64}}(missing, 3, 3)
    slices[1,1], slices[1,2], slices[2,2] = 0.2, 0.3, 0.4
    @test !iszero(count(ismissing, slices))

    irsample!(slices, knots, coefs)
    @test iszero(count(ismissing, slices))
    @test isa(convert(Matrix{Float64}, slices), Matrix{Float64})

    @test isa(irsample(ranges, coefs_ranges), NTuple{3,Float64})

    samp = irsample(ranges, coefs_ranges, 2)
    @test isa(samp, Matrix{Float64})
    @test size(samp) == (3,2)

    slices = Matrix{Union{Missing,Float64}}(missing, 3, 3)
    slices[1,1], slices[1,2], slices[2,2] = 0.2, 0.3, 0.4
    @test !iszero(count(ismissing, slices))

    irsample!(slices, ranges, coefs_ranges)
    @test iszero(count(ismissing, slices))
    @test isa(convert(Matrix{Float64}, slices), Matrix{Float64})
end


text:
module ICapeThermoPropertyRoutine

    export calcandgetlnphi!,calcsinglephaseprop,calctwophaseprop,checksinglephasepropspec!
    export checktwophasepropspec!,getsinglephaseproplist!,gettwophaseproplist!

    """
      This method is used to calculate the natural logarithm of the fugacity coefficients (and
      optionally their derivatives) in a single Phase mixture. The values of temperature, pressure
      and composition are specified in the argument list and the results are also returned through
      the argument list.
    """
    function calcandgetlnphi!(
        #= [in] =# phaseLabel::AbstractString,
        #= [in] =# temperature::Float64,
        #= [in] =# pressure::Float64,
        #= [in] =# moleNumbers::Vector{Float64},
        #= [in] =# fFlags::Int32,
        #= [out][in] =# lnPhi!::Vector{Float64},
        #= [out][in] =# lnPhiDT!::Vector{Float64},
        #= [out][in] =# lnPhiDP!::Vector{Float64},
        #= [out][in] =# lnPhiDn!::Vector{Float64})
		end

    """
      CalcSinglePhaseProp is used to calculate properties and property derivatives of a mixture in
      a single Phase at the current values of temperature, pressure and composition set in the
      Material Object. CalcSinglePhaseProp does not perform phase Equilibrium Calculations.
    """
    function calcsinglephaseprop(
        #= [in] =# props::Vector{AbstractString},
        #= [in] =# phaseLabel::AbstractString)
		end

    """
      CalcTwoPhaseProp is used to calculate mixture properties and property derivatives that depend on
      two Phases at the current values of temperature, pressure and composition set in the Material Object.
      It does not perform Equilibrium Calculations.
    """
    function calctwophaseprop(
        #= [in] =# props::Vector{AbstractString} ,
        #= [in] =# phaseLabels::Vector{AbstractString} )
		end

    """
      Checks whether it is possible to calculate a property with the CalcSinglePhaseProp method for a given Phase.
    """
    function checksinglephasepropspec!(
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabel::AbstractString,
        #= [retval][out] =# valid::Bool)
		end

    """
      Checks whether it is possible to calculate a property with the CalcTwoPhaseProp method for a given set of Phases.
    """
    function checktwophasepropspec!(
        #= [in] =# property::AbstractString,
        #= [in] =# phaseLabels::Vector{AbstractString},
        #= [retval][out] =# valid::Bool)
		end

    """
      Returns the list of supported non-constant single-phase Physical Properties.
    """
    function getsinglephaseproplist!(
        #= [retval][out] =# props::Vector{AbstractString})
		end

    """
      Returns the list of supported non-constant two-phase properties.
    """
    function gettwophaseproplist!(
        #= [retval][out] =# props::Vector{AbstractString})
		end

end


text:
import Images
using Base.Test, Color, FixedPointNumbers

gray = linspace(0.0,1.0,5)
ovr = Images.Overlay((2gray, 2gray), (RGB(1,0,0), RGB(0,0,1)), (Images.Clamp{Float64}(), Images.Clamp{Float64}()))
@test ovr[1] == RGB(0,0,0)
@test ovr[2] == RGB(0.5,0,0.5)
@test ovr[3] == ovr[4] == ovr[5] === RGB(1,0,1)
@test eltype(ovr) == RGB{Float64}
@test length(ovr) == 5
@test size(ovr) == (5,)
@test size(ovr,1) == 5
@test size(ovr,2) == 1
@test Images.nchannels(ovr) == 2
iob = IOBuffer()
show(iob, ovr)

ovr = Images.Overlay((gray, 0*gray), (RGB{Ufixed8}(1,0,1), RGB{Ufixed8}(0,1,0)), ((0,1),(0,1)))
@test eltype(ovr) == RGB{Ufixed8}
s = similar(ovr)
@test isa(s, Vector{RGB{Ufixed8}}) && length(s) == 5
s = similar(ovr, RGB{Float32})
@test isa(s, Vector{RGB{Float32}}) && length(s) == 5
s = similar(ovr, RGB{Float32}, (3,2))
@test isa(s, Matrix{RGB{Float32}}) && size(s) == (3,2)
buf = Images.uint32color(ovr)
gray8 = uint8(255*gray)
nogreen = [uint32(g)<<16 | uint32(g) for g in gray8]
@test buf == nogreen

ovr = Images.Overlay((gray, 0*gray), [RGB(1,0,1), RGB(0,1,0)], ([0,1],[0,1]))
@test_throws ErrorException Images.Overlay((gray, 0*gray), (RGB(1,0,1), RGB(0,1,0)), ((0,),(0,1)))
@test_throws ErrorException ovr[1] = RGB(0.2,0.4,0.6)

img1 = Images.Image(gray, Dict{ASCIIString,Any}([("colorspace","Gray"), ("spatialorder",["x"])]))
ovr = Images.OverlayImage((2gray,img1), (RGB{Ufixed8}(1,0,1), RGB{Ufixed8}(0,1,0)), ((0,1),(0,1)))
@test isa(ovr, Images.Image)
@test !haskey(ovr, "colorspace")
@test Images.colorspace(ovr) == "RGB"
@test ovr[2] == RGB{Float32}(0.5,0.25,0.5)
a = float32(rand(3,2))
b = float32(rand(3,2))
ovr = Images.OverlayImage((a,b), (RGB{Ufixed8}(1,0,1), RGB{Ufixed8}(0,1,0)), ((0,1),(0,1)))
@test isa(ovr, Images.Image)
@test abs(ovr[1,2] - RGB{Float32}(a[1,2],b[1,2],a[1,2])) < 1e-5


text:
import Observables
import Makie

scene = Makie.Scene()
width = 1.0
color = Makie.RGBAf(1.0, 0.0, 0.0, 1.0)
transparency = false
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrow = Arrow(tail,
              head,
              scene,
              width = width,
              color = color)
tail1 = Observables.to_value(arrow.tail)[1]
head1 = Observables.to_value(arrow.head)[1]
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
update(arrow, tail, head)
tail2 = Observables.to_value(arrow.tail)[1]
head2 = Observables.to_value(arrow.head)[1]

@test isapprox(tail1, tail2) == false
@test isapprox(head1, head2) == false

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
function runMode()
    players = getPlayerIds()

    while true
        sleep(0.05) #Don't overload the API & wait before removal
        for pl in players
            b1 = getTile(pl) .+ (0, -1, 0) #first block under the player
            b2 = b1 .+ (0, -1, 0) #second block under the player

            if getBlock(b1) == PiCraft.SAND && getBlock(b2) == PiCraft.TNT
                setBlock(b1, PiCraft.AIR)
                setBlock(b2, PiCraft.AIR)
            end
        end
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 35b500b2-9047-11eb-2570-338858be8b6a
html"""

Section 2.4

 Random Variables as Types

Lecture Video

"""

# ╔═╡ 5d62e16c-8fd9-11eb-1c44-0b0232614011
TableOfContents(aside=true)

# ╔═╡ 103cd2f4-903c-11eb-1116-a51dc540175c
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="Symbolics", version="0.1"),
    ])

    using Statistics, Plots, PlutoUI, LinearAlgebra, Symbolics
end

# ╔═╡ 6bbfa37e-8ffe-11eb-3031-19ea76a6a8d2
md"""
# Concepts for today

We won't be introducing any new Julia functions in this lecture, but we are going to show off a key way in which Julia really shines: its **type system**.

With a few carefully-chosen definitions, we can gain a great deal of power!

Concepts for today:

- Types for code organisation of abstract concepts
- Abstract types
- Subtypes
- Building up expressions "under the hood" using types

- `sum` over an iterator / generator expression
"""

# ╔═╡ 8a125ca0-8ff5-11eb-0607-45f993fb5ca7
md"""
# Random variables as types
"""

# ╔═╡ b2971770-8ff7-11eb-002c-f9dc9d6d0d70
md"""
This lecture might appear to be about random variables (and it is). But we would be thrilled if you see it rather as a more general discussion on software engineering and abstractions, since the principles involved extend to *many* different contexts.
"""

# ╔═╡ cdd4497c-903d-11eb-03be-abf6002e75e7
md"""
## Random variables
"""

# ╔═╡ ae1b3a26-8fd3-11eb-3746-ad48301ff96e
md"""
Recall that a **random variable** $X$ is an object which has different possible **outcomes** $x$, to which we assign **probabilities** $\mathbb{P}(X = x)$.

The correspondence of probabilities to outcomes is called the **probability distribution** of the random variable.
"""

# ╔═╡ ae0008ee-8fd3-11eb-38bd-f52598a97dce
md"""
## Gaussian distributions
"""

# ╔═╡ c6c3cf54-8fd4-11eb-3b4f-415f1a2da18e
md"""
Let's remind ourselves of a very important type of distribution, namely the **Gaussian distribution** with mean $\mu$ and standard deviation $\sigma$ (or variance $\sigma^2$).

We can sample from a Gaussian distribution with mean $0$ and variance $1$ with the `randn` function (short for "random normal"). We can then shift and scale to get the Gaussian distribution that we want:
"""

# ╔═╡ b11d964e-8fd4-11eb-3f6a-43e8d2fa462c
data = μ .+ σ .* randn(10^5)   # transform standard normal

# ╔═╡ d8b74772-8fd4-11eb-3943-f98c29d02171
md"""
μ = $(@bind μ Slider(-3:0.01:3, show_value=true, default=0.0))
σ = $(@bind σ Slider(0.01:0.01:3, show_value=true, default=1.0))
"""

# ╔═╡ b8c253f8-8fd4-11eb-304a-b1be962687e4
begin
	histogram(data, alpha=0.2, norm=true, bins=100, leg=false, title="μ=$(μ), σ=$(σ)", size=(500, 300))
	xlims!(-6, 6)
	ylims!(0, 0.7)

	xs = [μ - σ, μ, μ + σ]

	plot!(-6:0.01:6, x -> bell_curve(x, μ, σ), lw=2)

	plot!((μ - σ):0.01:(μ + σ), x -> bell_curve(x, μ, σ), fill=true, alpha=0.5, c=:purple)

	plot!([μ, μ], [0.05, bell_curve(μ, μ, σ)], ls=:dash, lw=2, c=:white)
	annotate!(μ, 0.03, text("μ", :white))
#	annotate!(μ + σ, 0.03, text("μ+σ", :yellow))
#	annotate!(μ, 0.03, text("μ", :white))

end

# ╔═╡ ad7b3bee-8fd5-11eb-06f6-b39738d4b1fd
bell_curve(x) = exp(-x^2 / 2) / √(2π)

# ╔═╡ c76cd760-8fd5-11eb-3500-5d15515c33f5
bell_curve(x, μ, σ) = bell_curve( (x - μ) / σ ) / σ

# ╔═╡ f31275fa-8fd5-11eb-0b76-7d0513705273
bell_curve(0, 3, 2)

# ╔═╡ 276a7c36-8fd8-11eb-25d8-3d4cfaa1f71c
md"""
### Sum of two Gaussians
"""

# ╔═╡ 11f3853c-903e-11eb-04cd-a125017ad5d8
md"""
Let's try to *add* two Gaussians. From the empirical (experimental) point of view the sum of two random variables is easy to compute: Just sample both of them and add the resulting values to get the value of the random variable that is their sum:
"""

# ╔═╡ 2be60570-8fd8-11eb-0bdf-951280dc6181
begin
	data1 = 4 .+ sqrt(0.3) .* randn(10^5)
	data2 = 6 .+ sqrt(0.7) .* randn(10^5)

	total = data1 + data2
end

# ╔═╡ 5e094b52-8fd8-11eb-2ac4-5797599ab013
histogram(total, alpha=0.5, leg=false, norm=true, size=(500, 300))

# ╔═╡ bd083794-8fd8-11eb-0155-e59fe27d64f2
Statistics.mean(total)

# ╔═╡ a32079ea-8fd8-11eb-01c4-81b603033b55
Statistics.std(total)

# ╔═╡ 79fb368c-8fd9-11eb-1c9c-bd0ceb122b11
md"""
We see that we get back a Gaussian again!
In fact, the sum of two Gaussians with means $\mu_1$ and $\mu_2$, and variances $\sigma_1^2$ and $\sigma_2^2$, is a Gaussian again! -- with mean $\mu_1 + \mu_2$ and variance $\sigma_1^2 + \sigma_2^2$. Below we will see how to encode this fact computationally.

"""

# ╔═╡ a2481afa-8fd3-11eb-1769-bf97f42ea79e
md"""
## Theoretical random variables vs. sampling

What should a theoretical random variable be able to do? What should we be able to do with sampling?

- **Naming**, e.g. `Gaussian`
- **Parameters**, e.g. $\mu$ and $\sigma^2$

### Theoretical

- Theoretical mean
- Theoretical variance

- Theoretical sum of  two random variables
- Theoretical product of two random variables

- Probability distribution

### Sampling

- **Sample** a  random variable
- Sample mean
- Sample variance

- Sample sum

- Histogram

"""

# ╔═╡ a9654334-8fd9-11eb-2ea8-8d308ea66963
md"""

# Why define a type at all?
"""

# ╔═╡ bb2132e0-8fd9-11eb-3bdd-594726c04859
md"""
How can we *represent a random variable* in software?

In some languages there are different names for the various functions associated to a random variable, but no name for the random variable itself!

For example, in R there is a standard naming convention in R, with `d` for density, etc., followed by the name like `norm` for normal. The indicators are:

- `d` for the density
- `p` for the distribution function
- `q` for the quantile function
- `r` for generating random variates

See e.g.
- [Normal distribution in R](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/Normal)

- [Chi-squared distribution in R](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/Chisquare)

What's wrong with this? All these functions are referring to an underlying random variable (or probability distribution), which you will find in any course in probability, and yet there's no way to refer to the underlying mathematical object!
"""

# ╔═╡ f307e3b8-8ff0-11eb-137e-4f9a03bb4d78
md"""
Instead, we would like to be able to refer to the random variable (or probability distribution) itself. We should be able to provide a type with the name and parameters of a random variable, but not yet specify how to generate random instances (variates)? It turns out that this is a very good example of thinking ahead by providing an **abstraction**.

We can *later* provide a means for random sampling -- and even, if a new algorithm comes along that is more efficient, we can replace it some day!
"""

# ╔═╡ e0ef47a6-903c-11eb-18aa-6ff06f0e28ac
md"""
## Defining abstract types for random variables
"""

# ╔═╡ 02051416-903d-11eb-0ade-3b20897989c5
md"""
We define an **abstract type** using `abstract type <Name> end`, where `<Name>` is replaced with the name of the type.

We can think of an abstract type as being a collection of types that share a particular property. In our case, we want to create a type to represent "any random variable", and also the sub-types "any continuous(-valued) random variable" and "any discrete(-valued) random variable".

This will allow us to specify later on whether a given concrete (i.e. particular) random variable is discrete or continuous.

We use `<:` to denote **sub-type**:
"""

# ╔═╡ 51ee3c3c-903d-11eb-1bfa-3bbcda98e977
begin
	abstract type RandomVariable end

	abstract type DiscreteRandomVariable <: RandomVariable end
	abstract type ContinuousRandomVariable <: RandomVariable end
end

# ╔═╡ 681429d8-8fff-11eb-0fa1-bbf9e05e6cea
md"""
## Defining a type for a Gaussian random variable
"""

# ╔═╡ 21236a86-8fda-11eb-1fcf-59d1de75470c
md"""
Let's start off by looking at **Gaussian** random variables, also called normal random variables. We reminded ourselves above what they look like.
"""

# ╔═╡ dd130aae-8ff2-11eb-2b15-2f5123b40d20
md"""
### Name and parameters
"""

# ╔═╡ 4771e8e6-8fd2-11eb-178c-419cbdb348f4
begin
	struct Gaussian <: ContinuousRandomVariable
		μ     # mean
		σ²    # variance
	end

	Gaussian() = Gaussian(0.0, 1.0)  # normalised Gaussian with mean 0 and variance 1
end

# ╔═╡ ece14f6c-8fda-11eb-279b-c18dc0fc46d7
G = Gaussian(1, 2)

# ╔═╡ 0be6c548-8fd3-11eb-07a2-5bb382614cab
md"""
Note that here we have created a Gaussian random variable with given parameter values, *without* sampling from it.
"""

# ╔═╡ 9c03814c-8ff2-11eb-189d-f3b0507507cb
md"""
### Theoretical mean and variance
"""

# ╔═╡ f24bcb1a-8fda-11eb-1828-0307c65b81cd
md"""
Now we can extend the `mean`, `var` (variance) and `std` (standard deviation) functions from the `Statistics` library to act on this object:
"""

# ╔═╡ 0c76f49e-8fdb-11eb-2c2c-59a4060e8e1d
begin
	Statistics.mean(X::Gaussian) = X.μ
	Statistics.var(X::Gaussian) = X.σ²
end

# ╔═╡ 4596edec-8fdb-11eb-1369-f3e98cb831cd
mean(G)

# ╔═╡ 255d51b6-8ff3-11eb-076c-31ba4c5ce10d
var(G)

# ╔═╡ 43292298-8ff7-11eb-2b9d-017acc9ef185
md"""
### Planning ahead: Standard deviation for *any* random variable, not just Gaussians
"""

# ╔═╡ 1bdc6068-8ff7-11eb-069c-6d0f1a83f373
md"""
Once we have defined the variance, we know how to calculate the standard deviation: it's just the square root of the variance. But, thinking ahead, this is true for *any* random variable, so we can define it to act on any random variable that we will define later!:
"""

# ╔═╡ 11ef1d18-8ff7-11eb-1645-113c1aae6e9b
Statistics.std(X::RandomVariable) = sqrt(var(X))

# ╔═╡ 592afd96-8ff7-11eb-1200-47eaf6b6a28a
md"""
This is an example of good software design.
"""

# ╔═╡ 4ef2153a-8fdb-11eb-1a23-8b1c28381e34
std(G)

# ╔═╡ 76f4424e-8fdc-11eb-35ee-dd09752b947b
md"""
### Sum of two Gaussian random variables
"""

# ╔═╡ dede8022-8fd2-11eb-22f8-a5614d703c01
md"""
Gaussians have a special property: the sum of two Gaussians is always a Gaussian. (We say that Gaussians are **stable** distributions; [there are others](https://en.wikipedia.org/wiki/Stable_distribution).) Note that we don't need random samples for this theoretical observation, embodied in the following code:
"""

# ╔═╡ b4f18188-8fd2-11eb-1950-ef3bee3e4724
Base.:+(X::Gaussian, Y::Gaussian) = Gaussian(X.μ + Y.μ, X.σ² + Y.σ²)

# ╔═╡ 292f0e44-8fd3-11eb-07cf-5733a5327630
begin
	G1 = Gaussian(0, 1)
	G2 = Gaussian(5, 6)
end

# ╔═╡ 6b50e9e6-8fd3-11eb-0ab3-6fc3efea7e37
G1 + G2

# ╔═╡ 5300b082-8fdb-11eb-03fe-55f511364ad9
mean(G1 + G2) == mean(G1) + mean(G2)

# ╔═╡ ac3a7c18-8ff3-11eb-1c1f-17939cff40f9
md"""
The theoretical product of two Gaussians is not Gaussian; we will do the general case later.
"""

# ╔═╡ b8e59db2-8ff3-11eb-1fa1-595fdc6652bc
md"""
### Probability distribution of a Gaussian
"""

# ╔═╡ cb167a60-8ff3-11eb-3885-d1d8a70e2a5e
md"""
A Gaussian random variable is a **continuous** random variable, i.e. it has a continuous range of possible outcomes. The possible range of outcomes is called the **support** of the distribution. For a Gaussian it is the whole real line, $(-\infty, \infty)$.
"""

# ╔═╡ fc0acee6-8ff3-11eb-13d1-1350364f03a9
md"""
One way to specify a continous random variable $X$ is via its **probability density function**, or **PDF**, $f_X$. The probability that $X$ lies in the interval $[a, b]$ is given by an area under the curve $f_X(x)$ from $a$ to $b$:

$$\mathbb{P}(X \in [a, b]) = \int_{a}^b f_X(x) \, dx.$$
"""

# ╔═╡ 2c8b75de-8ff4-11eb-1bc6-cde7b67a2007
md"""
For a Gaussian distribution with mean $\mu$ and variance $\sigma^2$, the PDF is given by

$$f_X(X) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left[ -\frac{1}{2} \left( \frac{x - \mu}{\sigma} \right)^2 \right]$$.
"""

# ╔═╡ 639e517c-8ff4-11eb-3ea8-07b73b0bff78
pdf(X::Gaussian) = x -> exp(-0.5 * ( (x - X.μ)^2 / X.σ²) ) / √(2π * X.σ²)

# ╔═╡ 3d78116c-8ff5-11eb-0ee3-71bf9413f30e
pdf(G)

# ╔═╡ b63c02e4-8ff4-11eb-11f6-d760e8760118
pdf(Gaussian())(0.0)

# ╔═╡ a899fc08-8ff5-11eb-1d00-e95b55c3be4b
begin
	plot(pdf(Gaussian(μμ, σσ)), leg=false)
	xlims!(-6, 6)
	ylims!(0, 0.5)
end

# ╔═╡ cd9b10f0-8ff5-11eb-1671-99c6738b8074
md"""
μ = $(@bind μμ Slider(-3:0.01:3, show_value=true, default=0.0))
σ = $(@bind σσ Slider(0.01:0.01:3, show_value=true, default=1.0))
"""

# ╔═╡ 0b901c34-8ff6-11eb-225b-511718412309
md"""
### Sampling from a Gaussian distribution
"""

# ╔═╡ 11cbf48a-903f-11eb-1e77-81eeb358ec24
md"""
We can also specify how to sample from a Gaussian distribution. We can re-purpose `rand` for this!
"""

# ╔═╡ 180a4746-8ff6-11eb-046f-ddf6bb938a35
Base.rand(X::Gaussian) = X.μ + √(X.σ²) * randn()

# ╔═╡ 27595716-8ff6-11eb-1d99-7bf13eedddf7
histogram!([rand(Gaussian(μμ, σσ)) for i in 1:10^4], alpha=0.5, norm=true)

# ╔═╡ bfb62f1a-8ff5-11eb-0f7a-cf725f3269c5
md"""
# More general distributions
"""

# ╔═╡ 0a150880-8ff6-11eb-0046-45fa2d4b476e
md"""
Let's recall the Bernoulli distribution from last lecture. This represents a weighted coin with probability $p$ to come up "heads" (1), and probability $1-p$ to come up "tails" (0).

Note that this is a **discrete** random variable: the possible outcomes are the discrete values $0$ and $1$.
"""

# ╔═╡ 26569eb4-8fd1-11eb-2df9-098a0792a09e
struct Bernoulli <: DiscreteRandomVariable
	p::Float64
end

# ╔═╡ baf1fe40-8ff6-11eb-1da1-cd43880db334
B = Bernoulli(0.25)

# ╔═╡ 38b10d94-903f-11eb-3e26-890382342dc1
md"""
Again we can specify the theoretical mean and variance:
"""

# ╔═╡ be5f9900-8ff6-11eb-1816-03447cabd9a9
begin
	Statistics.mean(X::Bernoulli) = X.p
	Statistics.var(X::Bernoulli) = X.p * (1 - X.p)
end

# ╔═╡ 52a2cac8-8ff8-11eb-2da4-0b113618c64b
mean(B), var(B), std(B)

# ╔═╡ 6a06b64e-903f-11eb-0d23-53f223ed1ed3
md"""
And again `std` just works!
"""

# ╔═╡ 6e9ad41a-903f-11eb-06ae-1b34a6674eaf
md"""
Finally we specify how to sample:
"""

# ╔═╡ 754b0a80-8ff8-11eb-364b-85fa49d6bb8e
Base.rand(X::Bernoulli) = Int(rand() < X.p)

# ╔═╡ bd484878-8ff8-11eb-0337-7730ab2b07d4
md"""
## Adding two random variables
"""

# ╔═╡ c3043eb6-8ff8-11eb-2cae-13d1c1613234
md"""
What happens if we add two Bernoulli random variables? There are two routes we could go: We could use the known theoretical sum, or we could write a general-purpose tool. Let's do the latter.
"""

# ╔═╡ 7e3cea0c-903f-11eb-0b41-0f381c1cce4b
md"""
When we add two Bernoulli random variables we do *not* get a Bernoulli back. To see this it's enough to observe that the sum can have the outcome 2, which is impossible for a Bernoulli.

So the result is just the random variable "the sum of these two given random variables". In general it won't even have a common name.

So we actually need to *define a new type* to represent the "sum of two given random variables", which itself will be a random variable!:

"""

# ╔═╡ eb555508-8ff8-11eb-1b70-e95290084742
struct SumOfTwoRandomVariables <: RandomVariable
	X1::RandomVariable
	X2::RandomVariable
end

# ╔═╡ 1e5047c4-8ff9-11eb-2628-c7824725678a
begin
	B1 = Bernoulli(0.25)
	B2 = Bernoulli(0.6)
end

# ╔═╡ cb2cd908-903f-11eb-09c0-99acde6d765c
md"""
Now we can define the sum of two random variables of *any* type:
"""

# ╔═╡ 44a5ef96-8ff9-11eb-06a0-d3a8dcf5c1aa
Base.:+(X1::RandomVariable, X2::RandomVariable) = SumOfTwoRandomVariables(X1, X2)

# ╔═╡ d8437e6c-903f-11eb-3ac2-5f7c380c0872
md"""
For example, let's sum two Bernoullis:
"""

# ╔═╡ 574744ec-8ff9-11eb-033c-a3dff07a292b
B1 + B2

# ╔═╡ e024377a-903f-11eb-316a-b5b7936e610f
md"""
However, for the special case of Gaussians we still get the correct result (we have *not* overwritten the previous definition):
"""

# ╔═╡ f37bb49c-903f-11eb-03fb-35d6ac35822d
G1 + G2

# ╔═╡ 318f5274-8ff9-11eb-1c88-5fde5b546099
md"""
Now we need to define the various functions on this type representing a sum
"""

# ╔═╡ 90ee7558-8ff9-11eb-3602-271890987ece
Statistics.mean(S::SumOfTwoRandomVariables) = mean(S.X1) + mean(S.X2)

# ╔═╡ a168ace6-8ff9-11eb-393d-45c34fdf577c
mean(B1 + B2)

# ╔═╡ a51b1538-8ff9-11eb-0635-81088e826bb3
md"""
To have a simple equation for the variance, we need to assume that the two random variables are **independent**. Perhaps the name should have been `SumOfTwoIndependentRandomVariables`, but it seems too long.
"""

# ╔═╡ d88c0830-8ff9-11eb-3e71-c1ac327f4e25
Statistics.var(S::SumOfTwoRandomVariables) = var(S.X1) + var(S.X2)

# ╔═╡ dd995120-8ff9-11eb-1a53-2d65f0b5585a
md"""
How can we sample from the sum? It's actually easy!
"""

# ╔═╡ f171c27c-8ff9-11eb-326c-6b2d8c38451d
Base.rand(S::SumOfTwoRandomVariables) = rand(S.X1) + rand(S.X2)

# ╔═╡ ab7c6b18-8ffa-11eb-2b6c-dde3dca1c6f7
md"""
Now it's easy to look at the sum of a Bernoulli and a Gaussian. This is an example of a [**mixture distribution**](https://en.wikipedia.org/wiki/Mixture_distribution).
"""

# ╔═╡ 79c55fc0-8ffb-11eb-229f-49198aee8245
md"""
Let's extend the `histogram` function to easily draw the histogram of a random variable:
"""

# ╔═╡ 4493ac9e-8ffb-11eb-15c9-a146091d7696
Plots.histogram(X::RandomVariable; kw...) = histogram([rand(X) for i in 1:10^6], norm=true, leg=false, alpha=0.5, size=(500, 300), kw...)

# ╔═╡ a591ac12-8ffb-11eb-1734-09657f875b1e
histogram(Bernoulli(0.25) + Bernoulli(0.75))

# ╔═╡ 6c5ef6d4-8ffb-11eb-2fe9-9dd87d92abe4
histogram(Bernoulli(0.25) + Gaussian(0, 0.1))

# ╔═╡ 1d4e236c-8ffb-11eb-209d-851e1af231d4
md"""
Now... What if we sum more random variables?
"""

# ╔═╡ 0c3cfb16-8ffb-11eb-3ef9-33ea9acbb8c0
mixture = Bernoulli(0.25) + Bernoulli(0.75) + Gaussian(0, 0.1)

# ╔═╡ 325560f4-8ffb-11eb-0e9b-53f6869bdd97
rand( mixture )

# ╔═╡ 89138e02-8ffb-11eb-2ad2-c32e663f57b0
histogram( mixture )

# ╔═╡ 71cf1724-9040-11eb-25c3-69ccde4abf0d
md"""
## Generic programming: `sum`
"""

# ╔═╡ 646ac706-9040-11eb-297c-7b6d8fdf3c16
md"""
Now we have defined `+`, Julia's generic definition of `sum` can kick in to define the sum of many random variables!:
"""

# ╔═╡ 34bcab72-8ffb-11eb-1d0c-29bd83de638b
S = sum(Bernoulli(0.25) for i in 1:30)

# ╔═╡ 77b54a34-9040-11eb-00d1-57f712ff6732
md"""
Note that we do not need the `[...]` in the following expression. There is no need to actually create an array of random variables; instead we are using an **iterator** or **generator expression**:
"""

# ╔═╡ bf9a4722-8ffb-11eb-1652-f1bfb4916d2a
histogram(S)

# ╔═╡ 611f770c-8ffc-11eb-2c23-0d5750afd7c8
mean(S)

# ╔═╡ 636b2ce0-8ffc-11eb-2411-571f78fb8a84
var(S)

# ╔═╡ c7613c44-8ffc-11eb-0a43-dbc8b1f62be9
rand(S)

# ╔═╡ 656999be-8ffc-11eb-0f51-51d5a162a004
md"""
This is a big deal! Everything just works.

By the way, the sum of $n$ Bernoulli random variables with the *same* probability $p$ is called a **binomial** random variable with parameters $(n, p)$.
"""

# ╔═╡ 8cb5373a-8ffc-11eb-175d-b11ec5fae1ab
md"""
If we were worried about performance, we would probably want to define a separate `Binomial` type, rather than using nested sums of Bernoullis. You will do this in the homework.
"""

# ╔═╡ 25b4ccca-8ffd-11eb-3fe0-2d14ce90d8b3
md"""
## χ₁² distribution
"""

# ╔═╡ 9913427a-8ffd-11eb-3c28-c12ee05ed9d3
md"""
Another well-known distribution is the chi-squared distribution, which is defined as the sum of squares of Gaussians. Chi-1 squared ($\chi_1^2$) is the square of a single Gaussian:
"""

# ╔═╡ aa9d94dc-8ffd-11eb-0537-c39344c224c2
struct ChiSquared1 <: ContinuousRandomVariable
end

# ╔═╡ b5300720-8ffd-11eb-2ca4-ed7a7c9b5516
Base.rand(X::ChiSquared1) = rand(Gaussian())^2

# ╔═╡ eea0c788-8ffd-11eb-06e0-cfb65dbba7f1
histogram(ChiSquared1())

# ╔═╡ c57fb088-9040-11eb-3881-dfbd82e72b35
md"""
We can get the $\chi_n^2$ distribution by summing several $\chi_1^2$:
"""

# ╔═╡ 028695a2-8ffe-11eb-3cf4-bd07bfe4df3a
histogram(sum(ChiSquared1() for i in 1:4))

# ╔═╡ 820af6dc-8fdc-11eb-1792-ad9f32eb915e
md"""
## Using symbolics
"""

# ╔═╡ ce330b78-8fda-11eb-144d-756275542eea
md"""
We can even do some of this symbolically!, e.g.:
"""

# ╔═╡ 668ade6a-8fda-11eb-0a5e-cb791b245ec0
@variables μ₁, σ₁², μ₂, σ₂²   # introduce symbolic variables from Symbolics.jl

# ╔═╡ 4e94d1da-8fda-11eb-3572-398fb4a12c3c
Gaussian(μ₁, σ₁²) + Gaussian(μ₂, σ₂²)

# ╔═╡ 9ca993a6-8fda-11eb-089c-4d7f89c81b94
Gaussian(17, 3) + Gaussian(μ₂, σ₂²)

# ╔═╡ Cell order:
# ╟─35b500b2-9047-11eb-2570-338858be8b6a
# ╠═5d62e16c-8fd9-11eb-1c44-0b0232614011
# ╠═103cd2f4-903c-11eb-1116-a51dc540175c
# ╟─6bbfa37e-8ffe-11eb-3031-19ea76a6a8d2
# ╟─8a125ca0-8ff5-11eb-0607-45f993fb5ca7
# ╟─b2971770-8ff7-11eb-002c-f9dc9d6d0d70
# ╟─cdd4497c-903d-11eb-03be-abf6002e75e7
# ╟─ae1b3a26-8fd3-11eb-3746-ad48301ff96e
# ╟─ae0008ee-8fd3-11eb-38bd-f52598a97dce
# ╟─c6c3cf54-8fd4-11eb-3b4f-415f1a2da18e
# ╠═b11d964e-8fd4-11eb-3f6a-43e8d2fa462c
# ╟─d8b74772-8fd4-11eb-3943-f98c29d02171
# ╟─b8c253f8-8fd4-11eb-304a-b1be962687e4
# ╠═ad7b3bee-8fd5-11eb-06f6-b39738d4b1fd
# ╠═c76cd760-8fd5-11eb-3500-5d15515c33f5
# ╠═f31275fa-8fd5-11eb-0b76-7d0513705273
# ╟─276a7c36-8fd8-11eb-25d8-3d4cfaa1f71c
# ╟─11f3853c-903e-11eb-04cd-a125017ad5d8
# ╠═2be60570-8fd8-11eb-0bdf-951280dc6181
# ╠═5e094b52-8fd8-11eb-2ac4-5797599ab013
# ╠═bd083794-8fd8-11eb-0155-e59fe27d64f2
# ╠═a32079ea-8fd8-11eb-01c4-81b603033b55
# ╟─79fb368c-8fd9-11eb-1c9c-bd0ceb122b11
# ╟─a2481afa-8fd3-11eb-1769-bf97f42ea79e
# ╟─a9654334-8fd9-11eb-2ea8-8d308ea66963
# ╟─bb2132e0-8fd9-11eb-3bdd-594726c04859
# ╟─f307e3b8-8ff0-11eb-137e-4f9a03bb4d78
# ╟─e0ef47a6-903c-11eb-18aa-6ff06f0e28ac
# ╟─02051416-903d-11eb-0ade-3b20897989c5
# ╠═51ee3c3c-903d-11eb-1bfa-3bbcda98e977
# ╟─681429d8-8fff-11eb-0fa1-bbf9e05e6cea
# ╟─21236a86-8fda-11eb-1fcf-59d1de75470c
# ╟─dd130aae-8ff2-11eb-2b15-2f5123b40d20
# ╠═4771e8e6-8fd2-11eb-178c-419cbdb348f4
# ╟─ece14f6c-8fda-11eb-279b-c18dc0fc46d7
# ╟─0be6c548-8fd3-11eb-07a2-5bb382614cab
# ╟─9c03814c-8ff2-11eb-189d-f3b0507507cb
# ╟─f24bcb1a-8fda-11eb-1828-0307c65b81cd
# ╠═0c76f49e-8fdb-11eb-2c2c-59a4060e8e1d
# ╠═4596edec-8fdb-11eb-1369-f3e98cb831cd
# ╠═255d51b6-8ff3-11eb-076c-31ba4c5ce10d
# ╟─43292298-8ff7-11eb-2b9d-017acc9ef185
# ╟─1bdc6068-8ff7-11eb-069c-6d0f1a83f373
# ╠═11ef1d18-8ff7-11eb-1645-113c1aae6e9b
# ╟─592afd96-8ff7-11eb-1200-47eaf6b6a28a
# ╠═4ef2153a-8fdb-11eb-1a23-8b1c28381e34
# ╟─76f4424e-8fdc-11eb-35ee-dd09752b947b
# ╟─dede8022-8fd2-11eb-22f8-a5614d703c01
# ╠═b4f18188-8fd2-11eb-1950-ef3bee3e4724
# ╠═292f0e44-8fd3-11eb-07cf-5733a5327630
# ╠═6b50e9e6-8fd3-11eb-0ab3-6fc3efea7e37
# ╠═5300b082-8fdb-11eb-03fe-55f511364ad9
# ╟─ac3a7c18-8ff3-11eb-1c1f-17939cff40f9
# ╟─b8e59db2-8ff3-11eb-1fa1-595fdc6652bc
# ╟─cb167a60-8ff3-11eb-3885-d1d8a70e2a5e
# ╟─fc0acee6-8ff3-11eb-13d1-1350364f03a9
# ╟─2c8b75de-8ff4-11eb-1bc6-cde7b67a2007
# ╠═639e517c-8ff4-11eb-3ea8-07b73b0bff78
# ╠═3d78116c-8ff5-11eb-0ee3-71bf9413f30e
# ╠═b63c02e4-8ff4-11eb-11f6-d760e8760118
# ╠═a899fc08-8ff5-11eb-1d00-e95b55c3be4b
# ╟─cd9b10f0-8ff5-11eb-1671-99c6738b8074
# ╟─0b901c34-8ff6-11eb-225b-511718412309
# ╟─11cbf48a-903f-11eb-1e77-81eeb358ec24
# ╠═180a4746-8ff6-11eb-046f-ddf6bb938a35
# ╠═27595716-8ff6-11eb-1d99-7bf13eedddf7
# ╟─bfb62f1a-8ff5-11eb-0f7a-cf725f3269c5
# ╟─0a150880-8ff6-11eb-0046-45fa2d4b476e
# ╠═26569eb4-8fd1-11eb-2df9-098a0792a09e
# ╠═baf1fe40-8ff6-11eb-1da1-cd43880db334
# ╟─38b10d94-903f-11eb-3e26-890382342dc1
# ╠═be5f9900-8ff6-11eb-1816-03447cabd9a9
# ╟─52a2cac8-8ff8-11eb-2da4-0b113618c64b
# ╟─6a06b64e-903f-11eb-0d23-53f223ed1ed3
# ╟─6e9ad41a-903f-11eb-06ae-1b34a6674eaf
# ╠═754b0a80-8ff8-11eb-364b-85fa49d6bb8e
# ╟─bd484878-8ff8-11eb-0337-7730ab2b07d4
# ╟─c3043eb6-8ff8-11eb-2cae-13d1c1613234
# ╟─7e3cea0c-903f-11eb-0b41-0f381c1cce4b
# ╠═eb555508-8ff8-11eb-1b70-e95290084742
# ╠═1e5047c4-8ff9-11eb-2628-c7824725678a
# ╟─cb2cd908-903f-11eb-09c0-99acde6d765c
# ╠═44a5ef96-8ff9-11eb-06a0-d3a8dcf5c1aa
# ╟─d8437e6c-903f-11eb-3ac2-5f7c380c0872
# ╠═574744ec-8ff9-11eb-033c-a3dff07a292b
# ╟─e024377a-903f-11eb-316a-b5b7936e610f
# ╠═f37bb49c-903f-11eb-03fb-35d6ac35822d
# ╟─318f5274-8ff9-11eb-1c88-5fde5b546099
# ╠═90ee7558-8ff9-11eb-3602-271890987ece
# ╠═a168ace6-8ff9-11eb-393d-45c34fdf577c
# ╟─a51b1538-8ff9-11eb-0635-81088e826bb3
# ╠═d88c0830-8ff9-11eb-3e71-c1ac327f4e25
# ╟─dd995120-8ff9-11eb-1a53-2d65f0b5585a
# ╠═f171c27c-8ff9-11eb-326c-6b2d8c38451d
# ╟─ab7c6b18-8ffa-11eb-2b6c-dde3dca1c6f7
# ╟─79c55fc0-8ffb-11eb-229f-49198aee8245
# ╠═4493ac9e-8ffb-11eb-15c9-a146091d7696
# ╠═a591ac12-8ffb-11eb-1734-09657f875b1e
# ╠═6c5ef6d4-8ffb-11eb-2fe9-9dd87d92abe4
# ╟─1d4e236c-8ffb-11eb-209d-851e1af231d4
# ╠═0c3cfb16-8ffb-11eb-3ef9-33ea9acbb8c0
# ╠═325560f4-8ffb-11eb-0e9b-53f6869bdd97
# ╠═89138e02-8ffb-11eb-2ad2-c32e663f57b0
# ╟─71cf1724-9040-11eb-25c3-69ccde4abf0d
# ╟─646ac706-9040-11eb-297c-7b6d8fdf3c16
# ╠═34bcab72-8ffb-11eb-1d0c-29bd83de638b
# ╟─77b54a34-9040-11eb-00d1-57f712ff6732
# ╠═bf9a4722-8ffb-11eb-1652-f1bfb4916d2a
# ╠═611f770c-8ffc-11eb-2c23-0d5750afd7c8
# ╠═636b2ce0-8ffc-11eb-2411-571f78fb8a84
# ╠═c7613c44-8ffc-11eb-0a43-dbc8b1f62be9
# ╟─656999be-8ffc-11eb-0f51-51d5a162a004
# ╟─8cb5373a-8ffc-11eb-175d-b11ec5fae1ab
# ╟─25b4ccca-8ffd-11eb-3fe0-2d14ce90d8b3
# ╟─9913427a-8ffd-11eb-3c28-c12ee05ed9d3
# ╠═aa9d94dc-8ffd-11eb-0537-c39344c224c2
# ╠═b5300720-8ffd-11eb-2ca4-ed7a7c9b5516
# ╠═eea0c788-8ffd-11eb-06e0-cfb65dbba7f1
# ╟─c57fb088-9040-11eb-3881-dfbd82e72b35
# ╠═028695a2-8ffe-11eb-3cf4-bd07bfe4df3a
# ╟─820af6dc-8fdc-11eb-1792-ad9f32eb915e
# ╟─ce330b78-8fda-11eb-144d-756275542eea
# ╠═668ade6a-8fda-11eb-0a5e-cb791b245ec0
# ╠═4e94d1da-8fda-11eb-3572-398fb4a12c3c
# ╠═9ca993a6-8fda-11eb-089c-4d7f89c81b94


text:
@testset "100.same-tree.jl" begin
    @test is_same_tree(TreeNode{Int}([1, 2, 3]), TreeNode{Int}([1, 2, 3])) == true
    @test is_same_tree(TreeNode{Int}([1, 2]), TreeNode([1, nothing, 2])) == false
    @test is_same_tree(TreeNode{Int}([1, 2, 1]), TreeNode([1, 1, 2])) == false
end


text:
using StitchDPI, Test

#nanplus
a = ones(10,10)
b = deepcopy(a)
a[:,1] .= NaN
c = StitchDPI.nanplus(a,b)
@test all(c[:,1] .== 1.0)
@test all(c[:,2:end] .== 2.0)

#crop_vals
cr = StitchDPI.crop_vals(a, NaN)
@test size(cr,2) == size(a,2) - 1
@test size(cr,1) == size(a,1)
@test all(cr.==a[:,2:end])
a[8:end,:] .= NaN
cr = StitchDPI.crop_vals(a, NaN)
@test size(cr,2) == size(a,2) - 1
@test size(cr,1) == size(a,1) - 3
@test all(cr.==a[1:7,2:end])

#ypad
r = rand(10,10)
rp = StitchDPI.ypad(r, 12; flip_y = false, fillval=NaN)
@test all(isnan.(rp[:,1]))
@test all(isnan.(rp[:,12]))
@test all(rp[:,2:11].==r)

#flipy
rpf = StitchDPI.ypad(r, 12; flip_y = true, fillval=NaN)
@test all(reverse(rpf[1,2:11]).==rp[1,2:11])
@test all(reverse(rpf[10,2:11]).==rp[10,2:11])


text:
function ci()
  @info ENV["GITHUB_EVENT_PATH"]
  f = JSON.parsefile(ENV["GITHUB_EVENT_PATH"])
  event = GitHub.event_from_payload!("commit", f)

  repo = event.repository

  # Fails if push to repo directly
  # issue => pull_request
  if haskey(event.payload, "pull_request")
	  reply_to = event.payload["pull_request"]["number"]
	  comment_kind = :issue
	else
		reply_to = event.payload["commits"] |> last
	  comment_kind = :commit
	  reply_to = reply_to["id"]
	end

	if comment_kind == :issue
	  files = [f.filename for f in GitHub.pull_request_files(repo, reply_to, auth = myauth)]
	  possible_models = [split(f, "/") for f in files]
	  possible_models = [length(f) > 1 ? f : [] for f in possible_models]
	  possible_models = [f[1:end-1] for f in possible_models]
	  possible_models = union(possible_models...)
	  possible_models = uppercase.(possible_models)

	  # meta = TOML.parsefile(joinpath(pwd(), "scripts", "Notebooks.toml"))
	  meta = ["MNIST", "CPPN", ]
	  available_models = keys(meta)

	  model = intersect(possible_models, available_models)
	  model = join(model, ' ')
	  @show model
	else
		model = "MNIST"
	end

	# if
  # Handle when model is not found
  if all(map(isspace, collect(model)))
  	# comment_kind == :commit && return
  	mas = "No matching models found, consider adding the appropriate models to the `Notebooks.toml`."
		# throw(ErrorException(msg))
		@info msg
  	return HTTP.Response(200)
  end

  # resp = trigger_pipeline(reply_to, model, event, fluxbot = false)
  resp = Dict("id" => "98138293",
              "web_url" => "https://gitlab.com/JuliaGPU/Flux.jl/pipelines/98138293",
              "sha" => "fbb377a7b436327c298c536ecb9d2ff5ee8e07d4")

  # GitHub.create_comment(event.repository, reply_to, comment_kind,
  #                       auth = myauth,
  #                       params = placeholder_resp(reply_to, resp, comment_kind = comment_kind))
end


text:
module FunctionWranglers

export FunctionWrangler, smap!, sfindfirst, sreduce, sindex, railway, Fail

import Base: length, show

struct FunctionWrangler{TOp, TNext}
    op::TOp
    next::TNext
end

function FunctionWrangler(fns)
    if length(fns) == 0
        return FunctionWrangler(nothing, nothing)
    end
    next = FunctionWrangler(fns[2:end])
    return FunctionWrangler(fns[1], next)
end

Base.length(w::FunctionWrangler) = isnothing(w.op) ? 0 : 1 + length(w.next)

Base.show(io::IO, ::MIME"text/plain", w::FunctionWrangler) = begin
    print(io, "FunctionWrangler with $(length(w)) items: ")
    iw = w
    while !isnothing(iw.op)
        print(io, string(iw.op) * ", ")
        iw = iw.next
    end
    return nothing
end

# Using variable argument definitions sometimes caused the compiler to "loose track" of type information
# and generate slow code (tested up to 1.5.1 with https://gist.github.com/tisztamo/1ce6632d7e6ffc4
# @generate-ing solves the issue
@inline @generated function _smap!(outputs, wrangler::FunctionWrangler{TOp, TNext}, myidx, args...) where {TNext, TOp}
    TOp === Nothing && return :(nothing)
    argnames = [:(args[$i]) for i = 1:length(args)]
    return quote
        outputs[myidx] = wrangler.op($(argnames...))
        return _smap!(outputs, wrangler.next, myidx + 1, $(argnames...))
    end
end

"""
    smap!(outputs, wrangler::FunctionWrangler, args...)

Map a single set of arguments using all the functions into a preallocated array.
"""
smap!(outputs, wrangler::FunctionWrangler, args...) = _smap!(outputs, wrangler, 1, args...)

@inline @generated function _sfindfirst(wrangler::FunctionWrangler{TOp, TNext}, myidx, args...) where {TNext, TOp}
    TOp === Nothing && return :(nothing)
    argnames = [:(args[$i]) for i = 1:length(args)]
    return quote
        if wrangler.op($(argnames...)) === true
            return myidx
        end
        return _sfindfirst(wrangler.next, myidx + 1, $(argnames...))
    end
end

"""
    sfindfirst(wrangler::FunctionWrangler, args...)

Look for the first function which returns `true` for the given arguments, and returns its index.
Return `nothing` if no function returned `true`.
"""
sfindfirst(wrangler::FunctionWrangler, args...) = _sfindfirst(wrangler, 1, args...)

@inline @generated function _sindex(wrangler::FunctionWrangler{TOp, TNext}, idx, args...) where {TNext, TOp}
    argnames = [:(args[$i]) for i = 1:length(args)]
    return quote
        if idx == 1
            return wrangler.op($(argnames...))
        end
        return _sindex(wrangler.next, idx - 1, $(argnames...))
    end
end

"""
    sindex(wrangler::FunctionWrangler, idx, args...)

Call the `idx`-th function with args.

Note that this call iterates the wrangler from 1 to `idx`. Try to
put the frequently called functions at the beginning to minimize overhead.
"""
function sindex(wrangler::FunctionWrangler, idx, args...)
    _sindex(wrangler, idx, args...)
end

@inline @generated function _sreduce(wrangler::FunctionWrangler{TOp, TNext}, myidx, prev, args...; init = nothing) where {TNext, TOp}
    TOp === Nothing && return :(prev)
    argnames = [:(args[$i]) for i = 1:length(args)]
    return quote
        current = wrangler.op(prev, $(argnames...))
        return _sreduce(wrangler.next, myidx + 1, current, $(argnames...))
    end
end

sreduce(wrangler::FunctionWrangler, args...; init = nothing) = _sreduce(wrangler, 1, init, args...)

abstract type Fail end

@inline function _railway(wrangler::FunctionWrangler, myidx, prev)
    isnothing(wrangler.op) && return prev
    current = wrangler.op(prev)
    current isa Fail && return current
    return _railway(wrangler.next, myidx + 1, current)
end

function railway(wrangler::FunctionWrangler, input)
    return _railway(wrangler, 1, input)
end

end # module


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

JLLWrappers.@generate_main_file_header("NOVAS")
JLLWrappers.@generate_main_file("NOVAS", UUID("0c6c9120-0c4d-51e7-92ed-b5d3e78608aa"))
end  # module NOVAS_jll


text:
using Documenter, ACEatoms

makedocs(sitename="ACEatoms.jl Documentation",
         pages = [
        "Home" => "index.md",
        "Getting Started" => "gettingstarted.md",
        "Linear Site E Models" => "linear.md",
        "Developer Docs" => "devel.md",
        "Types & Functions" => "docs.md"
         ])

deploydocs(
    repo = "github.com/ACEsuit/ACEatoms.jl.git",
    devbranch = "main"
)


text:
_getfields(obj) = map(n -> getfield(obj, n), fieldnames(typeof(obj))) :: Tuple

prefixof(f) = join(fullname(parentmodule(f)), '.')

function print_apply(io, f, args)
    if !get(io, :limit, false)
        # Don't show full name in REPL etc.:
        print(io, prefixof(f), '.')
    end
    print(io, nameof(f))
    if length(args) == 1
        print(io, '(')
        show(io, args[1])
        print(io, ')')
    else
        show(io, args)
    end
    return
end

_default_show(io, obj) = print_apply(io, typeof(obj), _getfields(obj))


text:
# Find blobs in frame
# Find connected components in a binary mask
# Label components 1 ... number of compoents

# function doclustering(f)
#     doclustering(copy(f))
# end
"""
    Find connected components in the 2d binary image given as a matrix of integers,
    which should only be 0 or 1.
    Returns the input matrix with values replaced by componentids 1:maxid, and maxid.
"""
function doclustering!(f, clusteridstart=one(eltype(f)))
    f1 = f
    clusterid, dict = findclusters1!(f1, clusteridstart)
    d, key = mergeclusters(dict)
    clusterid = changeclusters!(f1, d, key)
    f1, clusterid
end

function findclusters1!(f1::AbstractArray{T}, clustermax::T) where T
    # Connect each node to its left or top neighbour
    # If left and top both exists then make a connection
    edges = Set{Tuple{T,T}}()
    clusterid = clustermax + one(T)
    i = 1
    for j in 2:size(f1, 2)
        cur = f1[i, j]
        if cur != 0
            top = f1[i, j-1]
            if top != 0
                cur = top
            else
                cur = clusterid
                clusterid += one(T)
            end
            f1[i, j] = cur
        end
    end
    j = 1
    for i in 2:size(f1, 1)
        cur = f1[i, j]
        if cur != 0
            left = f1[i-1, j]
            if left != 0
                cur = left
            else
                cur = clusterid
                clusterid += one(T)
            end
            f1[i, j] = cur
        end
    end
    for j in 2:size(f1, 2)
        for i in 2:size(f1, 1)
            cur = f1[i, j]
            if cur != zero(T)
                left = f1[i-1, j]
                top = f1[i, j-1]
                if left != 0 && top != 0
                    cur = left
                    if left != top
                        push!(edges, (top, cur))
                        #push!(edges, (cur, top))
                    end
                elseif top != 0
                    cur = top
                elseif left != 0
                    cur = left
                else
                    cur = clusterid
                    clusterid += one(T)
                end
                f1[i, j] = cur
            end
        end
    end
    clusterid, edges
end
function compresskeys(edges, vertexindex, vertexindex2)
    key = 1
    for (a,b) in edges
        for k in (a,b)
            if !haskey(vertexindex, k)
                vertexindex[k] = key
                push!(vertexindex2, k)
                @assert vertexindex2[key] == k
                key += 1
            end
        end
    end
end
function mergeclusters(edges::Set{Tuple{T, T}}) where T
    vertexindex = Dict{T, Int}() # used to compress arbitrary keys to 1:maxid
    vertexindex2 = Vector{T}() # 1:maxid to old keys
    compresskeys(edges, vertexindex, vertexindex2) # to use with SimpleGraph
    nnodes = length(vertexindex2)
    # Maybe there is some more complicated graph that supports Int32
    g = SimpleGraph(nnodes)
    for (a,b) in edges
        add_edge!(g, vertexindex[a], vertexindex[b])
    end
    components = connected_components(g)

    # Make a dictionary (old key => merged clusterid)
    dict = Dict{T, T}( 0 => 0)
    key = one(T)
    for set in components
        for elem in set
            dict[vertexindex2[elem]] = key
        end
        key += one(T)
    end
    dict, key
end
function changeclusters!(f1, dict, newclusterstart::T) where T
    key = newclusterstart
    for i in eachindex(f1)
        # f1[i] = get!(dict, f1[i]) do
        #     key += one(T)
        #     key
        # end
        f1i = f1[i]
        # Accessing the dict is costly
        # Usually there are a lot of zeros and then there is no need for the dict.
        # Should maybe change dict to Vector
        # Or Vector for small + dict for others?
        if f1i != 0
            if !haskey(dict, f1i)
                key += one(T)
                dict[f1i] = key
            end
            f1[i] = dict[f1i]
        end
    end
    key
end


text:
using FastTransforms, Makie, Random

include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")
include("sphcesaro.jl")

Random.seed!(0)

n = 64
U0 = zeros(4n, 8n-1), zeros(4n, 8n-1)
U0[1][1:n, 1:2n-1] = sphrandn(Float64, n, 2n-1)/n
U0[2][1:n, 1:2n-1] = sphrandn(Float64, n, 2n-1)/3n
U0[2][1] += 0.5sqrt(4π)
n = size(U0[1], 1)

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = Float32[cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = Float32[sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = Float32[cospi(θ) for θ in θ, φ in φ]

P = plan_sph2fourier(U0[1])
PS = plan_sph_synthesis(U0[1])
F = deepcopy(U0[1])
Ft = zeros(Float32, size(F, 1)+2, size(F, 2)+1)
lmul!(P, F)
lmul!(PS, F)
fillF!(Ft, F)

G = deepcopy(U0[2])
Gt = zeros(Float32, size(G, 1)+2, size(G, 2)+1)
lmul!(P, G)
lmul!(PS, G)
fillF!(Gt, G)
println("The extrema of v: ", extrema(G))

scene = Scene(resolution = (1200, 1200));
surf = surface!(scene, x, y, z, color = Ft, colormap = :viridis, colorrange = extrema(Ft))[end];
update_cam!(scene, Vec3f0(2.5), Vec3f0(0), Vec3f0(0, 0, 1))
scene.center = false
scene
Makie.save("plots/LGMu0.jpeg", scene)
Makie.save("plots/NGMu0.jpeg", scene)

T = 4.0
n = round(Int, 10*T)

SLPDE = SemiLinearPDE((LaplaceBeltrami(0.01),LaplaceBeltrami(0.1)), (NonlinearOperator((u,v)->(u-v)*u/v),NonlinearOperator((u,v)->u^2-v)), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io; colorrange = :notfixed)
Makie.save("plots/LGM.gif", io)

SLPDE = SemiLinearPDE((NonlocalLaplaceBeltrami(0.01, 0.0, 1.0),NonlocalLaplaceBeltrami(0.1, -0.5, 0.05)), (NonlinearOperator((u,v)->(u-v)*u/v),NonlinearOperator((u,v)->u^2-v)), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io; colorrange = :notfixed)
Makie.save("plots/NGM.gif", io)


text:
function Base.print(io::IO, scheduler::AbstractScheduler)
    print(io, nameof(typeof(scheduler)), " with ")
    nworkers = length(scheduler.workers)
    print(io, nworkers, " worker", nworkers == 1 ? "" : "s")
end

function Base.show(io::IO, ::MIME"text/plain", scheduler::AbstractScheduler)
    print(io, scheduler)
end

function Base.show(io::IO, ::MIME"text/plain", multiq::MultiQueue)
    print(io, "MultiQueue with ", length(multiq.heaps), " heaps and")
    print(io, " #elements = ")
    isfirst = true
    for h in multiq.heaps
        if isfirst
            isfirst = false
        else
            print(io, ", ")
        end
        if trylock(h)
            n = try
                length(h.heap)
            finally
                unlock(h)
            end
            print(io, n)
        else
            print(io, "???")
        end
    end
end