text:
# by lazarusA
using Legendre, GLMakie
using GeometryBasics, LinearAlgebra, StatsBase
using Makie: get_dim, surface_normals
ϵ = 1e-7
# thanks to @jkrumbiegel for the lift
let
    function Y(θ, ϕ, l, m)
        if m < 0
            return (-1)^m * √2 * Nlm(l, abs(m)) * Plm(l, abs(m), cos(θ)) * sin(abs(m)*ϕ)
        elseif m == 0
            return sqrt((2*l+1)/4π)*Plm(l, m, cos(θ))
        else
            return (-1)^m * √2 * Nlm(l, m) * Plm(l, m, cos(θ)) * cos(m*ϕ)
        end
    end
    function getMesh(x,y,z)
        positions = vec(map(CartesianIndices(z)) do i
        GeometryBasics.Point{3, Float32}(
            get_dim(x, i, 1, size(z)),
            get_dim(y, i, 2, size(z)),
            z[i])
        end)
        faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
        normals = surface_normals(x, y, z)
        vertices = GeometryBasics.meta(positions; normals=normals)
        meshObj = GeometryBasics.Mesh(vertices, faces)
        meshObj
    end
    function Π(v, n)
        v - dot(v, n) / norm(n)^2 * n
    end
    function multiply(q₁, q₂)
        z₁ = Complex(q₁[1], q₁[2])
        z₂ = Complex(q₁[3], q₁[4])
        A = [z₁ -conj(z₂); z₂ conj(z₁)]
        w₁ = Complex(q₂[1], q₂[2])
        w₂ = Complex(q₂[3], q₂[4])
        B = [w₁ -conj(w₂); w₂ conj(w₁)]
        C = A * B
        x₁, x₂ = C[1,1], C[2,1]
        [real(x₁); imag(x₁); real(x₂); imag(x₂)]
    end
    d(Φ, p) = begin
        rows = length(p)
        D = Array{Float64,1}(undef, rows)
        E = Array{Float64,2}(I, rows, rows) .* ϵ
        P = repeat(p, 1, rows)
        P′ = P + E
        for i in 1:rows
            D[i] = (Φ(P′[:, i]) - Φ(P[:, i])) ./ ϵ
        end
        D
    end
    directional(Φ, p, x) = begin
        rows = length(x)
        D = Array{Float64,1}(undef, rows)
        E = Array{Float64,2}(I, rows, rows) .* ϵ
        P = repeat(x, 1, rows)
        P′ = P + E
        for i in 1:rows
            D[i] = (Φ(p, P′[:, i]) - Φ(p, P[:, i])) ./ ϵ
        end
        D
    end
    getz₀(longitude, latitude, phase) = ℯ^(im * (longitude + phase)) * √((1 - sin(latitude)) / 2)
    getz₁(longitude, latitude, phase) = ℯ^(im * phase) * √((1 + sin(latitude)) / 2)
    TS³(z₀, z₁) =  begin
        # initialization
        X₀ = Complex(rand() + im * rand())
        X₁ = Complex(rand() + im * rand())
        getloss(z, x) = begin
            z₀ = Complex(z[1] + im * z[2])
            z₁ = Complex(z[3] + im * z[4])
            X₀ = Complex(x[1] + im * x[2])
            X₁ = Complex(x[3] + im * x[4])
            abs(conj(z₀) * X₀ + conj(z₁) * X₁)
        end
        x = [real(X₀); imag(X₀); real(X₁); imag(X₁)]
        z = [real(z₀); imag(z₀); real(z₁); imag(z₁)]
        loss = getloss(z, x)
        threshold = 1e-4
        η = 1e-2
        iterations = 150
        i = 0
        while loss > threshold && i < iterations
            g = directional(getloss, z, x)
            X₀ = Complex(x[1] - η * g[1] + im * (x[2] - η * g[2]))
            X₁ = Complex(x[3] - η * g[3] + im * (x[4] - η * g[4]))
            x = [real(X₀); imag(X₀); real(X₁); imag(X₁)]
            loss = getloss(z, x)
            i += 1
        end
        X₀, X₁
    end
    connection1form(z₀, z₁) = begin
        X₀, X₁ = TS³(z₀, z₁)
        1/2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁))
    end
    tangent_tail(x, f) = [x[1]; x[2]; f(x)]
    tangent_head(x, f, ξ) = begin
        dΦ = d(f, x)
        ξΦ = dot(dΦ, ξ)
        dξΦ = dΦ .* ξ
        [dξΦ[1]; dξΦ[2]; (f([x[1] + dξΦ[1] * ϵ; x[2] + dξΦ[2] * ϵ]) - f(x)) ./ ϵ] .* ξΦ
    end
    f₁(i) = begin
        √(1 - (i[1])^2 - (i[2])^2)
    end
    f₂(i) = -f₁(i)

    segments = 200
    radius = 0.05
    # Grids of polar and azimuthal angles
    θ = LinRange(0, π, segments)
    ϕ = LinRange(0, 2π, segments)
    xx = [sin(θ)*sin(ϕ) for θ in θ, ϕ in ϕ]
    yy = [sin(θ)*cos(ϕ) for θ in θ, ϕ in ϕ]
    zz = [cos(θ)        for θ in θ, ϕ in ϕ]

    l = Node(4)
    m = Node(1)

    lon = Node(1.0)
    lat = Node(1.0)
    ang = Node(1.0)
    vector = Node([1.0; 0.0; 0.0])

    #w₀ = @lift(ℯ^(im * $ang) * sqrt((1 + sin($lat)) / 2))
    #w₁ = @lift(ℯ^(im * ($lon + $ang)) * sqrt((1 - sin($lat)) / 2))
    #w = @lift($w₁ / $w₀)
    w = lift(lon, lat) do lon, lat
        a₁ = cos(lat) * cos(lon)
        a₂ = cos(lat) * sin(lon)
        a₃ = sin(lat)
        a₁ / (1 - a₃) + im * a₂ / (1 - a₃)
    end
    fiber = lift(lon, lat, ang) do lon, lat, ang
        points = [[lon; lat] + [cos(i); sin(i)] .* radius for i in range(0, stop = 2π, length = segments)]
        ψ = range(0, stop = ang, length = segments)
        xyz(ψ, point) = begin
            v₀ = getz₀(point[1], point[2], ψ)
            v₁ = getz₁(point[1], point[2], ψ)
            q₁, q₂ = real(v₀), imag(v₀)
            q₃, q₄ = real(v₁), imag(v₁)
            q = [q₁; q₂; q₃; q₄]
            p = q[1:3] ./ (1 - q₄)
            magnitude = norm(p)
            normalize(p) * tanh(magnitude)
        end
        a = [xyz(i, j) for i in ψ, j in points]
        a₁ = map(i -> i[1], a)
        a₂ = map(i -> i[2], a)
        a₃ = map(i -> i[3], a)
        a₁, a₂, a₃
    end
    fiber₁ = @lift($fiber[1])
    fiber₂ = @lift($fiber[2])
    fiber₃ = @lift($fiber[3])

    fiber2 = lift(lon, lat, ang) do lon, lat, ang
        points = [[lon; lat] + [cos(i); sin(i)] .* radius for i in range(0, stop = 2π, length = segments)]
        ψ = range(ang, stop = 2π, length = segments)
        xyz(ψ, point) = begin
            v₀ = getz₀(point[1], point[2], ψ)
            v₁ = getz₁(point[1], point[2], ψ)
            q₁, q₂ = real(v₀), imag(v₀)
            q₃, q₄ = real(v₁), imag(v₁)
            q = [q₁; q₂; q₃; q₄]
            p = q[1:3] ./ (1 - q₄)
            magnitude = norm(p)
            normalize(p) * tanh(magnitude)
        end
        a = [xyz(i, j) for i in ψ, j in points]
        a₁ = map(i -> i[1], a)
        a₂ = map(i -> i[2], a)
        a₃ = map(i -> i[3], a)
        a₁, a₂, a₃
    end
    fiber2₁ = @lift($fiber2[1])
    fiber2₂ = @lift($fiber2[2])
    fiber2₃ = @lift($fiber2[3])

    tail = lift(w) do w
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        tangent_tail([u; v], f)
    end
    xconst = lift(w, tail) do w, tail
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [1; 0]))
    end
    yconst = lift(w, tail) do w, tail
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [0; 1]))
    end
    head = lift(w, ang, tail, xconst, yconst) do w, ang, tail, x, y
        # parallel transport
        #v = to_value(vector)
        #n = normalize(cross(x, y))
        #if -1 < dot(x, y) < 1
        #    projected = Π(v, n)
        #else
        #    projected = v
        #end
        #projected = normalize(projected)
        #vector[] = projected
        #projected
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [cos(ang); sin(ang)]))

        #p₁ = [lon; lat]
        #p₂ = p₁ + [0.1; 0.1]
        #v₀ = ℯ^(im * (p₁[1] + ψ)) * sqrt((1 - sin(p₁[2])) / 2)
        #v₁ = ℯ^(im * ψ) * sqrt((1 + sin(p₁[2])) / 2)

    end
    connectiontail = lift(lon, lat, ang) do lon, lat, ang
        point = [lon; lat]
        v₀ = getz₀(point[1], point[2], ang)
        v₁ = getz₁(point[1], point[2], ang)
        q₁, q₂ = real(v₀), imag(v₀)
        q₃, q₄ = real(v₁), imag(v₁)
        q = [q₁; q₂; q₃; q₄]
        p = q[1:3] ./ (1 - q₄)
        magnitude = norm(p)
        normalize(p) * tanh(magnitude)
    end
    connectionhead = lift(lon, lat, ang) do lon, lat, ang
        point = [lon; lat]
        v₀ = getz₀(point[1], point[2], ang)
        v₁ = getz₁(point[1], point[2], ang)
        oneform = connection1form(v₀, v₁)
        v₀, v₁ = oneform * v₀, oneform * v₁
        q₁, q₂ = real(v₀), imag(v₀)
        q₃, q₄ = real(v₁), imag(v₁)
        q = [q₁; q₂; q₃; q₄]
        p = q[1:3] ./ (1 - q₄)
        magnitude = norm(p)
        h = normalize(p) * tanh(magnitude)
        t = to_value(connectiontail)
        h - t
    end

    tails₁ = @lift([$tail[1], $tail[1], $tail[1]])
    tails₂ = @lift([$tail[2], $tail[2], $tail[2]])
    tails₃ = @lift([$tail[3], $tail[3], $tail[3]])
    heads₁ = @lift([$head[1], $xconst[1], $yconst[1]])
    heads₂ = @lift([$head[2], $xconst[2], $yconst[2]])
    heads₃ = @lift([$head[3], $xconst[3], $yconst[3]])

    conn_tail₁ = @lift([$connectiontail[1]])
    conn_tail₂ = @lift([$connectiontail[2]])
    conn_tail₃ = @lift([$connectiontail[3]])
    conn_head₁ = @lift([$connectionhead[1]])
    conn_head₂ = @lift([$connectionhead[2]])
    conn_head₃ = @lift([$connectionhead[3]])

    x₁ = @lift(xx .* 0.1 .+ $tail[1])
    y₁ = @lift(yy .* 0.1 .+ $tail[2])
    z₁ = @lift(zz .* 0.1 .+ $tail[3])

    ambient =  Vec3f0(0.75, 0.75, 0.75)
    cmap = (:dodgerblue, :white) # how to include this into menu options?
    with_theme(theme_black()) do
        fig = Figure(resolution = (1400, 800))
        menu = Menu(fig, options = ["Spectral_11", "viridis", "heat", "plasma", "magma", "inferno"])
        Ygrid = lift(l, m) do l, m
            [Y(θ, ϕ, l, m) for θ in θ, ϕ in ϕ]
        end
        Ylm = @lift(abs.($Ygrid))
        Ygrid2 = @lift(vec($Ygrid))

        ax1 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        ax2 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        pltobj1 = mesh!(ax1, getMesh(xx, yy, zz), color = Ygrid2, colormap = cmap, ambient = ambient)
        pltobj3 = mesh!(ax1, @lift(getMesh($x₁, $y₁, $z₁)), color = :gold, colormap = cmap, ambient = ambient)
        pltobj4 = arrows!(ax1, tails₁, tails₂, tails₃, heads₁, heads₂, heads₃, color = [:red, :silver, :gray], lengthscale = 0.3f0)
        pltobj5 = mesh!(ax2, @lift(getMesh($fiber₁, $fiber₂, $fiber₃)), color = :red, colormap = cmap, ambient = ambient)
        pltobj6 = mesh!(ax2, @lift(getMesh($fiber2₁, $fiber2₂, $fiber2₃)), color = Ygrid2, colormap = cmap, ambient = ambient)
        pltobj7 = arrows!(ax2, conn_tail₁, conn_tail₂, conn_tail₃, conn_head₁, conn_head₂, conn_head₃, color = [:gold], lengthscale = 1.0f0)
        cbar = Colorbar(fig, pltobj1, label = "Yₗₘ(θ,ϕ)", width = 11, tickalign = 1, tickwidth = 1)
        fig[1,1] = ax1
        fig[1,2] = ax2
        fig[1,3] = cbar
        round3(n) = round(n, sigdigits = 3)
        degree(n) = round3(180(n / π))
        fig[0,1:2] = Label(fig, @lift("Hopf fiber: longitude = $(degree($lon)), latitude = $(degree($lat)), S¹action = $(degree($ang))."), textsize = 20)
        fig[2, 0] = vgrid!(
            Label(fig, "Colormap", width = nothing),
            menu; tellheight = false, width = 100)
        on(menu.selection) do s
            pltobj1.colormap = s
            pltobj3.colormap = s
            pltobj4.colormap = s
            pltobj5.colormap = s
            pltobj6.colormap = s
            pltobj7.colormap = s
        end
        sl3 = Slider(fig[end+1, 2:3], range = range(-π, stop = π, length = 100))
        sl4 = Slider(fig[end+1, 2:3], range = range(-π/2, stop = 0.99 * π/2, length = 100))
        sl5 = Slider(fig[end+1, 2:3], range = range(0, stop = 2π, length = 100))
        set_close_to!(sl3, to_value(lon))
        set_close_to!(sl4, to_value(lat))
        set_close_to!(sl5, to_value(ang))
        connect!(lon, sl3.value)
        connect!(lat, sl4.value)
        connect!(ang, sl5.value)
        tight_ticklabel_spacing!(cbar)
        display(fig)
        framerate = 30
        totaltime = 15
        timestamps = range(0, totaltime, step=1/framerate)

#=         record(fig, joinpath("gallery", "connection1form.mp4"), timestamps;
            framerate = framerate) do t
            p = t / totaltime
            longitude = sin(p * 2π - π)
            latitude = cos(p * π - π/2)
            action = p * 2π
            set_close_to!(sl3, longitude)
            set_close_to!(sl4, latitude)
            set_close_to!(sl5, action)
        end =#
    end
end


text:
mutable struct Trie
    children::Dict{String, Trie}
    terminal::Bool
    metadata::Dict{String, String}

    function Trie()
        self = new()
        self.children = Dict{String, Trie}()
        self.terminal = false
        self
    end
end

function add_to_trie(trie::Trie, words::Array{SubString{String},1}, metadata::Dict{String, String})
    if length(words) > 0
        word = words[1]
        if haskey(trie.children, word)
            add_to_trie(trie.children[word], words[2:end], metadata)
        else
            trie.children[word] = Trie()
            add_to_trie(trie.children[word], words[2:end], metadata)
        end
    else
        trie.terminal = true
        trie.metadata = metadata
    end
    return trie
end

function add_to_trie(trie::Trie, phrase::String, metadata::Dict{String, String})
    return add_to_trie(trie, split(phrase, " "), metadata)
end

ext: jl


text:
# primer81.jl
# inicjalizujemo podatke u grafu
include("labyrinth.jl")
include("primer81BFS.jl")

nbRows = 5;
nbColumns = 5;
G = emptyLabyrinth(nbRows, nbColumns)

defineNodesType!(G,[1 2 3 5 6 8 10 11 12 13 14 15 17 21 22 23 24 25],'P'); #PUT
defineNodesType!(G, [4 7 9 16 18 19 20], 'Z'); #ZID

println("Izgled lavirinta pre obilaska")
println()
printLabyrinthType(G)
println()
println("Ispis lavirinta pre obilaska")
println()
printLabyrinthColor(G)
println()

println("Ispis lavirinta posle obilaska")
println()
BFS!(G, 21);
printLabyrinthColor(G)
println()

println("Izgled lavirinta sa putanjom")
println()
path = getPath(G, 21, 5)
printLabyrinthPath(G, path)
println()


text:
using Distributions

module IsingModels

using Distributions

export Ising

type Ising <: Distribution
    state::Array{Int8}
    isferromagnetic::Bool
    isperiodic::Bool
    temperature::Float64
    applied_field::Array{Int8}
    iterations::Int
end

function Ising(n::Integer;
               m::Integer = -1,
               isferromagnetic::Bool = true,
               isperiodic::Bool = true,
               temperature::Real = 1.0,
               iterations::Integer = -1,
               applied_field::Matrix = Array(Int8, 0, 0))

    if m == -1
        m = n
    end

    if iterations == -1
        iterations = n * n * 1_000
    end

    state = Array(Int8, n, m)
    fill!(state, -1)

    Ising(state,
          isferromagnetic,
          isperiodic,
          temperature,
          applied_field,
          iterations)
end

function determine_j(d::Ising)
    if d.isferromagnetic
        return int8(1)
    else
        return int8(-1)
    end
end

function determine_local_field(d::Ising, r::Integer, c::Integer, n::Integer, m::Integer, j::Int8)
    # Sum neighbors
    s = int8(0)

    # Upper neighbor
    if r > 1
        s += d.state[r - 1, c] * j
    else
        if d.isperiodic
            s += d.state[n, c] * j
        end
    end

    # Lower neighbor
    if r + 1 <= n
        s += d.state[r + 1, c] * j
    else
        if d.isperiodic
            s += d.state[1, c] * j
        end
    end

    # Left neighbor
    if c > 1
        s += d.state[r, c - 1] * j
    else
        if d.isperiodic
            s += d.state[r, m] * j
        end
    end

    # Right neighbor
    if c + 1 <= m
        s += d.state[r, c + 1] * j
    else
        if d.isperiodic
            s += d.state[r, c] * j
        end
    end

    # Apply external field
    if !isempty(d.applied_field)
        s += d.applied_field[r, c]
    end

    return s
end

function flipbit!(d::Ising, n::Integer, m::Integer, j::Int8)
    r, c = rand(1:n), rand(1:m)

    b = determine_local_field(d, r, c, n, m, j)

    p = 1 / (1 + exp(-2 * b / d.temperature))

    if !isfinite(p)
        d.state[r, c] = int8(1)
    else
        if rand() < p
            d.state[r, c] = int8(1)
        else
            d.state[r, c] = int8(-1)
        end
    end

    return
end

function Distributions.rand!(d::Ising)
    # Determine constants
    j = determine_j(d)
    n, m = size(d.state)

    # Randomly initialize state
    for r in 1:n
        for c in 1:m
            if rand() < 0.5
                d.state[r, c] = int8(-1)
            else
                d.state[r, c] = int8(1)
            end
        end
    end

    # Flip bits for a fixed number of iterations
    for iteration in 1:d.iterations
        flipbit!(d, n, m, j)
    end

    return
end

function Distributions.rand(d::Ising)
    rand!(d)
    return copy(d.state)
end

Distributions.mean(d::Ising) = NaN
Distributions.var(d::Ising) = NaN

end # module


text:
export ExploringStartPolicy

using Random
using MacroTools: @forward

mutable struct ExploringStartPolicy{P<:AbstractPolicy, A, R<:AbstractRNG} <: AbstractPolicy
    policy::P
    actions::A
    is_start::Bool
    rng::R
end

ExploringStartPolicy(;policy, actions, is_start=true,seed=nothing) = ExploringStartPolicy(policy, actions, is_start, MersenneTwister(seed))

@forward ExploringStartPolicy.policy RLBase.get_prob, RLBase.update!

(p::ExploringStartPolicy)(obs) = p.is_start ? rand(p.rng, p.actions) : p.policy(obs)

function (agent::Agent{<:ExploringStartPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::PreEpisodeStage,
    obs,
)
    empty!(agent.trajectory)
    agent.policy.is_start = true
    nothing
end

function (agent::Agent{<:ExploringStartPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::PostActStage,
    obs,
)
    push!(agent.trajectory; reward = get_reward(obs), terminal = get_terminal(obs))
    agent.policy.is_start = false
    nothing
end


text:
module LifeTable

using DataFrames, LsqFit, GLM

export periodlifetable, causelife, mortsurvfit, predict, mortsurvparamsfit

function periodlifetable(inframe, sex, openend=true, intype="count")
	if sex == 1
		aint = 0.045
		acoef = 2.684
		aconst = 0.330
	elseif sex == 2
		aint = 0.053
		acoef = 2.8
		aconst = 0.350
	end

	nrows = size(inframe, 1)
	age = inframe[1]

	if openend
		i = age[2:nrows] .- age[1:nrows-1]
	else
		ifirst = age[2:nrows] .- age[1:nrows-1]
		iend = ifirst[nrows-1]
		i = [ifirst; iend]
	end

	if intype == "count"
		m = inframe[3]./inframe[2]
	elseif intype == "rate"
		m = inframe[2]
	end

	if age[1] == 0
		if m[1] >= 0.107
			a0 = aconst
		else
			a0 = aint + acoef * m[1]
		end
	else
		a0 = 0.5
	end

	q(m, a, i) = (i*m) / (1+(1-a)*i*m)
	ld(l, a, d, i) = (l - (1-a) * d) * i

	if openend
		aend = 1 / m[nrows]
		aoth = fill(0.5, nrows-2)
		aclosed = [a0; aoth]
		a = [aclosed; aend]
		qclosed = map(q, m[1:nrows-1], aclosed, i)
		q = [qclosed; 1]
	else
		aoth = fill(0.5, nrows-1)
		a = [a0; aoth]
		q = map(q, m, a, i)
	end

	p = 1 .- q
	l = cumprod([1; p[1:nrows-1]])
	d = q .* l

	if openend
		ldclosed = map(ld, l[1:nrows-1], aclosed, d[1:nrows-1], i)
		ld = [ldclosed; l[nrows]*aend]
		revl = reverse(ld)
		tclosed = reverse(cumsum(revl))[1:nrows-1]
		t = [tclosed; ld[nrows]]
	else
		ld = map(ld, l, a, d, i)
		revl = reverse(ld)
		t = reverse(cumsum(revl))
	end

	e = t./l

	outframe = DataFrame(age = age, m = m, a = a, q = q, p = p, l = l, d = d,
	ld = ld, t = t, e = e)

	return outframe

end

function causelife(lifetable, causefreq)
	dca = causefreq .* lifetable[:d]
	l0ca = reverse(cumsum(reverse(dca)))
	fca = l0ca ./ lifetable[:l]
	mca = (causefreq .* lifetable[:m])
	mcacoh = mca ./ fca

	outframe = DataFrame(age = lifetable[:age], m = mcacoh, f = fca)
	return outframe
end

gomprate(x,p) = p[1] * exp.(p[2].*x)
gompsurv(x,p) = exp.(-(p[1]/p[2]) .* (exp.(p[2].*x)-1))
weibrate(x,p) = (p[1]/p[2]) .* (x./p[2]).^(p[1]-1)
weibsurv(x,p) = exp.(-(x./p[2]).^p[1])
# Functions for blocks with redundant components with exponentially distributed lifespans,
# see Gavrilov and Gavrilova, Journal of Theoretical Biology (2001).
gavrblrate(x,p) = (p[1]*p[2]*exp.(-p[2].*x).*(1-exp.(-p[2].*x)).^(p[1]-1))./
	(1-(1-exp.(-p[2].*x)).^p[1])
gavrblsurv(x,p) = 1-(1-exp.(-p[2].*x)).^p[1]

models = Dict(
	"gompertz" => Dict("rate" => gomprate, "surv" => gompsurv,
		"p" => [exp(-18); 0.1]),
	"weibull" => Dict("rate" => weibrate, "surv" => weibsurv,
		"p" => [10.0; 80.0]),
	"gavrblock" => Dict("rate" => gavrblrate, "surv" => gavrblsurv,
		"p" => [10.0; 0.03]),
	)

function mortsurvfit(lifetable, numbdeaths, func, functype)
	ycols = Dict("rate" => :m, "surv" => :l)

	xarr = convert(Array{Int}, lifetable[:age])
	yarr = convert(Array{Float64}, lifetable[ycols[functype]])
	warr = sqrt.(convert(Array{Int}, numbdeaths))
	model = models[func][functype]

	fit = curve_fit(model, xarr, yarr, warr, models[func]["p"])
	return Dict("fit" => fit, "func" => func, "functype" => functype, "model" => model)
end

predict(x, msfit) = msfit["model"](x, msfit["fit"].param)

gompratetr_x(p1, p2) = p2
gompratetr_y(p1, p2) = log(p1)
gompsurvtr_x(p1, p2) = p2
gompsurvtr_y(p1, p2) = log(p1/p2)
weibratetr_x(p1, p2) = p1 - 1
weibratetr_y(p1, p2) = log(p1/p2) - (p1-1) * log(p2)
weibsurvtr_x(p1, p2) = p1
weibsurvtr_y(p1, p2) = -p1 * log(p2)

transmodels = Dict(
	"gompertz" =>  Dict("rate" => [gompratetr_x; gompratetr_y],
		"surv" => [gompsurvtr_x; gompsurvtr_y]),
	"weibull" =>  Dict("rate" => [weibratetr_x; weibratetr_y],
		"surv" => [weibsurvtr_x; weibsurvtr_y]),
	)

function mortsurvparamsfit(msfits)
	p1arr = Float64[]
	p2arr = Float64[]

	for msfit in msfits
		push!(p1arr, msfit["fit"].param[1])
		push!(p2arr, msfit["fit"].param[2])
	end

	func = msfits[1]["func"]
	functype = msfits[1]["functype"]

	trans_x = transmodels[func][functype][1]
	trans_y = transmodels[func][functype][2]

	xarr = map(trans_x, p1arr, p2arr)
	yarr = map(trans_y, p1arr, p2arr)
	trans_params = DataFrame(p1 = p1arr, p2 = p2arr, X = xarr, Y = yarr)

	fit = glm(@formula(Y~X), trans_params, Normal(), IdentityLink())
	return Dict("trans_params" => trans_params, "fit" => fit)
end

end


text:
# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Synthetic (Tetrad SDE) setup][Synthetic (Tetrad SDE) setup:1]]
function generateSyntheticTetradLearningProblem(data)
    #--------------------begin problem setup--------------------#
    trainingData = data
    initialConditions = trainingData[:,1,:];

    #split into training/test
    validationData = trainingData[:,:,:];
    trainingData = trainingData[:,:,:];
    numDims, numSteps, numTrajs = size(trainingData);
    dt = 3e-4;
    tsteps = Array{Float64, 1}(range(0.0, step=dt, length=numSteps));

    x,y = parseTrainingData(trainingData, dt);
    miniBatchSize = 1000; #see HyperParameterTuning for justification
    maxiters = 5000; #plateau's around 4k
    #optimizer = Flux.Optimise.Optimiser(ExpDecay(0.001, 0.5, 1000, 1e-7), ADAM(0.001));
    optimizer = ADAM(0.01);
    # see https://fluxml.ai/Flux.jl/stable/training/optimisers/

    hiddenDim = 80;
    θ_f, drift_nn = Flux.destructure(Chain(Dense(18,hiddenDim,tanh),
					   Dense(hiddenDim,9)));

    function encode(u,p,args...)
	return u;
    end
    function decode(u,p,args...)
	return u;
    end

    function trueDrift(u,p,t)
	return nothing;
    end

    function truePrecisionMat(u,p,t)
	return nothing;
    end

    function trueDiff(u, p, t)
	return nothing;
    end

    function parameterizedDriftClosure(u,p)
	return drift_nn(p)(u);
    end

    function parameterizedPrecisionMat(u,p)
	return nothing;
    end

    driftLength = length(θ_f);
    diffLength = 0;

    driftStart = 1;
    if (driftLength > 0)
	driftEnd = driftStart + (driftLength-1);
	diffStart = driftEnd+1;
    else
	driftEnd = 0;
	diffStart = 1;
    end
    diffEnd = diffStart + (diffLength-1);

    function modelDrift(u,p,t)
	numSamples = size(u)[end];
	ρ = [reshape(u[1:9,i], (3,3)) for i in 1:numSamples];
	M = [reshape(u[10:18,i], (3,3)) for i in 1:numSamples];
	closure = parameterizedDriftClosure(u,p);
	placeholder = [ reshape(-M[i]^2 + reshape(closure[:,i],(3,3)), 9) for i in 1:numSamples];
	return hcat(placeholder...);
    end

    function fullSystemDrift(u,p,t)
	numTrajs = size(u)[2];
	vgtPiece = modelDrift(u,p,t);
	rhoPiece = zeros(9,numTrajs);
	for i in 1:numTrajs
	    M = reshape(u[10:18,i], (3,3));
	    rho = reshape(u[1:9,i], (3,3));
	    rhoPiece[:,i] = reshape(M'*rho, 9);
	end
	return cat(rhoPiece, vgtPiece, dims=1);
    end

    function createNoiseRateForcingHIT(D_a, D_s)
	c = zeros(3,3,3,3);
	delta(i,j) = (i==j) ? 1 : 0;
	for i in 1:3
	    for j in 1:3
		for k in 1:3
		    for l in 1:3
			term1 = -(1/3)*sqrt(D_s/5)*delta(i,j)*delta(k,l);
			term2 = (1/2)*(sqrt(D_s/5) + sqrt(D_a/3))*delta(i,k)*delta(j,l);
			term3 = (1/2)*(sqrt(D_s/5) - sqrt(D_a/3))*delta(i,l)*delta(j,k);
			c[i,j,k,l] = term1 + term2 + term3;
		    end
		end
	    end
	end
	return c;
    end

    c = reshape(createNoiseRateForcingHIT(15,15), (9,9));
    function modelDiff(u,p,t)
	return [Matrix(I,9,9)  zeros(9,9)
		zeros(9,9)     c];
    end

    #statically allocate for training
    driftGrad = zeros(driftLength);
    driftGradHolder = zeros(driftLength, miniBatchSize);
    predictionErrorGrad = zeros(diffLength);
    predictionErrorGradHolder = zeros(diffLength, miniBatchSize);
    detPiece = zeros(diffLength);
    detPieceHolder = zeros(diffLength, miniBatchSize);

    precisionMat = zeros(9,miniBatchSize);
    lossPerStep = zeros(miniBatchSize);
    debiasWeight = zeros(miniBatchSize);

    function predict_(x, p)
	#assume autonomous
	return prob.modelDrift_(x, p, 0.0);
    end

    function lossAndGrads(x,y,p,_drift,_precisionMat)
	eps = 64.0; #approximately median(\dot M)
	δ = _drift(x,p,0.0) .- y[10:18,:];

	for i in 1:miniBatchSize
	    debiasWeight[i] = 1.0/(norm(reshape(y[10:18,i], (3,3)))^2 + eps);
	    lossPerStep[i] = debiasWeight[i]*dot(δ[:,i],δ[:,i]);
	end
	sumLoss = (dt*eps/miniBatchSize)*sum(lossPerStep);

	#-----------------------------drift grad-----------------------------#

	if (driftLength > 0)
	    #x is state variable, w is weight vector, pullback(p->f(x,p), p)[2](w)[1] performs the tensor contraction
	    # (∂/∂p_i f_k) w^k
	    ∂f(x,w) = pullback(p->parameterizedDriftClosure(x,p), p[driftStart:driftEnd])[2](w)[1];
	    Threads.@threads for i in 1:miniBatchSize
		driftGradHolder[:,i] .= ∂f(x[:,i], debiasWeight[i]*δ[:,i]);
	    end
	    driftGrad = ((2*dt*eps)/miniBatchSize)*sum(driftGradHolder, dims=2);
	end

	#-----------------------------diff grad-----------------------------#
	diffGrad = zeros(diffLength);

	if (diffLength > 0)
	    #x is state variable, w is weight vector, pullback(p->f(x,p), p)[2](w)[1] performs the tensor contraction
	    # (∂/∂p_i f_k) w^k
	    ∂Π(x,w) = pullback(p_->_precisionMat(x,p_), p[diffStart:diffEnd])[2](w)[1];

	    #Threads.@threads for i in 1:miniBatchSize
	    for i in 1:miniBatchSize
		predictionErrorGradHolder[:,i] .= ∂Π(x[:,i], abs2.(δ[:,i]));
		detPieceHolder[:,i] .= ∂Π(x[:,i], precisionMat[:,i].^(-1));
	    end
	    predictionErrorGrad = dt*sum(predictionErrorGradHolder, dims=2);
	    detPiece = sum(detPieceHolder, dims=2);

	    diffGrad = (predictionErrorGrad - detPiece)/miniBatchSize;

	    #add in bias regularization gradient piece:
	    regularizationGrad = zeros(18);
	    for i in 1:size(regularizationGrad)[1]
		if ((-p[diffStart+17+i] + biasRegularizationOffset) > 0)
		    regularizationGrad[i] = -1.0;
		end
	    end
	    diffGrad[19:36] .+= biasRegularizationWeight*regularizationGrad;
	end

	grads = vcat(driftGrad, diffGrad);
	grads = reshape(grads, size(grads)[1]);
	return sumLoss,grads;
    end

    #covariance matrix is length(u)^2
    initialParams = [];
    if (driftLength > 0)
	initialParams = vcat(initialParams, θ_f)
    end
    if (diffLength > 0)
	initialParams = vcat(initialParams, ones(diffLength)); #currently using diagonal covariance
    end
    initialParams = Array{Float64}(initialParams);

    @assert length(initialParams) == (driftLength + diffLength)

    #--------------------end problem setup--------------------#

    lProb = LearningProblem(initialConditions,
			    trainingData,
			    x,
			    y,
			    validationData,
			    tsteps,
			    miniBatchSize,
			    maxiters,
			    optimizer,
			    trueDrift,
			    trueDiff,
			    encode,
			    decode,
			    modelDrift,
			    modelDiff,
			    fullSystemDrift,
			    nothing,
			    lossAndGrads,
			    initialParams,
			    driftLength,
			    diffLength,
			    parameterizedDriftClosure,
			    parameterizedPrecisionMat);

    println("Setup complete, beginning learning");

    return lProb;
end
# Synthetic (Tetrad SDE) setup:1 ends here

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
using Printf

function do_reduce()
    comm = MPI.COMM_WORLD

    MPI.Barrier(comm)

    root = 0
    r = MPI.Comm_rank(comm)

    sr = MPI.Reduce(r, MPI.SUM, root, comm)

    if(MPI.Comm_rank(comm) == root)
        @printf("sum of ranks: %s\n", sr)
    end
end


text:
## LEARNING NETWORK MACHINES

# ***
surrogate(::Type{<:Deterministic})  = Deterministic()
surrogate(::Type{<:Probabilistic})  = Probabilistic()
surrogate(::Type{<:Unsupervised}) = Unsupervised()
surrogate(::Type{<:Static}) = Static()

"""
    model_supertype(signature)

Return, if this can be deduced, which of `Deterministic`,
`Probabilistic` and `Unsupervised` is the appropriate supertype for a
composite model obtained by exporting a learning network with the
specified `signature`.

A learning network *signature* is a named tuple, such as
`(predict=yhat, transfrom=W)`, specifying what nodes of the network
are called to produce output of each operation represented by the
keys, in an exported version of the network.

If a supertype cannot be deduced, `nothing` is returned.

If the network with given `signature` is not exportable, this method
will not error but it will not a give meaningful return value either.

"""
function model_supertype(signature)

    operations = keys(signature)

    length(intersect(operations, (:predict_mean, :predict_median))) == 1 &&
        return Deterministic

    if :predict in operations
        node = signature.predict
        if node isa Source
            return Deterministic
        end
        if node.machine !== nothing
            model = node.machine.model
            model isa Deterministic && return Deterministic
            model isa Probabilistic && return Probabilistic
        end
    end

    return nothing

end

function machine(model::Surrogate, _sources::Source...; pair_itr...)

    # named tuple, such as `(predict=yhat, transform=W)`:
    signature = (; pair_itr...)
    for op in keys(signature)
        op in OPERATIONS || throw(ArgumentError(
            "`$op` is not an admissible operation. "))
    end
    if isempty(signature)
        model isa Supervised &&
            throw(ArgumentError(
                "You must specify at least `predict=<some node>`. "))
        model isa Unsupervised &&
            throw(ArgumentError(
                "You must specify at least `transform=<some node>`. "))
    end

    if model isa Supervised
        length(_sources) in [2, 3] ||
            error("Incorrect number of source nodes specified.\n"*
                  "Use  `machine(model, X, y; ...)` or "*
                  "`machine(model, X, y, w; ...)` when "*
                  "`model isa Supervised`. ")
    elseif model isa Unsupervised
        length(_sources) == 1 ||
            error("Incorrect number of source nodes specified.\n"*
                  "Use `machine(model, X; ...)` when "*
                  "`model isa Unsupervised. ` (even if `Static`). ")
    else
        throw(DomainError)
    end

    mach = Machine(model, _sources...)

    mach.fitresult = signature

    return mach

end

function machine(sources::Source...; pair_itr...)

    signature = (; pair_itr...)

    isempty(signature) &&
        throw(ArgumentError(
            "You must specify at least `predict=<some node>` or "*
            "or `transform=<some node>`. "))

    T = model_supertype(signature)
    if T == nothing
        @warn "Unable to infer surrogate model type. \n"*
        "Using Deterministic(). To override use:\n "*
        "`machine(Probabilistic(), ...)` or `machine(Interval(), ...)`"
        model = Deterministic()
    else
        model = surrogate(T)
    end

    return machine(model, sources...; pair_itr...)

end

"""
    anonymize!(sources)

Returns a named tuple `(sources=..., data=....)` whose values are the
provided source nodes and their contents respectively, and clears the
contents of those source nodes.

"""
function anonymize!(sources)
    data = Tuple(s.data for s in sources)
    [MLJBase.rebind!(s, nothing) for s in sources]
    return (sources=sources, data=data)
end

"""
    fit!(mach::Machine{<:Surrogate};
         rows=nothing,
         acceleration=CPU1(),
         verbosity=1,
         force=false))

Train the complete learning network wrapped by the machine
`mach`.

More precisely, if `s` is the learning network signature used to
construct `mach`, then call `fit!(N)`, where `N = glb(values(s)...)`
is a greatest lower bound on the nodes appearing in the signature. For
example, if `s = (predict=yhat, transform=W)`, then call
`fit!(glb(yhat, W))`. Here `glb` is `tuple` overloaded for nodes.

See also [`machine`](@ref)

"""
function fit!(mach::Machine{<:Surrogate}; kwargs...)

    signature = mach.fitresult
    glb_node = glb(values(signature)...) # greatest lower bound node
    fit!(glb_node; kwargs...)

    # mach.cache = anonymize!(mach.args)
    mach.state += 1
    mach.report = report(glb_node)
    return mach

end

# make learning network machines callable for use in manual export of
# learning networks:
function (mach::Machine{<:Surrogate})()
    # anonymize sources:
    mach.cache = anonymize!(mach.args)
    return mach.fitresult, mach.cache, mach.report
end

MLJModelInterface.fitted_params(mach::Machine{<:Surrogate}) =
    fitted_params(glb(values(mach.fitresult)...))

## DUPLICATING AND REPLACING PARTS OF A LEARNING NETWORK MACHINE

"""
    replace(mach, a1=>b1, a2=>b2, ...; empty_unspecified_sources=false)

Create a deep copy of a learning network machine `mach` but replacing
any specified sources and models `a1, a2, ...` of the original
underlying network with `b1, b2, ...`.

If `empty_unspecified_sources=true` then any source nodes not
specified are replaced with empty source nodes.

"""
function Base.replace(mach::Machine{<:Surrogate},
                      pairs::Pair...; empty_unspecified_sources=false)

    signature = mach.fitresult
    interface_nodes = values(signature)

    W = glb(interface_nodes...)

    # Note: We construct nodes of the new network as values of a
    # dictionary keyed on the nodes of the old network. Additionally,
    # there are dictionaries of models keyed on old models and
    # machines keyed on old machines. The node and machine
    # dictionaries must be built simultaneously.

    # build model dict:
    model_pairs = filter(collect(pairs)) do pair
        first(pair) isa Model
    end
    models_ = models(W)
    models_to_copy = setdiff(models_, first.(model_pairs))
    model_copy_pairs = [model=>deepcopy(model) for model in models_to_copy]
    newmodel_given_old = IdDict(vcat(model_pairs, model_copy_pairs))

    # build complete source replacement pairs:
    sources_ = mach.args # sources(W)
    specified_source_pairs = filter(collect(pairs)) do pair
        first(pair) isa Source
    end
    unspecified_sources = setdiff(sources_, first.(specified_source_pairs))
    unspecified_sources_wrapping_something =
        filter(s -> !isempty(s), unspecified_sources)
    if !isempty(unspecified_sources_wrapping_something) &&
        !empty_unspecified_sources
        @warn "No replacement specified for one or more non-empty source "*
        "nodes. Contents will be duplicated. "
    end
    if empty_unspecified_sources
        unspecified_source_pairs = [s => source() for
                                    s in unspecified_sources]
    else
        unspecified_source_pairs = [s => deepcopy(s) for
                                    s in unspecified_sources]
    end

    all_source_pairs = vcat(specified_source_pairs, unspecified_source_pairs)

    # drop source nodes from all nodes of network terminating at W:
    nodes_ = filter(nodes(W)) do N
        !(N isa Source)
    end
    isempty(nodes_) && error("All nodes in network are source nodes. ")
    # instantiate node and machine dictionaries:
    newnode_given_old =
        IdDict{AbstractNode,AbstractNode}(all_source_pairs)
    newsources = [newnode_given_old[s] for s in sources_]
    newinterface_node_given_old =
        IdDict{AbstractNode,AbstractNode}()
    newmach_given_old = IdDict{Machine,Machine}()

    # build the new network:
    for N in nodes_
       args = [newnode_given_old[arg] for arg in N.args]
         if N.machine === nothing
             newnode_given_old[N] = node(N.operation, args...)
         else
             if N.machine in keys(newmach_given_old)
                 m = newmach_given_old[N.machine]
             else
                 train_args = [newnode_given_old[arg] for arg in N.machine.args]
                 m = Machine(newmodel_given_old[N.machine.model],
                                train_args...)
                 newmach_given_old[N.machine] = m
             end
             newnode_given_old[N] = N.operation(m, args...)
         end
        if N in interface_nodes
            newinterface_node_given_old[N] = newnode_given_old[N]
        end
    end

    newinterface_nodes = Tuple(newinterface_node_given_old[N] for N in
                          interface_nodes)
    newsignature = NamedTuple{keys(signature)}(newinterface_nodes)

    return machine(mach.model, newsources...; newsignature...)

end


text:
using Main.CaesarM: caesarshift
using Main.ComparisonM: faComparison
using Main.FAM: fa

while true

    println("\nText to shift [leave blank to break]: ")
    human_input = readline()
    human_input == "" ? break :
    println("\nAlphabet [leave blank for default]: ")
    human_alphabet = readline()
    human_alphabet == "" ? human_alphabet = "abcdefghijklmnopqrstuvwxyz" :
    println()

    shifts = caesarshift.(uppercase(human_input), uppercase(human_alphabet), 1:length(human_alphabet))

    best = ""
    bestscore = -10000
    for shift in shifts
        score = faComparison(fa(uppercase(shift), uppercase(human_alphabet)))
        if score > bestscore
            best = shift
            bestscore = score
        end
    end

    println("We think the decryption is this (Assuming the text is in English): \n\n", best, "\n\n[Press enter to see them all]")
    readline()

    for item in range(1, length = length(shifts))
        println("\nShift of ", item % length(shifts), ": ", shifts[item])
    end

end

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
#EasyPlot: docstring definitions
#-------------------------------------------------------------------------------

#==
===============================================================================#

@doc """`fill(fn::Function d::DataRS, sweeplist)`

Construct a DataRS structure storing results from parametric sweeps using recursive data structures.

# Examples
```julia-repl
signal = fill(DataRS, PSweep("A", [1, 2, 4] .* 1e-3)) do A
    fill(DataRS, PSweep("phi", [0, 0.5, 1] .* (π/4))) do 𝜙
    fill(DataRS{DataF1}, PSweep("freq", [1, 4, 16] .* 1e3)) do 𝑓
       𝜔 = 2π*𝑓; T = 1/𝑓
       Δt = T/100 #Define resolution from # of samples per period
       Tsim = 4T #Simulated time
       t = DataF1(0:Δt:Tsim) #DataF1 creates a t:{y, x} container with y == x
       sig = A * sin(𝜔*t + 𝜙) #Still a DataF1 sig:{y, x=t} container
       return sig
end; end; end
```

Note that inner-most sweep needs to specify element type (DataF1).
Other (scalar) element types include: DataInt/DataFloat/DataComplex.
""" Base.fill(::Function, ::Type{DataRS}, args...)

@doc """`fill(value, [DataRS/DataF1], sweeplist)`

Construct a filled DataRS or DataHR structure with the provided list of parametric sweeps.

# Examples
```julia-repl
sweeplist = PSweep[
	PSweep("A", [1, 2, 4])
	PSweep("freq", [1, 2, 4, 8, 16] .* 1e3)
	PSweep("phi", π .* [1/4, 1/3, 1/2])
]

all32s = fill(32.0, DataRS, sweeplist)
```
""" Base.fill(::DF1_Num, ::Type{DataMD}, args...)

@doc """`zeros([DataRS/DataF1], sweeplist)`

Construct a 0.0-filled DataRS or DataHR structure with the provided list of parametric sweeps.

# Examples
```julia-repl
sweeplist = PSweep[
	PSweep("A", [1, 2, 4])
	PSweep("freq", [1, 2, 4, 8, 16] .* 1e3)
	PSweep("phi", π .* [1/4, 1/3, 1/2])
]

all0s = zeros(DataRS, sweeplist)
```
""" Base.zeros(::Type{DataMD}, args...)

@doc """`ones([DataRS/DataF1], sweeplist)`

Construct a 1.0-filled DataRS or DataHR structure with the provided list of parametric sweeps.

# Examples
```julia-repl
sweeplist = PSweep[
	PSweep("A", [1, 2, 4])
	PSweep("freq", [1, 2, 4, 8, 16] .* 1e3)
	PSweep("phi", π .* [1/4, 1/3, 1/2])
]

all1s = ones(DataRS, sweeplist)
```
""" Base.ones(::Type{DataMD}, args...)

#Last line


text:
function triple(::Type{T}, x::BigFloat) where {T<:IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    hi = T(x)
    md = T(x - hi)
    lo = T(x - hi - md)
    setprecision(BigFloat, prec)
    return hi, md, lo
end

function triple(::Type{T}, x::String) where {T<:IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    z = parse(BigFloat, x)
    hi = T(z)
    md = T(z - hi)
    lo = T(z - hi - md)
    setprecision(BigFloat, prec)
    return hi,md,lo
end

triple64(x::BigFloat) = triple(Float64, x)
triple32(x::BigFloat) = triple(Float32, x)
triple16(x::BigFloat) = triple(Float16, x)

triple64(x::String) = triple(Float64, x)
triple32(x::String) = triple(Float32, x)
triple16(x::String) = triple(Float16, x)

@inline function clean0s(hi::T, md::T, lo::T) where {T}
    if !iszero(hi)
        hi, md, lo
    elseif iszero(md)
        lo, md, hi
    else
        md, lo, hi
    end
end

@inline function clean0s(hi::T, lo::T) where {T}
    if !iszero(hi)
        hi, lo
    else
        lo, hi
    end
end

function triple_inv(::Type{T}, x::BigFloat) where {T<:IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    x = inv(x)
    hi = T(x)
    md = T(x - hi)
    lo = T(x - hi - md)
    setprecision(BigFloat, prec)
    return hi,md,lo
end

triple64inv(x::BigFloat) = triple_inv(Float64, x)
triple32inv(x::BigFloat) = triple_inv(Float32, x)
triple16inv(x::BigFloat) = triple_inv(Float16, x)

#=
Basic building blocks for a triple-double intermediate format
Christoph Quirin Lauter
Thème NUM — Systèmes numériques, Projet Arénaire
Rapport de recherche n 5702 — Septembre 2005 — 67 pages
=#

# Algorithm 3.3 [Triple-Double Building Blocks]

function renorm(hi::T, md::T, lo::T) where {T<:AbstractFloat}
    md, lo = two_sum(md, lo)
    hi, m  = two_sum(hi, md)
    md, lo = two_sum(m,  md)
    return hi, md, lo
end

function renorm_hilo(hi::T, md::T, lo::T) where {T<:AbstractFloat}
    md, lo = two_hilo_sum(md, lo)
    hi, m  = two_hilo_sum(hi, md)
    md, lo = two_hilo_sum(m,  md)
    return hi, md, lo
end

# Algorithm 4.1  (relerr 4.5 .. 16 u^2) [Triple-Double Building Blocks]

function add222(ahilo, bhilo)
    ahi, alo = ahilo
    bhi, blo = bhilo
    t1 = ahi + bhi
    if abs(ahi) >= abs(bhi)
        t2 = ahi - t1
        t3 = t2  + bhi
        t4 = t3  + blo
        t5 = t4  + alo
    else
        t2 = bhi - t1
        t3 = t2  + ahi
        t4 = t3  + alo
        t5 = t4  + blo
    end
    zhi, zlo = two_sum(t1, t5)
    zhi, zlo = clean0s(zhi, zlo)
    return zhi, zlo
end

"""
    _Basic building blocks for a triple-double intermediate format_

Christoph Quirin Lauter, 2005
""" TripleDouble_BuildingBlocks

# Algorithm 4.6  (relerr 16 u^2) [Triple-Double Building Blocks]
"""
    mul222( (ahi, alo), (bhi, blo) ) ↦ (chi, clo)

### Preconditions
- abs(alo) ≤ abs(ahi) * pow(2, -53)
- abs(blo) ≤ abs(bhi) * pow(2, -53)

### Postconditions
- abs(clo) ≤ abs(ahi) * pow(2, -53)
- chi + clo == ((ahi + alo) * (bhi + blo)) * (1 + ε), abs(ε) ≤ pow(2, -102)
- chi + clo == chi

#### [TripleDouble_BuildingBlocks](@ref)
"""
function mul222(ahilo, bhilo)
    ahi, alo = ahilo
    bhi, blo = bhilo
    t1, t2 = two_prod(ahi, bhi)
    t3 = ahi * blo
    t4 = alo * bhi
    t5 = t3 + t4
    t6 = t2 + t5
    zhi, zlo = two_sum(t1, t6)
    return zhi, zlo
end

# addition

# Algorithm 5.1 (Add33) [Triple-Double Building Blocks]
# (ahi,amd,alo) + (bhi,bmd,blo) :: (zhi,zmd,zlo)

function add333(ahi::T, amd::T, alo::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, bhi)
    t2, t3 = two_sum(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t6 = alo + blo
    t5 = t3 + t4
    t8 = t5 + t6
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function add333(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return add333(a[1], a[2], a[3], b[1], b[2], b[3])
end

function add332(ahi::T, amd::T, alo::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, bhi)
    t2, t3 = two_sum(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t6 = alo + blo
    t5 = t3 + t4
    zmd = t5 + t6
    zmd += t7
    zhi, zmd = clean0s(zhi, zmd)
    return zhi, zmd
end

@inline function add332(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return add332(a[1], a[2], a[3], b[1], b[2], b[3])
end

function add323(ahi::T, amd::T, alo::T, bhi::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, bhi)
    t2, t3 = two_sum(amd, blo)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    t8 = t5 + alo
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function add323(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return add323(a[1], a[2], a[3], b[1], b[2])
end

function add322(ahi::T, amd::T, alo::T, bhi::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, bhi)
    t2, t3 = two_sum(amd, blo)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    zmd = t5 + alo
    zmd += t7
    zhi, zmd = clean0s(zhi, zmd)
    return zhi, zmd
end

@inline function add322(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return add322(a[1], a[2], a[3], b[1], b[2])
end

function add223(ahi::T, amd::T, bhi::T, bmd::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, bhi)
    t2, t3 = two_sum(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    zmd, zlo = two_sum(t7, t5)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function add223(a::Tuple{T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return add223(a[1], a[2], b[1], b[2])
end

function add313(ahi::T, amd::T, alo::T, b::T) where {T<:AbstractFloat}
    zhi, t1 = two_sum(ahi, b)
    t7, t4 = two_sum(t1, amd)
    t8 = t4 + alo
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

add133(b::T, ahi::T, amd::T, alo::T) where {T<:AbstractFloat} = add313(ahi, amd, alo, b)

add313(a::NTuple{3,T}, b::NTuple{1,T}) where {T<:AbstractFloat} =
    add313(a[1], a[2], a[3], b[1])

add133(b::NTuple{1,T}, a::NTuple{3,T}) where {T<:AbstractFloat} =
    add313(a[1], a[2], a[3], b[1])

# subtraction

function sub333(ahi::T, amd::T, alo::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t6 = alo - blo
    t5 = t3 + t4
    t8 = t5 + t6
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function sub333(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return sub333(a[1], a[2], a[3], b[1], b[2], b[3])
end

function sub332(ahi::T, amd::T, alo::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t6 = alo - blo
    t5 = t3 + t4
    t8 = t5 + t6
    zmd = t7 + t8
    zhi, zmd = clean0s(zhi, zmd)
    return zhi, zmd
end

@inline function sub332(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return sub332(a[1], a[2], a[3], b[1], b[2], b[3])
end

function sub323(ahi::T, amd::T, alo::T, bhi::T, bmd::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    t8 = t5 + alo
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function sub323(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return sub323(a[1], a[2], a[3], b[1], b[2])
end

function sub322(ahi::T, amd::T, alo::T, bhi::T, bmd::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    zmd = t5 + alo
    zmd += t7
    zhi, zmd = clean0s(zhi, zmd)
    return zhi, zmd
end

@inline function sub322(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return sub322(a[1], a[2], a[3], b[1], b[2])
end

function sub233(ahi::T, amd::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    t8 = t5 - blo
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function sub233(a::Tuple{T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return sub233(a[1], a[2], b[1], b[2], b[3])
end

function sub232(ahi::T, amd::T, bhi::T, bmd::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    zmd = t5 - blo
    zmd += t7
    zhi, zmd = clean0s(zhi, zmd)
    return zhi, zmd
end

@inline function sub232(a::Tuple{T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    return sub232(a[1], a[2], b[1], b[2], b[3])
end

function sub223(ahi::T, amd::T, bhi::T, bmd::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, bhi)
    t2, t3 = two_diff(amd, bmd)
    t7, t4 = two_sum(t1, t2)
    t5 = t3 + t4
    zmd, zlo = two_sum(t7, t5)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

@inline function sub223(a::Tuple{T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return sub223(a[1], a[2], b[1], b[2])
end

function sub313(ahi::T, amd::T, alo::T, b::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, b)
    t7, t4 = two_sum(t1, amd)
    t8 = t4 + alo
    zmd, zlo = two_sum(t7, t8)
    zhi, zmd, zlo = clean0s(zhi,zmd,zlo)
    return zhi, zmd, zlo
end

function sub312(ahi::T, amd::T, alo::T, b::T) where {T<:AbstractFloat}
    zhi, t1 = two_diff(ahi, b)
    t7, t4 = two_sum(t1, amd)
    zlo = t4 + alo
    zlo += t7
    zhi, zlo = clean0s(zhi, zlo)
    return zhi, zlo
end

sub313(a::NTuple{3,T}, b::NTuple{1,T}) where {T<:AbstractFloat} =
    sub313(a[1], a[2], a[3], b)

sub312(a::NTuple{3,T}, b::NTuple{1,T}) where {T<:AbstractFloat} =
    sub312(a[1], a[2], a[3], b)

# Algorithm 6.1 [Triple-Double Building Blocks]
# (ahi,alo) * (bhi,blo) :: (zhi,zmd,zlo)

function mul223(ahi::T, alo::T, bhi::T, blo::T) where {T<:AbstractFloat}
    zhi, t1 = two_prod(ahi, bhi)
    t2, t3 = two_prod(ahi, blo)
    t4, t5 = two_prod(alo, bhi)
    t6 = alo * blo
    t7, t8 = two_sum(t2, t3, t4, t5)
    t9, t10 = two_sum(t1, t6)
    zmd, zlo = two_sum(t7, t8, t9, t10)
    return zhi, zmd, zlo
end

@inline function mul223(a::Tuple{T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
    return mul223(a[1], a[2], b[1], b[2])
end

function mul323(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
  ahi, amd, alo = a
  bhi, blo = b
  p0,q0 = two_prod(ahi, bhi)
  p1,q1 = two_prod(ahi, blo)
  p2,q2 = two_prod(amd, bhi)
  p4,q4 = two_prod(amd, blo)
  p5,q5 = two_prod(alo, bhi)

  # Start Accumulation
  p1,p2,q0 = three_sum(p1, p2, q0)

  # Six-Three Sum  of p2, q1, q2, p3, p4, p5
  p2,q1,q2 = three_sum(p2, q1, q2)
  p3,p4 = two_sum(p4, p5)
  # compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5)
  s0,t0 = two_sum(p2, p3)
  s1,t1 = two_sum(q1, p4)
  s2 = q2
  s1,t0 = two_sum(s1, t0)
  s2 += (t0 + t1)

  # O(eps^3) order terms
  s1 += alo*blo + q0 + q4 + q5
  #p0,p1,s0 = renormAs3(p0, p1, s0, s1+s2)
  s1 += s2
  s0,s1 = two_hilo_sum(s0,s1)
  p1,s0 = two_hilo_sum(p1,s0)
  p0,p1 = two_hilo_sum(p0,p1)

  return p0,p1,s0
end

mul323(a::Tuple{Float64, Float64, Float64}, b::Tuple{Float32, Float32}) =
    mul323(a, (Float64(b[1]), Float64(b[2])) )
mul323(a::Tuple{Float32, Float32, Float32}, b::Tuple{Float64, Float64}) =
    mul323((Float64(a[1]), Float64(a[2]), Float64(a[3])), b)

function mul322(a::Tuple{T,T,T}, b::Tuple{T,T}) where {T<:AbstractFloat}
  ahi, amd, alo = a
  bhi, blo = b
  p0,q0 = two_prod(ahi, bhi)
  p1,q1 = two_prod(ahi, blo)
  p2,q2 = two_prod(amd, bhi)
  p4,q4 = two_prod(amd, blo)
  p5,q5 = two_prod(alo, bhi)

  # Start Accumulation
  p1,p2,q0 = three_sum(p1, p2, q0)

  # Six-Three Sum  of p2, q1, q2, p3, p4, p5
  p2,q1,q2 = three_sum(p2, q1, q2)
  p3,p4 = two_sum(p4, p5)
  # compute (s0, s1, s2) = (p2, q1, q2) + (p3, p4, p5)
  s0,t0 = two_sum(p2, p3)
  s1,t1 = two_sum(q1, p4)
  s2 = q2
  s1,t0 = two_sum(s1, t0)
  s2 += (t0 + t1)

  # O(eps^3) order terms
  s1 += alo*blo + q0 + q4 + q5
  #p0,p1,s0 = renormAs3(p0, p1, s0, s1+s2)
  s1 += s2
  s0 += s1
  p1 += s0
  p0,p1 = two_hilo_sum(p0,p1)

  return p0,p1
end

mul322(a::Tuple{Float64, Float64, Float64}, b::Tuple{Float32, Float32}) =
    mul322(a, (Float64(b[1]), Float64(b[2])) )
mul322(a::Tuple{Float32, Float32, Float32}, b::Tuple{Float64, Float64}) =
    mul322((Float64(a[1]), Float64(a[2]), Float64(a[3])), b)

mul233(a::Tuple{T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat} =
    mul323(b,a)

mul232(a::Tuple{T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat} =
    mul322(b,a)

function mul333(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    ahi, amd, alo = a
    bhi, bmd, blo = b

    hi,t1 = two_prod(ahi, bhi)
    t2,t3 = two_prod(ahi, bmd)
    t4,t5 = two_prod(amd, bhi)
    t6,t7 = two_prod(amd, bmd)

    t8  = ahi * blo
    t9  = alo * bhi
    t10 = amd * blo
    t11 = alo * bmd
    t12 = t8  + t9
    t13 = t10 + t11

    t14, t15 = two_hilo_sum(t1, t6)

    t16 = t7  + t15
    t17 = t12 + t13
    t18 = t16 + t17

    t19, t20 = two_hilo_sum(t14, t18)
    t21, t22 = two_hilo_sum(t2, t3, t4, t5)
    md,  lo  = two_hilo_sum(t21, t22, t19, t20)

    return hi, md, lo
end

mul333(a::Tuple{Float64, Float64, Float64}, b::Tuple{Float32, Float32, Float32}) =
    mul333(a, (Float64(b[1]), Float64(b[2]), Float64(b[3])) )
mul333(a::Tuple{Float32, Float32, Float32}, b::Tuple{Float64, Float64, Float64}) =
    mul333((Float64(a[1]), Float64(a[2]), Float64(a[3])), b)

function mul332(a::Tuple{T,T,T}, b::Tuple{T,T,T}) where {T<:AbstractFloat}
    ahi, amd, alo = a
    bhi, bmd, blo = b

    hi,t1 = two_prod(ahi, bhi)
    t2,t3 = two_prod(ahi, bmd)
    t4,t5 = two_prod(amd, bhi)
    t6,t7 = two_prod(amd, bmd)

    t8  = ahi * blo
    t9  = alo * bhi
    t10 = amd * blo
    t11 = alo * bmd

    t14, lo = two_hilo_sum(t1, t6)

    lo += t7
    lo += t8  + t9
    lo += t10 + t11
    lo += t16 + t17
    lo += t14
    lo += t5  + t4
    lo += t3  + t2

    return hi, lo
end

mul332(a::Tuple{Float64, Float64, Float64}, b::Tuple{Float32, Float32, Float32}) =
    mul332(a, (Float64(b[1]), Float64(b[2]), Float64(b[3])) )
mul332(a::Tuple{Float32, Float32, Float32}, b::Tuple{Float64, Float64, Float64}) =
    mul332((Float64(a[1]), Float64(a[2]), Float64(a[3])), b)

function mul313(a::Tuple{T,T,T}, b::Tuple{T}) where {T<:AbstractFloat}
  ahi, amd, alo = a
  bhi = b[1]
  p0,q0 = two_prod(ahi, bhi)
  p2,q2 = two_prod(amd, bhi)
  p5,q5 = two_prod(alo, bhi)

  # Start Accumulation
  p1,p2 = two_sum(p2, q0)

  p2,q1 = two_sum(p2, q2)

  s0,t0 = two_sum(p2, p5)
  s1,t0 = two_sum(q1, t0)

  # O(eps^3) order terms
  s1 += q5
  #p0,p1,s0 = renormAs3(p0, p1, s0, s1+s2)
  s1 += t0
  s0,s1 = two_hilo_sum(s0,s1)
  p1,s0 = two_hilo_sum(p1,s0)
  p0,p1 = two_hilo_sum(p0,p1)

  return p0,p1,s0
end

mul133(a::NTuple{1,T}, b::NTuple{3,T}) where {T<:AbstractFloat} = mul313(b, a)

function mul312(a::Tuple{T,T,T}, b::Tuple{T}) where {T<:AbstractFloat}
  ahi, amd, alo = a
  bhi = b[1]
  p0,q0 = two_prod(ahi, bhi)
  p2,q2 = two_prod(amd, bhi)
  p5,q5 = two_prod(alo, bhi)

  # Start Accumulation
  p1,p2 = two_sum(p2, q0)

  p2,q1 = two_sum(p2, q2)

  s0,t0 = two_sum(p2, p5)
  s1,t0 = two_sum(q1, t0)

  # O(eps^3) order terms
  s1 += q5
  #p0,p1,s0 = renormAs3(p0, p1, s0, s1+s2)
  s1 += t0
  p1 = p1 + (s0 + s1)
  p0,p1 = two_hilo_sum(p0,p1)

  return p0,p1
end

mul132(a::NTuple{1,T}, b::NTuple{3,T}) where {T<:AbstractFloat} = mul312(b, a)

function muladd222(a::Tuple{T,T}, b::Tuple{T,T}, c::Tuple{T,T}) where {T<:AbstractFloat}
    m = mul223(a, b)
    hi, lo = add322(m, c)
    return hi, lo
end


text:
in_path = ARGS[1]

out_path = try
                ARGS[2]
            catch
                replace(in_path, ".csv" => "_formatted.csv")
            end

file = open(in_path, "r")
file_content = read(file, String)
file_check = eachmatch(r"(?<day>\d{1,2})/(?<month>\d{1,2})/(?<year>\d{4}) (?<hour>\d{1,2}):(?<minute>\d{1,2}),(?<wl19>\d{0,9}.?\d{0,9}),(?<wl2>\d{0,9}.?\d{0,9}),(?<wl21>\d{0,9}.?\d{0,9})", file_content)
timesteps = [replace(x[:year]*"-"*x[:month]*"-"*x[:day]*" "*x[:hour]*":"*x[:minute], r"\D(\d):(\d{2})" => s" 0\1:\2") for x in file_check]

data = Dict()
data["19"] = [each[:wl19] for each in file_check]
data["2"] = [each[:wl2] for each in file_check]
data["21"] = [replace(each[:wl21], r"\r" => "") for each in file_check]

names = Dict()
names["19"] = "Lagoon"
names["2"] = "Cabbage_in"
names["21"] = "Cabbage"

file_out = open(out_path, "w")
write(file_out, "name,id,"*join(timesteps, ",")*"\n")
for level_id in ["19", "2", "21"]
    write(file_out, names[level_id]*","*level_id*","*join(data[level_id], ",")*"\n")
end

max_stars_repo_name: smart-facility/smart-faware

max_issues_repo_name: smart-facility/smart-faware

max_forks_repo_name: smart-facility/smart-faware


text:
module TestNItems

using DataTools
using Test
using Transducers
using Transducers: IdentityTransducer

@testset "_pop_innermost_maplikes" begin
    pop(args...) = DataTools._pop_innermost_maplikes(opcompose(args...))
    @test pop(Map(inv)) === IdentityTransducer()
    @test pop(MapSplat(tuple), Map(inv)) === IdentityTransducer()
    @test pop(Filter(isodd), MapSplat(tuple), Map(inv)) === Filter(isodd)
    @test pop(Map(isodd), Filter(isodd), MapSplat(tuple), Map(inv)) ===
          opcompose(Map(isodd), Filter(isodd))
end

@testset "nitems" begin
    @test nitems(1:10) == 10
    @test nitems(error(x) for x in 1:10) == 10
    @test 1:10 |> Map(error) |> MapSplat(error) |> Scan(+) |> nitems == 10
    @test 1:10 |> Filter(isodd) |> Map(error) |> MapSplat(error) |> nitems == 5
    @test 1:10 |> Filter(isodd) |> Map(x -> x ÷ 3) |> Filter(isodd) |> nitems == 3
end

end  # module


text:
# convert padding etc. size to an Int array of the right dimension
function psize(p, x)
  nd = ndims(x)-2
  if isa(p,Number)
    fill(Int(p),nd)
  elseif length(p)==nd
    collect(Int,p)
  else
    throw(DimensionMismatch("psize: $p $nd"))
  end
end

function im2col_2d!{T}(img::AbstractArray{T,3}, col::AbstractArray{T,2}, width::Int, height::Int, channels::Int,
  kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int, stride_h::Int, mode::Int)

  height_col = div((height + 2pad_h - kernel_h), stride_h) + 1
  width_col = div((width + 2pad_w - kernel_w), stride_w) + 1
  channels_col = channels * kernel_h * kernel_w

  #pragma omp parallel for
  for c = 1:channels_col
    w_offset = (c - 1) % kernel_w
    h_offset = div(c - 1, kernel_w) % kernel_h
    c_im = div(c - 1, kernel_h * kernel_w)
    if mode == 0
      w_offset = kernel_w - 1 - w_offset
      h_offset = kernel_h - 1 - h_offset
    end
    for h = 1:height_col
      for w = 1:width_col
        h_pad = (h - 1) * stride_h - pad_h + h_offset
        w_pad = (w - 1) * stride_w - pad_w + w_offset
        if h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width
          col[((c - 1)*height_col+h-1) * width_col + w] =
           img[(c_im  * height + h_pad) * width + w_pad + 1]
        else
          col[((c - 1)*height_col+h - 1) * width_col + w] = 0
        end
      end
    end
  end
end

function col2im_2d!{T}(col::AbstractArray{T,2}, img::AbstractArray{T,3}, width::Int, height::Int,
  channels::Int, kernel_w::Int, kernel_h::Int, pad_w::Int, pad_h::Int, stride_w::Int,
  stride_h::Int, mode::Int)

  height_col = div(height + 2pad_h - kernel_h, stride_h) + 1
  width_col = div(width + 2pad_w - kernel_w, stride_w) + 1
  channels_col = channels * kernel_h * kernel_w

  fill!(img, 0)
  #pragma omp parallel for
  for c = 1:channels_col
    w_offset = (c - 1) % kernel_w
    h_offset = div(c - 1,  kernel_w) % kernel_h
    c_im = div(c - 1, kernel_h * kernel_w)
    if mode == 0
      w_offset = kernel_w - 1 - w_offset
      h_offset = kernel_h - 1 - h_offset
    end
    for h = 1:height_col, w = 1:width_col
      h_pad = (h - 1) * stride_h - pad_h + h_offset
      w_pad = (w - 1) * stride_w - pad_w + w_offset
      if h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width
        cval::T = col[((c - 1) * height_col + h - 1) * width_col + w]
        img[(c_im * height + h_pad) * width + w_pad + 1] += cval
      end
    end
  end
end

function im2col_3d!{T}(img::AbstractArray{T,4}, col::AbstractArray{T,2}, width::Int, height::Int, depth::Int,
  channels::Int, kernel_w::Int, kernel_h::Int, kernel_d::Int, pad_w::Int, pad_h::Int, pad_d::Int,
  stride_w::Int, stride_h::Int, stride_d::Int, mode::Int)

  height_col = div((height + 2pad_h - kernel_h), stride_h) + 1
  width_col = div((width + 2pad_w - kernel_w), stride_w) + 1
  depth_col = div((depth + 2pad_d - kernel_d), stride_d) + 1
  channels_col = channels * kernel_h * kernel_w * kernel_d

  #pragma omp parallel for
  for c = 1:channels_col
    w_offset = (c - 1) % kernel_w
    h_offset = div(c - 1, kernel_w) % kernel_h
    d_offset = div(c - 1, kernel_w * kernel_h) % kernel_d
    c_im = div(c - 1, kernel_w * kernel_h * kernel_d)
    if mode == 0
      w_offset = kernel_w - 1 - w_offset
      h_offset = kernel_h - 1 - h_offset
      d_offset = kernel_d - 1 - d_offset
    end
    for d = 1:depth_col, h = 1:height_col, w = 1:width_col
      d_pad = (d - 1) * stride_d - pad_d + d_offset
      h_pad = (h - 1) * stride_h - pad_h + h_offset
      w_pad = (w - 1) * stride_w - pad_w + w_offset
      if d_pad >= 0 && d_pad < depth && h_pad >= 0 && h_pad < height &&
        w_pad >= 0 && w_pad < width
        col[(((c - 1) * depth_col + d - 1) * height_col + h - 1) * width_col + w] =
    	    img[((c_im * depth + d_pad) * height + h_pad) * width + w_pad + 1]
    	else
    	  col[(((c - 1) * depth_col + d - 1) * height_col + h - 1) * width_col + w] = 0
      end
    end
  end
end

function col2im_3d!{T}(col::AbstractArray{T,2}, img::AbstractArray{T,4}, width::Int, height::Int,
  depth::Int, channels::Int, kernel_w::Int, kernel_h::Int, kernel_d::Int,
  pad_w::Int, pad_h::Int, pad_d::Int, stride_w::Int, stride_h::Int, stride_d::Int, mode::Int)

  depth_col = div(depth + 2pad_d - kernel_d, stride_d) + 1
  height_col = div(height + 2pad_h - kernel_h, stride_h) + 1
  width_col = div(width + 2 * pad_w - kernel_w, stride_w) + 1
  channels_col = channels * kernel_h * kernel_w * kernel_d

  fill!(img, 0)
  #pragma omp parallel for
  for c = 1:channels_col
    w_offset = (c - 1) % kernel_w;
    h_offset = div(c - 1, kernel_w) % kernel_h
    d_offset = div(c - 1, kernel_w * kernel_h) % kernel_d
    c_im = div(c - 1, kernel_h * kernel_w * kernel_d)

    if mode == 0
      w_offset = kernel_w - 1 - w_offset
      h_offset = kernel_h - 1 - h_offset
      d_offset = kernel_d - 1 - d_offset
    end

    for d = 1:depth_col, h = 1:height_col, w = 1:width_col
      d_pad = (d - 1) * stride_d - pad_d + d_offset
    	h_pad = (h - 1) * stride_h - pad_h + h_offset
    	w_pad = (w - 1) * stride_w - pad_w + w_offset
    	if h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width &&
        d_pad >= 0 && d_pad < depth
    	  cval::T = col[(((c - 1) * depth_col + d - 1) * height_col + h - 1) * width_col + w]
    	  iidx = ((c_im * depth + d_pad) * height + h_pad) * width + w_pad + 1
              #pragma omp atomic
    	  img[iidx] += cval
    	end
    end
  end
end

function im2col_dims(w,y)
    N = ndims(y)
    r,c = 1,1
    for i=1:N-2
        r *= size(y,i)
        c *= size(w,i)
    end
    c *= size(w,N-1)
    return (r, c)
end

function conv2d!{T}(y::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
                  padding=0, stride=1, mode=0, alpha=T(1))
    if mode != 0 && mode != 1; throw(ArgumentError("conv2d only supports mode=0 or 1.")); end
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    if Cx!=C1; throw(DimensionMismatch()); end
    Wy,Hy,Cy,Ny = size(y)
    x2dims = im2col_dims(w,y)
    x2 = similar(x, x2dims)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    M,N,K,Y = Wy*Hy,Cy,Ww*Hw*Cx,Wy*Hy*Cy
    yidx = 1
    @inbounds for n in 1:Nx
        im2col2d!(w, x, x2, n, p1, p2, s1, s2, mode)
        gemm!('N','N',M,N,K,alpha,pointer(x2),pointer(w),T(0),pointer(y,yidx))
        yidx += Y
    end
    return y
end

function conv2d_grad_w!{T}(dw::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4}, dy::AbstractArray{T,4};
                   padding=0, stride=1, mode=0, alpha=1)
    # dw = x'*dy
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    Wy,Hy,Cy,Ny = size(dy)
    # if mode != 0 && mode != 1; throw(ArgumentError("conv2d only supports mode=0 or 1.")); end
    # @assert Cx==C1 && Cy==C2 && Ny==Nx
    x2dims = im2col_dims(w,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Cy,Ww*Hw*Cx,Cy,Wy*Hy
    alpha,beta = T(alpha),T(1)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        im2col2d!(w, x, x2, n, p1, p2, s1, s2, mode)
        gemm!('T','N',M,N,K,alpha,pointer(x2),pointer(dy,dyi),beta,pointer(dw))
        dyi += Y
    end
    return dw
end

function conv2d_grad_x!{T}(dx::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4}, dy::AbstractArray{T,4};
                   padding=0, stride=1, mode=0, alpha=1)
    # dx = dy*w'
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    Wy,Hy,Cy,Ny = size(dy)
    # if mode != 0 && mode != 1; throw(ArgumentError("conv2d only supports mode=0 or 1.")); end
    @assert Cx==C1 && Cy==C2 && Ny==Nx
    x2dims = im2col_dims(w,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Cy,Wy*Hy,Ww*Hw*Cx,Cy
    alpha,beta = T(alpha),T(0)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        gemm!('N','T',M,N,K,alpha,pointer(dy,dyi),pointer(w),beta,pointer(x2))
        col2im2d!(w,dx,x2,n,p1,p2,s1,s2,mode)
        dyi += Y
    end
    return dx
end

function im2col2d!(w::AbstractArray{T,4}, x::AbstractArray{T,4}, x2::AbstractArray{T,2},
                 n::Int, p1::Int, p2::Int, s1::Int, s2::Int, mode::Int) where T
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    xn = x[:, :, :, n]
    im2col_2d!(xn,x2,Wx,Hx,Cx,Ww,Hw,p1,p2,s1,s2,mode)
    return x2
end

function col2im2d!(w::AbstractArray{T,4}, x::AbstractArray{T,4}, x2::AbstractArray{T,2},
                 n::Int, p1::Int, p2::Int, s1::Int, s2::Int, mode::Int) where T
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    xn = x[:, :, :, n]
    col2im_2d!(x2,xn,Wx,Hx,Cx,Ww,Hw,p1,p2,s1,s2,mode)
    x[:, :, :, n] = xn
    return x
end

function conv3d!{T}(y::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5};
                  padding=0, stride=1, mode=0, alpha=T(1))
    if mode != 0 && mode != 1; throw(ArgumentError("conv3d only supports mode=0 or 1.")); end
    Wx,Hx,Dx,Cx,Nx = size(x)
    Ww,Hw,Dw,C1,C2 = size(w)
    if Cx!=C1; throw(DimensionMismatch()); end
    Wy,Hy,Dy,Cy,Ny = size(y)
    # @assert Cy==C2 && Ny==Nx
    x2dims = im2col_dims(w,y)
    x2 = similar(x, x2dims)
    (p1,p2,p3) = psize(padding,x)
    (s1,s2,s3) = psize(stride,x)
    M,N,K,Y = Wy*Hy*Dy,Cy,Ww*Hw*Dw*Cx,Wy*Hy*Dy*Cy
    yidx = 1
    W = reshape(w, (Ww,:,C1,C2))
    @inbounds for n in 1:Nx
        im2col3d!(w, x, x2, n, p1, p2, p3, s1, s2, s3, mode)
        gemm!('N','N',M,N,K,alpha,pointer(x2),pointer(W),T(0),pointer(y,yidx))
        yidx += Y
    end
    return y
end

function conv3d_grad_w!{T}(dw::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5}, dy::AbstractArray{T,5};
                   padding=0, stride=1, mode=0, alpha=1)
    # dw = x'*dy
    Wx,Hx,Dx,Cx,Nx = size(x)
    Ww,Hw,Dw,C1,C2 = size(w)
    Wy,Hy,Dy,Cy,Ny = size(dy)
    # if mode != 0 && mode != 1; throw(ArgumentError("conv2d only supports mode=0 or 1.")); end
    # @assert Cx==C1 && Cy==C2 && Ny==Nx
    x2dims = im2col_dims(w,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Dy*Cy,Ww*Hw*Dw*Cx,Cy,Wy*Hy*Dy
    alpha,beta = T(alpha),T(1)
    (p1,p2,p3) = psize(padding,x)
    (s1,s2,s3) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        im2col3d!(w, x, x2, n, p1, p2, p3, s1, s2, s3, mode)
        gemm!('T','N',M,N,K,alpha,pointer(x2),pointer(dy,dyi),beta,pointer(dw))
        dyi += Y
    end
    return dw
end

function conv3d_grad_x!{T}(dx::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5}, dy::AbstractArray{T,5};
                   padding=0, stride=1, mode=0, alpha=1)
    # dx = dy*w'
    Wx,Hx,Dx,Cx,Nx = size(x)
    Ww,Hw,Dw,C1,C2 = size(w)
    Wy,Hy,Dy,Cy,Ny = size(dy)
    # if mode != 0 && mode != 1; throw(ArgumentError("conv2d only supports mode=0 or 1.")); end
    @assert Cx==C1 && Cy==C2 && Ny==Nx
    x2dims = im2col_dims(w,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Dy*Cy,Wy*Hy*Dy,Ww*Hw*Dw*Cx,Cy
    alpha,beta = T(alpha),T(0)
    (p1,p2,p3) = psize(padding,x)
    (s1,s2,s3) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        gemm!('N','T',M,N,K,alpha,pointer(dy,dyi),pointer(w),beta,pointer(x2))
        col2im3d!(w,dx,x2,n,p1,p2,p3,s1,s2,s3,mode)
        dyi += Y
    end
    return dx
end

function im2col3d!(w::AbstractArray{T,5}, x::AbstractArray{T,5}, x2::AbstractArray{T,2},
                 n::Int, p1::Int, p2::Int, p3::Int, s1::Int, s2::Int,
                 s3::Int, mode::Int) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Ww,Hw,Dw,C1,C2 = size(w)
    xn = x[:, :, :, :, n]
    im2col_3d!(xn,x2,Wx,Hx,Dx,Cx,Ww,Hw,Dw,p1,p2,p3,s1,s2,s3,mode)
    return x2
end

function col2im3d!(w::AbstractArray{T,5}, x::AbstractArray{T,5}, x2::AbstractArray{T,2},
                 n::Int, p1::Int, p2::Int, p3::Int, s1::Int, s2::Int,
                 s3::Int, mode::Int) where T
    Wx,Hx,Dx,Cx,Nx = size(x)
    Ww,Hw,Dw,C1,C2 = size(w)
    xn = x[:, :, :, :, n]
    col2im_3d!(x2,xn,Wx,Hx,Dx,Cx,Ww,Hw,Dw,p1,p2,p3,s1,s2,s3,mode)
    x[:, :, :, :, n] = xn
    return x
end


text:
using ..Bodies

using ReferenceFrameRotations: angle_to_dcm, ddcm

for name in ALL_NAMES
    endswith(name, "Barycenter") && continue

    frame = Symbol("IAU", name)
    cname = Symbol("iau_", lowercase(name))
    @eval begin
        struct $frame <: RotatingFrame end
        const $cname = $frame()
        from_sym(::Val{$(Meta.quot(frame))}) = $cname
        export $frame, $cname

        function Rotation(::ICRF, ::$frame, ep::Epoch)
            tdb = TDBEpoch(ep)
            body = $(Symbol(name))()
            euler = euler_angles(body, tdb)
            euler′ = collect(euler_rates(body, tdb))

            m = angle_to_dcm(euler..., :ZXZ)
            m′ = ddcm(m, euler′)

            Rotation{icrf, $cname}(m, m′)
        end

        Rotation(::$frame, ::ICRF, ep::Epoch) = inv(Rotation(icrf, $cname, ep))
    end
end


text:
include("header.jl")

#=
# Dimensional analysis example - pump flow
In this example, we will find an appropriate set of dimensionless parameters that
relate pressure increase across a pump to the flow rate through the pump. We start
with the basic expectation that pressure increase $\Delta p$ should depend on $D$
(the diameter of the pump), $\rho$ (the fluid density), $\omega$ (the rotational speed),
and $Q$ (the volume flow rate)

$$\Delta p = f(D,\rho,\omega,Q)$$

=#

# ### Set up the module
using MAE103
#-
using Plots
using Serialization # to read/write data

#=
### Read in some raw data
First, we will read in some raw experimental data for a pump of diameter 12 in.
These data give the measured pressure rise (as *head increase*, $h_p$, units of meters)
vs volume flow rate (in units of cubic meters per second).

We will load in data for pump speeds of 250, 500 and 1000 rpm. For later convenience,
we will put all of these data into one data structure called `pumpdata`.
=#
pumpdata = deserialize("pumpdata.jls")

#=
The data are in the form of a `Dict` (short for dictionary). For example,
to get the flow rates for speed 500 rpm,
=#
Q500 = pumpdata["500rpm"]["Q"]

# and the associated head increases are
h500 = pumpdata["500rpm"]["head"]

#=
To convert the flow rate data to liters/minute (lpm), for example
=#
value.(VolumeFlowRate.(Q500),u"L/minute")
#=
or head increase to ft
=#
value.(Head.(h500),u"ft")

#=
Let's plot this data
=#
scatter(xlabel=L"$Q$",ylabel=L"$h$",legend=true,xlims=(0,Inf),ylims=(0,Inf))
scatter!(pumpdata["250rpm"]["Q"],pumpdata["250rpm"]["head"],label="250 rpm")
scatter!(pumpdata["500rpm"]["Q"],pumpdata["500rpm"]["head"],label="500 rpm")
scatter!(pumpdata["1000rpm"]["Q"],pumpdata["1000rpm"]["head"],label="1000 rpm")
#=
The data don't appear very similar to each other, except for a rough shape.
And more importantly, these are only curves for three pump speeds. What if
we want to predict the performance at, say, 700 rpm?
=#
#-
#=
### Find a suitable set of dimensionless parameters ($\Pi$ groups)
Here, we should think of $\Delta p$ as the dependent variable and $Q$ as the
independent variable -- this is what defines a pump's characteristics. The other
quantities are the reference quantities. It is straightforward to come up with the
$\Pi$ groups:

$$\dfrac{\Delta p}{\rho \omega^2 D^2} = f \left( \dfrac{Q}{\omega D^3} \right)$$
=#
#-
#=
### Rewrite the first $\Pi$ group in terms of the pump's head gain, $h_p$.
From the energy equation, $h_p + p_1/\rho g = p_2/\rho g$, where $\Delta p = p_2 - p_1$.
Thus, we can rewrite the first $\Pi$ group in terms of $h_p$, so that the relationship is

$$\dfrac{g h_p}{\omega^2 D^2} = f \left( \dfrac{Q}{\omega D^3} \right)$$

The dimensionless flow rate is called the *flow coefficient*, and the
dimensionless head is called the *head rise coefficient*.
=#

#=
### Calculate the $\Pi$ groups and plot them on the same plot
We will now use the raw dimensional data in `pumpdata` to calculate the
$\Pi$ groups. We set the basic parameters here:
=#
D = 12u"inch"
g = 1u"ge" # This is Earth's gravitation constant

#=
So, let's say we want to calculate the flow coefficient at 250 rpm:
=#
omega = 250u"rpm"
pumpdata["250rpm"]["Q"]/(omega*D^3)

#=
It looks like the data have units! But on closer inspection, they are
actually units that will cancel (with appropriate conversions). We can make that
happen by using the `DimensionlessParameter` function. (We use this in
vectorized form on each element in the array.) We'll assign the result to
a new key in the `Dict`, called `Q_nd`:
=#
pumpdata["250rpm"]["Q_nd"] = DimensionlessParameter.(pumpdata["250rpm"]["Q"]/(omega*D^3))

#=
Similarly for the head rise coefficient, which we call head_nd here:
=#
pumpdata["250rpm"]["head_nd"] = DimensionlessParameter.(pumpdata["250rpm"]["head"]*g/(omega^2*D^2))

#=
Let's do the same at 500 and 1000 rpm
=#
omega = 500u"rpm"
pumpdata["500rpm"]["Q_nd"] = DimensionlessParameter.(pumpdata["500rpm"]["Q"]/(omega*D^3))
pumpdata["500rpm"]["head_nd"] = DimensionlessParameter.(pumpdata["500rpm"]["head"]*g/(omega^2*D^2));
#-
omega = 1000u"rpm"
pumpdata["1000rpm"]["Q_nd"] = DimensionlessParameter.(pumpdata["1000rpm"]["Q"]/(omega*D^3))
pumpdata["1000rpm"]["head_nd"] = DimensionlessParameter.(pumpdata["1000rpm"]["head"]*g/(omega^2*D^2));

#=
Now plot these
=#

scatter(xlabel=L"$Q/\omega D^3$",ylabel=L"$gh/\omega^2 D^2$",legend=:bottomleft,xlims=(0,0.1),ylims=(0,0.25))
scatter!(pumpdata["250rpm"]["Q_nd"],pumpdata["250rpm"]["head_nd"],label="250 rpm")
scatter!(pumpdata["500rpm"]["Q_nd"],pumpdata["500rpm"]["head_nd"],label="500 rpm")
scatter!(pumpdata["1000rpm"]["Q_nd"],pumpdata["1000rpm"]["head_nd"],label="1000 rpm")

#=
The data collapse onto a single curve. This reveals the power of dimensional
analysis! One would follow up by generating a statistical correlation
between the two dimensionless parameters. Then, we could use this correlation
to predict performance at any pump speed or pump diameter!
=#

#savefig("pumpplot.pdf")


text:
const u8vec4 = Vec{4, UInt8}
struct Vertex
    position::Vec3f0
    color::u8vec4
end

function testVertexInputState(inputStateRef::Ref)
    attributes = Vector{vk.VkVertexInputAttributeDescription}()
    bindings = Vector{vk.VkVertexInputBindingDescription}()
    lava.binding(attributes, bindings, UInt32(0), Vertex, :position)
    lava.binding(attributes, bindings, UInt32(0), Vertex, :color)
    vertexInputState = lava.PipelineVertexInputStateCreateInfo(attributes = attributes, bindings = bindings)
    inputStateRef[] = vertexInputState
    return true
end


text:
using Cropbox
using Garlic
using Test

import Dates

@testset "garlic" begin
    r = simulate(Garlic.Model;
        config=Garlic.Examples.AoB.KM_2014_P2_SR0,
        stop="calendar.count",
        snap=s -> Dates.hour(s.calendar.time') == 12,
    )
    @test r.leaves_initiated[end] > 0
    visualize(r, :DAP, [:leaves_appeared, :leaves_mature, :leaves_dropped]) |> display # Fig. 3.D
    visualize(r, :DAP, :green_leaf_area) |> display # Fig. 4.D
    visualize(r, :DAP, [:leaf_mass, :bulb_mass, :total_mass]) |> display
end


text:
"""
    absorption_spectrum(α, len)

Computes the absorption spectrum for the given length `len` in centimeters.
The vector `α` should be calculated using the corresponding [`α`](@ref) function.
"""
absorption_spectrum(α::AbstractVector, len) = @. 1 - exp(-α*len)

"""
    transmittance_spectrum(α, len)

Computes the transmittance spectrum for the given length `len` in centimeters.
The vector `α` should be calculated using the corresponding [`α`](@ref) function.
"""
transmittance_spectrum(α::AbstractVector, len) = @. exp(-α*len)

"""
    optical_depth(α, len)

Computes the optical depth for the given length `len` in centimeters.
The vector `α` should be calculated using the corresponding [`α`](@ref) function.
"""
optical_depth(α::AbstractVector, len) = @. α*len


text:
"""
    wetted!(area_size, drop_pos, maxheight, height, t; hthresh = 0.055)

Measures the wetted area and maximal height of the film at time step `t`.
"""
function wetted!(area_size, maxheight, height, t; hthresh = 0.055)
    area_size[t] = length(findall(height .> hthresh))
    maxheight[t] = maximum(height)

    return nothing
end

function wetted!(area_size, state::LBM_state; hthresh = 0.055)
    push!(area_size, length(findall(state.height .> hthresh)))

    return nothing
end
"""
    fluid_dry!(fluid, dummy, height, t; hthresh = 0.055)

Tracks the location of the thin film as a boolean field.
"""
function fluid_dry!(fluid, dummy, height, t; hthresh = 0.055)
    dummy .= false
    dummy[height .> hthresh] .= true
    fluid[t, :] .= vec(dummy)

    return nothing
end

"""
    t0(;hᵦ=0.07, γ=0.01, μ=1/6, θ=1/6)

Computes a characteristic time scale for an spinodally dewetting film.

# Arguments

- `hᵦ :: Float64`: height at which the disjoining pressure vanishes
- `γ :: Float64`: surface tension value
- `μ :: Float64`: kinematic viscosity, same as dynamic for ρ=1
- `θ :: Float64`: highest contact angle given as radiant, e.g. θ=π/9 for 20 degrees

# Mathematics

The charateristic time scale ``t_0`` is set using the surface tension as well as the disjoining pressure ``\\Pi(h)`` as

`` t_0 = \\frac{3}{\\gamma h_0^3 q_0^4}, \\quad q_0 = \\frac{\\Pi'(h_0)}{2\\gamma}, ``

where ``\\Pi'(h_0) = \\frac{\\partial\\Pi}{\\partial h}\\bigg|_{h_0}`` is the derivative of the disjoining pressure with respect to some characteristic height.

# References

- [Mecke, Rauscher](https://iopscience.iop.org/article/10.1088/0953-8984/17/45/042/meta)
"""
function t0(;hᵦ=0.07, γ=0.01, μ=1/6, θ=1/6)
    qsq = hᵦ * (1 - cospi(θ)) * (2 - 3 * hᵦ)
    charT = 3 * μ / (γ * qsq^2)

    return charT
end

"""
    snapshot!(snap, field, t; dumping)

Makes a copy of the current state of an input array `in` and saves the vectorized values as a column to `out`.

Function that fills a preallocated array out with a time series of system snap shots of e.g. the height field `h`.

# Arguments

- `snap :: Array{Number,2}`: Array that stores the snap shots as columns
- `field :: Array{Number,2}`: Input argument, e.g. ``h(\\mathbf{x},t)``
- `t :: Int`: The current time step
- `dumping :: Int`: Sampling frequency

# Examples

```jldoctest
julia> using Swalbe, Test

julia> h1 = reshape(collect(1:25),5,5); h2 = reshape(collect(5:5:125),5,5);

julia> snapshot = zeros(2, 25);

julia> Swalbe.snapshot!(snapshot,h1,10,dumping=10)

julia> Swalbe.snapshot!(snapshot,h2,20,dumping=10)

julia> @test all(h1 .== reshape(snapshot[1,:],5,5))
Test Passed
  Expression: all(h1 .== reshape(snapshot[1, :], 5, 5))

julia> @test all(h2 .== reshape(snapshot[2,:],5,5))
Test Passed
  Expression: all(h2 .== reshape(snapshot[2, :], 5, 5))
```
# References

See also: The [scripts folder](https://github.com/Zitzeronion/Swalbe.jl/tree/master/scripts)
"""
function snapshot!(snap, field, t; dumping = 1000)
    if t % dumping == 0
        snap[t÷dumping, :] .= vec(Array(field))
    end

    return nothing
end

"""
    surfacearea!(area_lv, red_energy, height, θ, ∇hx, ∇hy, dgrad, surface,)

Measures the surface area of the liquid vapor interface and the reduced surface energy.

# Arguments

- `area_lv :: Vector{Float64}`: array to store the computed liquid vapor area
- `red_energy :: Vector{Float64}`: array the stores the computed reduced surface energy
- `height :: Matrix{Float64}`: current height configuration
- `θ :: Matrix{Float64}`: contact angle field distribution
- `∇hx :: Matrix{Float64}`: height gradient x-component
- `∇hy :: Matrix{Float64}`: height gradient y-component
- `dgrad :: Array{Float64,3}`: dummy array to store derivatives
- `surface :: Matrix{Float64}`: array that computes locally the liquid vapor surface area
- `t :: Int`: current time step
- `hthresh :: Float64`: height threshold below which the substrate is considered *dry*

"""
function surfacearea!(area_lv, red_energy, height, θ::Float64, ∇hx, ∇hy, dgrad, surface, t; htresh = 0.055)
    ∇f_simple!(∇hx, ∇hy, height, dgrad)
    surf = 0.0
    surface .= sqrt.(∇hx.^2 .+ ∇hy.^2 .+ 1)
    surf = sum(surface[height .> htresh])
    area_lv[t] = surf
    red_energy[t] = surf - length(height[height .> htresh]) * cospi(θ)

    return nothing
end

function surfacearea!(area_lv, red_energy, height, θ::Matrix, ∇hx, ∇hy, dgrad, surface, t; htresh = 0.055)
    ∇f_simple!(∇hx, ∇hy, height, dgrad)
    surf = 0.0
    surface .= sqrt.(∇hx.^2 .+ ∇hy.^2 .+ 1)
    surf = sum(surface[height .> htresh])
    area_lv[t] = surf
    red_energy[t] = surf - sum(cospi.(θ[height .> htresh]))

    return nothing
end

"""
    ∇f_simple!(outputx, outputy, f, dgrad)

Simple gradient calculation for the differential surface area.
"""
function ∇f_simple!(outputx, outputy, f, dgrad)
    fip, fjp, fim, fjm, fipjp, fimjp, fimjm, fipjm = Swalbe.viewneighbors(dgrad)
    # Straight elements  j+1, i+1, i-1, j-1
    circshift!(fip, f, (1,0))
    circshift!(fjp, f, (0,1))
    circshift!(fim, f, (-1,0))
    circshift!(fjm, f, (0,-1))
    # Diagonal elements
    circshift!(fipjp, f, (1,1))
    circshift!(fimjp, f, (-1,1))
    circshift!(fimjm, f, (-1,-1))
    circshift!(fipjm, f, (1,-1))
    # In the end it is just a weighted sum...
    outputx .= -1/3 .* (fip .- fim) .- 1/12 .* (fipjp .- fimjp .- fimjm .+ fipjm)
    outputy .= -1/3 .* (fjp .- fjm) .- 1/12 .* (fipjp .+ fimjp .- fimjm .- fipjm)

    return nothing
end


text:
using FunProject
using Documenter

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

makedocs(;
    modules=[FunProject],
    authors="Nivedita Rethnakar et al.",
    repo="https://github.com/nivupai/FunProject.jl/blob/{commit}{path}#{line}",
    sitename="FunProject.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nivupai.github.io/FunProject.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nivupai/FunProject.jl",
    devbranch="main",
)


text:
macro checkerr(expr)
  quote
    ex = $(esc(expr))
    if haskey(ex, "error")
      error("API Error: $(ex["error"]["message"])")
    else
      ex["result"]
    end
  end
end

function showres(results::Vector)
  println("Result List:")
  for res in results
    showres(res; offset=1)
  end
end
function showres(result::Dict; offset=0)
  println(repeat(" ", offset), "Result:")
  for (key,val) in result
    println(repeat(" ", offset+1), key, " => ", val)
  end
end

macro showres(expr)
  :(showres($(esc(expr))))
end


text:
isprime(n, primes) = isempty(filter(p -> mod(n, p) == 0 ,primes))
function sieve(limit)
    limit < 2 && return []
    primes = []
    foreach(i -> isprime(i, primes) && push!(primes, i), 2:limit)
    primes
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
# using Pkg
# Pkg.add(["DataFrames", "DataFramesMeta", "CSV"])
# Pkg.add(PackageSpec(url="https://github.com/FNj/Hose.jl.git", rev="indexing"))

lib_load_time = @elapsed using DataFrames, CSV, DataFramesMeta, Hose, Statistics, StatsBase
# import Base.Meta.@dump
const dir = @__DIR__

optionally_convert_to_df(df::DataFrame) = df
optionally_convert_to_df(df) = df |> DataFrame

file_load_compilation_time = @elapsed CSV.File("data/hflights.csv", missingstring="NA") |> optionally_convert_to_df
file_load_time = @elapsed const hflights = CSV.File("data/hflights.csv", missingstring="NA") |> optionally_convert_to_df
file_load_compilation_time -= file_load_time

# size(hflights)
# names(hflights)
# head(hflights)
# describe(hflights)

# macro sample(df, n::Int)
#     esc(:($df[sample(1:nrow($df), $n), :]))
# end

# macro sample(df, q::Real)
#     esc(:($df[sample(1:nrow($df), round(Int, $q * nrow($df))), :]))
# end

function process(df)
    @hose df |>
    @transform(Speed = :Distance ./ :AirTime .* 60) |>
    @select(:Month, :ArrDelay, :Speed) |>
    dropmissing |>
    @by(:Month, AvgDelay = mean(:ArrDelay), MaxSpeed = maximum(:Speed))
end

compile_time = @elapsed process(hflights[1:5, :])
const N = 20
times = zeros(N)
# GC.enable(false)
# for i=1:N
#     GC.gc()
#     sleep(0.01)
#     times[i] = @elapsed process(hflights)
# end
# GC.enable(true)
#
# mean_time_no_GC = mean(times)

for i=1:N
    times[i] = @elapsed process(hflights)
end

mean_time = mean(times)
compile_time -= mean_time

process_strrev(df) = @hose df |> @transform(TailNum2=reverse.(:TailNum))

string_reversal_compile_time = @elapsed process_strrev(hflights)
string_reversal_time = @elapsed process_strrev(hflights)

string_reversal_compile_time -= string_reversal_time

println("""Stats:
- Libraries loading time: $lib_load_time s,
- File loading compilation time: $file_load_compilation_time s,
- File loading time: $file_load_time s,
- Group by compilation time: $compile_time s,
- Group by execution time: $mean_time s,
- String reversal compilation time: $string_reversal_compile_time s,
- String reversal execution time: $string_reversal_time s.""")

using Dates
task_names = vcat(["Library loading"], fill("File loading", 2), fill("Split-apply-combine", 2), fill("String reversal", 2))
task_types = vcat(["Loading"], repeat(["Compilation", "Execution"], 3))
task_exec_times = [lib_load_time, file_load_compilation_time, file_load_time, compile_time, mean_time, string_reversal_compile_time, string_reversal_time]
stats_df = DataFrame(DateTime=fill(Dates.now(), length(task_names)), Language=fill("Julia", length(task_names)), Dataset=fill("HFlights", length(task_names)), TaskNames=task_names, TaskTypes=task_types, TaskExecTimes=task_exec_times)

try
    vcat(CSV.read("logs/log.csv") |> DataFrame, stats_df)
catch
    stats_df
end |> CSV.write("logs/log.csv")

max_stars_repo_name: RevoltBI/Julia-Meetup

max_issues_repo_name: RevoltBI/Julia-Meetup

max_forks_repo_name: RevoltBI/Julia-Meetup


text:
import ..UncertainValues.CertainValue

Base.truncate(v::CertainValue) = v
function Base.truncate(v::CertainValue, constraint::TruncateMaximum)
    if v.value > constraint.max
        msg = "Truncating $v with $constraint failed\n"
        msg2 = "Need value < constraint.max, got $v < $(constraint.max)"
        throw(ArgumentError(msg * msg2))
    else
        return v
    end
end

function Base.truncate(v::CertainValue, constraint::TruncateMinimum)
    if v.value < constraint.min
        msg = "Truncating $v with $constraint failed\n"
        msg2 = "Need value > constraint.min, got $v > $(constraint.min)"
        throw(ArgumentError(msg * msg2))
    else
        return v
    end
end

function Base.truncate(v::CertainValue, constraint::TruncateRange)
    if v.value < constraint.min
        msg = "Truncating $v with $constraint failed\n"
        msg2 = "Need value > constraint.min, got $v > $(constraint.min)"
        throw(ArgumentError(msg * msg2))
    elseif v.value > constraint.max
        msg = "Truncating $v with $constraint failed\n"
        msg2 = "Need value < constraint.max, got $v < $(constraint.max)"
        throw(ArgumentError(msg * msg2))
    else
        return v
    end
end

truncate(v::CertainValue, s::TruncateLowerQuantile) = v
truncate(v::CertainValue, s::TruncateUpperQuantile) = v
truncate(v::CertainValue, s::TruncateQuantiles) = v
truncate(v::CertainValue, s::TruncateStd) = v


text:
######################################################################
# printing.jl
######################################################################
using Printf

function print_header()
	print_dashes()
	println(" iter  |ergodic score |control score |total score |direc deriv |step size")
	print_dashes()
end
function print_dashes()
	println("--------------------------------------------------------------------------")
end

function step_report(i::Int, es::Float64, cs::Float64, ts::Float64, dd::Float64, step_size::Float64)
	@printf " %-7i" i
	@printf " %-14.7f" es
	@printf " %-14.7f" cs
	@printf " %-12.7f" ts
	@printf " %-12.7f" dd
	@printf " %-7.5f" step_size
	println()
end


text:
# ---
# title: 1466. Reorder Routes to Make All Paths Lead to the City Zero
# id: problem1466
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/reorder-routes-to-make-all-paths-lead-to-the-city-zero/description/>
# hidden: true
# ---
#
# There are `n` cities numbered from `0` to `n-1` and `n-1` roads such that
# there is only one way to travel between two different cities (this network
# form a tree). Last year, The ministry of transport decided to orient the roads
# in one direction because they are too narrow.
#
# Roads are represented by `connections` where `connections[i] = [a, b]`
# represents a road from city `a` to `b`.
#
# This year, there will be a big event in the capital (city 0), and many people
# want to travel to this city.
#
# Your task consists of reorienting some roads such that each city can visit the
# city 0\. Return the **minimum** number of edges changed.
#
# It's **guaranteed** that each city can reach the city 0 after reorder.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/13/sample_1_1819.png)**
#
#
#
#     Input: n = 6, connections = [[0,1],[1,3],[2,3],[4,0],[4,5]]
#     Output: 3
#     Explanation: Change the direction of edges show in red such that each node can reach the node 0 (capital).
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/13/sample_2_1819.png)**
#
#
#
#     Input: n = 5, connections = [[1,0],[1,2],[3,2],[3,4]]
#     Output: 2
#     Explanation: Change the direction of edges show in red such that each node can reach the node 0 (capital).
#
# **Example 3:**
#
#
#
#     Input: n = 3, connections = [[1,0],[2,0]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 5 * 10^4`
#   * `connections.length == n-1`
#   * `connections[i].length == 2`
#   * `0 <= connections[i][0], connections[i][1] <= n-1`
#   * `connections[i][0] != connections[i][1]`
#
#
## @lc code=start
using LeetCode

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