text:
# `@doc` needed for using a `raw` string as a docstring; raw needed because
# of the escape characters in the regex.
@doc raw"""
    Document{T<:NamedTuple}

Represents a single string document. This object has two fields,

* `text::String`
* `metadata::T`

The `text` is automatically processed by applying the replacements
from [`AUTOMATIC_REPLACEMENTS`](@ref) and
adding a space to the end of the document.
"""
struct Document{T<:NamedTuple}
    text::String
    metadata::T

    function Document(text::AbstractString, metadata::T) where {T}
        check_keys(T)

        # Add a final space to ensure that the last word is recognized
        # as a word boundary.
        new_text = apply_replacements(text) * " "
        return new{T}(new_text, metadata)
    end
end

Document(text::AbstractString) = Document(text, NamedTuple())

function apply_replacements(str::AbstractString)
    # Apply automatic replacements
    # using https://github.com/JuliaLang/julia/issues/29849#issuecomment-449535743
    return foldl(replace, AUTOMATIC_REPLACEMENTS; init=str)
end

abstract type AbstractQuery end

"""
    QueryMatch{Q<:AbstractQuery,Doc<:Document,D,I}

Represents a match for an `AbstractQuery`, with four fields:

* `query::Q`: the query itself
* `document::Doc`: the [`Document`](@ref) which was matched to
* `distance::D`: the distance of the match
* `indices::I`: the indices of where in the `document` the match occurred.

"""
struct QueryMatch{Q<:AbstractQuery,Doc<:Document,D,I}
    query::Q
    document::Doc
    distance::D
    indices::I
end

"""
    Query <: AbstractQuery

A query to search for an exact match of a string,
with one field:

* `text::String`

The `text` is automatically processed by applying the replacements
from [`AUTOMATIC_REPLACEMENTS`](@ref).
"""
struct Query <: AbstractQuery
    text::String
    Query(str::AbstractString) = new(apply_replacements(str))
end

"""
    Base.match(query::AbstractQuery, document::Document)

Looks for a match for `query` in `document`. Returns either `nothing`
if no match is found, or a [`QueryMatch`](@ref) object.
"""
Base.match(query::AbstractQuery, document::Document)

"""
    match_all(query::AbstractQuery, document::Document)

Looks for all matches for `query` in the document. Returns a
`Vector` `QueryMatch` objects corresponding to all of the matches found.
"""
match_all(query::AbstractQuery, document::Document)

function Base.match(Q::Query, R::Document)
    if (length(R.text) < length(Q.text)) || (length(Q.text) == 0)
        return nothing
    end
    inds = findfirst(Q.text, R.text)
    inds === nothing && return nothing
    return QueryMatch(Q, R, 0, inds)
end

function match_all(Q::Query, R::Document)
    if (length(R.text) < length(Q.text)) || (length(Q.text) == 0)
        all_inds = UnitRange[]
    else
        all_inds = findall(Q.text, R.text)
        if all_inds === nothing
            all_inds = UnitRange[]
        end
    end
    return [QueryMatch(Q, R, 0, inds) for inds in all_inds]
end

struct Or{S<:Tuple} <: AbstractQuery
    subqueries::S
end

function Base.match(Q::Or, R::Document)
    for subquery in Q.subqueries
        m = match(subquery, R)
        m !== nothing && return m
    end
    return nothing
end

function match_all(Q::Or, R::Document)
    return reduce(vcat, (match_all(subquery, R) for subquery in Q.subqueries))
end

# Specializations to combine Or's
function Or(q1::Or, q2::AbstractQuery)
    return Or((q1.subqueries..., q2))
end

Or(q1::AbstractQuery, q2::Or) = Or(q2, q1)

function Or(q1::Or, q2::Or)
    return Or((q1.subqueries..., q2.subqueries...))
end

Or(q1::AbstractQuery, q2::AbstractQuery) = Or((q1, q2))

Base.:(|)(q1::AbstractQuery, q2::AbstractQuery) = Or(q1, q2)

"""
    FuzzyQuery{D,T} <: AbstractQuery

A query to search for an fuzzy match of a string,
with three fields:

* `text::String`: the text to match
* `dist::D`: the distance measure to use; defaults to `DamerauLevenshtein()`
* `threshold::T`: the maximum threshold allowed for a match; defaults to 2.

The `text` is automatically processed by applying the replacements
from [`AUTOMATIC_REPLACEMENTS`](@ref).
"""
struct FuzzyQuery{D,T} <: AbstractQuery
    text::String
    dist::D
    threshold::T
    function FuzzyQuery(str::AbstractString, dist::D, threshold::T) where {D,T}
        return new{D,T}(apply_replacements(str), dist, threshold)
    end
end

FuzzyQuery(str::String) = FuzzyQuery(str, DamerauLevenshtein(), 2)

function dist_with_threshold(dist::DamerauLevenshtein, str1, str2, max_dist)
    return DamerauLevenshtein(max_dist)(str1, str2)
end

function dist_with_threshold(dist::Levenshtein, str1, str2, max_dist)
    return Levenshtein(max_dist)(str1, str2)
end

# Maybe to be upstreamed?
# https://github.com/matthieugomez/StringDistances.jl/issues/29

"""
    _findmin(s1, s2, dist::Partial; max_dist) -> d, inds

StringDistances' `Partial(dist)(s1, s2, max_dist)` returns
the value `d`, the closest partial match between the two strings, up to a maximum
distance `max_dist` (if no match is found less than `max_dist`, then
`max_dist+1` is returned). `_findmin` returns the same value, but also
returns the first set of indices at which an optimal partial match was found.
"""
function _findmin(s1, s2, dist::Partial; max_dist)
    s1, s2 = StringDistances.reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 == len2 && return dist_with_threshold(dist.dist, s1, s2, max_dist),
           firstindex(s2):lastindex(s2)
    out = max_dist + 1
    len1 == 0 && return out, 1:0
    out_idx = 0
    for (i, x) in enumerate(qgrams(s2, len1))
        curr = dist_with_threshold(dist.dist, s1, x, max_dist)
        out_idx = ifelse(curr < out, i, out_idx)
        out = min(out, curr)
        max_dist = min(out, max_dist)
    end
    return out, nextind(s2, 0, out_idx):nextind(s2, 0, out_idx + len1 - 1)
end

function Base.match(Q::FuzzyQuery, R::Document)
    # We assume the document text is longer than the query text
    length(R.text) < length(Q.text) && return nothing

    dist, inds = _findmin(Q.text, R.text, Partial(Q.dist); max_dist=Q.threshold)
    return dist <= Q.threshold ? QueryMatch(Q, R, dist, inds) : nothing
end

"""
    _findall(s1, s2, dist::Partial; max_dist) -> Vector{Tuple{Int,UnitRange}}

Returns all of the matches within `max_dist` found by `dist`, returning tuples
of the distance along with the indices of the match.
"""
function _findall(s1, s2, dist::Partial; max_dist)
    s1, s2 = StringDistances.reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    matches = Tuple{Int,UnitRange}[]

    len1 == 0 && return matches

    if len1 == len2
        curr = dist_with_threshold(dist.dist, s1, s2, max_dist)
        if curr <= max_dist
            push!(matches, (curr, firstindex(s2):lastindex(s2)))
        end
        return matches
    end

    for (i, x) in enumerate(qgrams(s2, len1))
        curr = dist_with_threshold(dist.dist, s1, x, max_dist)
        if curr <= max_dist
            inds = nextind(s2, 0, i):nextind(s2, 0, i + len1 - 1)
            push!(matches, (curr, inds))
        end
    end
    return matches
end

# This takes as input all matches and partitions them into runs
# of overlapping matches. Then chooses the best match from each run.
# Let's say there are matches at 1:3, 2:3, 2:4, 4:5 and 6:7
# Then we only have 2 runs, one involving indices 1:5 and one from 6:7
# Note that 1:3 and 4:5 do not intersect, yet they are in the same run.
# (A math way to say this is that we choose the best representative from the equivalence class of connected components).
# Why do it this way? Well, it's the best I thought of so far.
# Another way would be to say "two matches are the same if they share > 50% of the same indices".
# That seems slightly more reasonable (avoding the situation in the example above where 1:3 and 4:5 have no overlap)
# But it's actually hard to do sensibly, since you have have say 1:3, 2:4, 3:5 and 4:6. Then 1:3 and 2:4 share 2 indices,
# and 2:4 and 3:5 share two indices, and 3:5 and 4:6 share two indices. So which do you keep? Could resolve them one at at time left-to-right
# but then you're actually doing the same thing as the "run of connected components" except requiring a 50% overlap and can run into
# the same issue where you discard disconnected matches. For exmaple what if 4:6 is the best match followed by 3:5 then 2:4?
# Then you have 1:3 and 2:4, keep 2:4. Then you have 2:4 and 3:5, keep 3:5. Then
# you have 3:5 and 4:6, keep 4:6. Now we've just chosen the best match from all of them, but didn't keep 1:3 even though it has no overlap with 4:6.
function non_overlapping_matches(matches)
    length(matches) <= 1 && return matches
    non_overlapping_matches = Tuple{Int,UnitRange}[]
    best_curr_match = first(matches)
    inds_so_far = best_curr_match[2]
    for m in matches
        dist, inds = m
        if first(inds) <= last(inds_so_far)
            inds_so_far = first(inds_so_far):last(inds)
            if dist < best_curr_match[1]
                best_curr_match = m
            end
        else
            push!(non_overlapping_matches, best_curr_match)
            best_curr_match = m
            inds_so_far = inds
        end
    end
    push!(non_overlapping_matches, best_curr_match)
    return non_overlapping_matches
end

function match_all(Q::FuzzyQuery, R::Document)
    matches = _findall(Q.text, R.text, Partial(Q.dist); max_dist=Q.threshold)
    matches_no_overlap = non_overlapping_matches(matches)
    return [QueryMatch(Q, R, m...) for m in matches_no_overlap]
end


text:
using Plots

function bracket_minimum(f, x=0; s=1e-2, k=2.0)
    a, ya = x, f(x)
    b, yb = a + s, f(a + s)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        s = -s
    end
    while true
        c, yc = b + s, f(b + s)
        if yc > yb
            return a < c ? (a, c) : (c, a)
        end
        a, ya, b, yb = b, yb, c, yc
        s *= k
    end
end

φ = MathConstants.φ
function golden_section_search(f, a, b, n)
    ρ = φ-1
    d = ρ * b + (1 - ρ)*a
    yd = f(d)
    for i = 1 : n-1
        c = ρ*a + (1 - ρ)*b
        yc = f(c)
        if yc < yd
            b, d, yd = d, c, yc
        else
            a, b = b, c
        end
    end
    return a < b ? (a, b) : (b, a)
end

#definicao da funcao
f = x-> -(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x    #multiplicada por -1 para encontrar o minimo
g = x-> -(-(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x) #funcao original

#ponto inicial
initial_point = 0
#busca do intervalo
a,c = bracket_minimum(f,initial_point)

#numero de iteracoes
n=(c-a)/(10^-3*log(φ))

am,bm = golden_section_search(f,a,c,n)
#(am,f(am))
#am,bm

println("x* = ", am)

plot(g,a-1.,c+1.,lw=2,draw_arrow="true",label="Function")
scatter!([a,c],[g(a),g(c)],label="interval")
scatter!([am,bm],[g(am),g(bm)],label="x*")

ext: jl


text:
abstract type           Factorization end
struct LeftFactored  <: Factorization end
struct RightFactored <: Factorization end

"""
    chomsky_normal_form(tree, factor=RightFactored(), labelf)

Convert a tree into chomsky normal form.

# Arguments

- `tree`: the constituency tree to convert
- `factor=RightFactored()`: can be `LeftFactored()` or `RightFactored()`
- `labelf`: function (called like `labelf(tree, children)`)

# Returns

- a tree of the same type as its argument
"""
function chomsky_normal_form(tree, factor::Factorization=RightFactored(), labelf=nltk_factored_label)
    T, node, ch = typeof(tree), label(tree), children(tree)
    N = length(ch)
    cnf = t -> chomsky_normal_form(t, factor, labelf)
    if isleaf(tree) || isterminal(tree)
        return tree
    elseif N == 1
        return cnf(ch[1])
    elseif N == 2
        return T(node, cnf.(ch))
    else
        return T(node, cnf.(factorize(tree, factor, labelf)))
    end
end

function factorize(tree, lf::LeftFactored, labelf)
    T, ch = typeof(tree), children(tree)
    N = length(ch)
    l_children = ch[1:N-1]
    l_label = labelf(tree, l_children)
    if length(l_children) > 2
        l_children = factorize(T(label(tree), l_children), lf, labelf)
    end
    r = last(ch)
    return [T(l_label, l_children), T(label(r), children(r))]
end
function factorize(tree, rf::RightFactored, labelf)
    T, ch = typeof(tree), children(tree)
    r_children = ch[2:end]
    r_label = labelf(tree, r_children)
    if length(r_children) > 2
        r_children = factorize(T(label(tree), r_children), rf, labelf)
    end
    l = first(ch)
    return [T(label(l), children(l)), T(r_label, r_children)]
end

"""
    collapse_unary(tree, labelf=unary_label; collapse_pos=false, collapse_root=false)

Transform a parse tree by collapsing all single-child nodes.

# Arguments

- `tree`: parse tree to transform
- `labelf`: function called to create new label representing the collapsed nodes.

# Keywords

- `collapse_pos=false`: whether to collapse `(POS word)` nodes
- `collapse_root=false`: whether to collapse the top-level root node

# Returns

- a transformed tree of the same type as its argument
"""
function collapse_unary(tree, labelf=unary_label; collapse_pos=false, collapse_root=false)
    isleaf(tree) && return tree
    T, n = typeof(tree), label(tree)
    kw = (collapse_pos=collapse_pos, collapse_root=collapse_root)
    if !collapse_root
        f = x -> collapse_unary(x, labelf; collapse_pos=collapse_pos, collapse_root=true)
        return T(n, f.(children(tree)))
    elseif collapse_pos && isterminal(tree)
        return children(tree)[1]
    elseif length(children(tree)) == 1
        if isterminal(tree)
            return tree
        else
            node = labelf(n, children(tree)[1])
            ch = children(children(tree)[1])
            return collapse_unary(T(node, ch), labelf; kw...)
        end
    else
        return T(n, [collapse_unary(c, labelf; kw...) for c in children(tree)])
    end
end

function nltk_factored_label(parent, children, child_sep = "|", parent_sep = "^")
    return label(parent) * child_sep * "<" * join(label.(children), "-") * ">"
end

parent_label(parent, args...) = parent
prime_label(parent, args...) = parent * "'"

unary_label(parent, children) = join([label(parent), label(children)], "+")


text:
# m5.2s.jl

using Pkg, DrWatson

using StanQuap
using StatisticalRethinking

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])

stan5_2 = "
data {
      int N;
      vector[N] divorce_s;
      vector[N] marriage_s;
}
parameters {
      real a;
      real bM;
      real<lower=0> sigma;
}
model {
      vector[N] mu = a + bM * marriage_s;
      a ~ normal( 0 , 0.2 );
      bM ~ normal( 0 , 0.5 );
      sigma ~ exponential( 1 );
      divorce_s ~ normal( mu , sigma );
}
";

data = (N=size(df, 1), divorce_s=df.Divorce_s, marriage_s=df.Marriage_s)
init = (a=1.0, bM=1.0, sigma=10.0)
q5_2s, m5_2s, o5_2s = stan_quap("m5.2s", stan5_2; data, init);

if !isnothing(q5_2s)
      quap5_2s_df = sample(q5_2s)
      precis(quap5_2s_df)
end

# End m5.2s.jl


text:
using Test

@testset "Examples" begin
  examples_path = joinpath(@__DIR__, "..", "examples")
  example_files = filter(starts_and_ends_with("ex_", ".jl"), readdir(examples_path))
  for file in example_files
    file_path = joinpath(examples_path, file)
    println("Testing file $(file_path)")
    empty!(ARGS)
    push!(ARGS, "false")
    res = include(file_path)
    @test isnothing(res) || all(isnothing, res)
    empty!(ARGS)
  end
end


text:
"""
    solve_karger(model, experiment, difftensors)
Solve the finite pulse Karger model (FPK) using precomputed effective diffusion tensors
`difftensors`.
"""
function solve_karger(model, experiment, difftensors)

    # Measure function evalutation time
    starttime = Base.time()

    # Extract input parameters
    @unpack mesh, T₂, ρ, κ = model
    @unpack directions, sequences, values, values_type = experiment.gradient
    @unpack odesolver, timestep = experiment.karger
    dt = timestep

    # Deduce sizes
    ncompartment, nboundary = size(mesh.facets)
    ndirection = size(directions, 2)
    nsequence = length(sequences)
    namplitude = length(values)

    qvalues, bvalues = get_values(experiment.gradient)

    # Volumes
    volumes = get_cmpt_volumes(mesh)

    # Compute surface areas
    surface_areas = spzeros(ncompartment, nboundary)
    for icmpt = 1:ncompartment
        # Finite elements
        points = mesh.points[icmpt]
        facets = mesh.facets[icmpt, :]

        # Surface area
        for iboundary = 1:nboundary
            if !isempty(facets[iboundary])
                surface_areas[icmpt, iboundary], =
                    get_mesh_surface(points, facets[iboundary])
            end
        end
    end

    # Associate surface areas to the compartment interfaces they measure
    tau_inv = spzeros(ncompartment, ncompartment)
    for iboundary = 1:nboundary
        inds = surface_areas[:, iboundary].nzind
        if length(inds) == 2
            i₁ = inds[1]
            i₂ = inds[2]
            tmp = κ[iboundary] * surface_areas[i₁, iboundary]
            tau_inv[i₁, i₂] = tmp
            tau_inv[i₂, i₁] = tmp
        end
    end
    tau_inv = tau_inv ./ volumes'
    A = tau_inv - spdiagm((tau_inv * volumes) ./ volumes)

    # Relaxation tensor
    R = spdiagm(1 ./ T₂)

    # Initial signal
    S₀ = volumes .* ρ

    # Update function for linear ODE operator
    function update_func(J, u, p, t)
        @unpack A, ADC_diag, R, f, q = p
        J .= A .- (integral(f, t)^2 * q^2) .* ADC_diag .- R
    end

    # Allocate output arrays
    signal = zeros(ncompartment, namplitude, nsequence, ndirection)
    signal_allcmpts = zeros(namplitude, nsequence, ndirection)
    itertimes = zeros(namplitude, nsequence, ndirection)

    # Iterate over gradient amplitudes, sequences and directions
    for idir = 1:ndirection, iseq = 1:nsequence, iamp = 1:namplitude

        # Measure iteration time
        itertime = Base.time()

        # Gradient amplitude
        q = qvalues[iamp, iseq]
        b = bvalues[iamp, iseq]

        # Time profile
        f = sequences[iseq]
        TE = echotime(f)

        # Gradient direction
        g = directions[:, idir]

        # Display state of iterations
        println("Computing Karger signal:")
        println("  Direction $idir of $ndirection: g = $g")
        println("  Sequence  $iseq of $nsequence: f = $f")
        println("  Amplitude $iamp of $namplitude: q = $q, b = $b")

        ADC_diag = spdiagm([g' * D * g for D ∈ difftensors[:, iseq]])

        p = (; A, ADC_diag, R, f, q)
        J_prototype = A - ADC_diag - R
        J = DiffEqArrayOperator(J_prototype, update_func = update_func)

        prob = ODEProblem(J, S₀, (0.0, TE), p)
        sol = solve(prob, odesolver, dt = dt)

        # Save final signal
        signal[:, iamp, iseq, idir] = sol.u[end]

        # Store timing
        itertimes[iamp, iseq, idir] = Base.time() - itertime
    end

    signal_allcmpts = sum(signal, dims = 1)[1, :, :, :]

    totaltime = Base.time() - starttime

    (; signal, signal_allcmpts, itertimes, totaltime)
end


text:
module Observations

using Oceananigans
using Oceananigans: short_show, fields
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z
using Oceananigans.Grids: pop_flat_elements, topology, halo_size, on_architecture
using Oceananigans.Fields
using Oceananigans.Utils: SpecifiedTimes
using Oceananigans.Architectures: arch_array, architecture
using JLD2

import Oceananigans.Fields: set!

abstract type AbstractObservation end

include("normalization.jl")

"""
    SyntheticObservations{F, G, T, P, M} <: AbstractObservation

A time series of synthetic observations generated by Oceananigans.jl's simulations
gridded as Oceananigans.jl fields.
"""
struct SyntheticObservations{F, G, T, P, M, N} <: AbstractObservation
     field_time_serieses :: F
                    grid :: G
                   times :: T
                    path :: P
                metadata :: M
           normalization :: N
end

observation_names(ts::SyntheticObservations) = keys(ts.field_time_serieses)

"""
    observation_names(ts_vector::Vector{<:SyntheticObservations})

Return a Set representing the union of all names in `ts_vector`.
"""
function observation_names(ts_vector::Vector{<:SyntheticObservations})
    names = Set()
    for ts in ts_vector
        push!(names, observation_names(ts)...)
    end

    return names
end

obs_str(ts::SyntheticObservations) = "SyntheticObservations of $(keys(ts.field_time_serieses)) on $(short_show(ts.grid))"
obs_str(ts::Vector{<:SyntheticObservations}) = "Vector of SyntheticObservations of $(keys(ts[1].field_time_serieses)) on $(short_show(ts[1].grid))"

tupleit(t) = try
    Tuple(t)
catch
    tuple(t)
end

const not_metadata_names = ("serialized", "timeseries")

read_group(group::JLD2.Group) = NamedTuple(Symbol(subgroup) => read_group(group[subgroup]) for subgroup in keys(group))
read_group(group) = group

function with_size(new_size, old_grid)

    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    new_size = pop_flat_elements(new_size, topo)
    halo = pop_flat_elements(halo_size(old_grid), topo)

    new_grid = RectilinearGrid(architecture(old_grid), eltype(old_grid);
        size = new_size,
        x = x, y = y, z = z,
        topology = topo,
        halo = halo)

    return new_grid
end

location_guide = Dict(:u => (Face, Center, Center),
                      :v => (Center, Face, Center),
                      :w => (Center, Center, Face))

function infer_location(field_name)
    if field_name in keys(location_guide)
        return location_guide[field_name]
    else
        return (Center, Center, Center)
    end
end

function SyntheticObservations(path; field_names, normalize = IdentityNormalization, times = nothing, field_time_serieses = nothing, grid_size = nothing)
    field_names = tupleit(field_names)

    field_time_serieses === nothing && (field_time_serieses = NamedTuple(name => FieldTimeSeries(path, string(name); times) for name in field_names))

    grid = first(field_time_serieses).grid
    times = first(field_time_serieses).times
    boundary_conditions = first(field_time_serieses).boundary_conditions

    if !isnothing(grid_size) # regrid! field_time_serieses on a new grid

        new_field_time_serieses = Dict()

        # Re-grid the data in `field_time_serieses`
        for (field_name, ts) in zip(keys(field_time_serieses), field_time_serieses)

            grid = with_size(grid_size, grid)

            LX, LY, LZ = infer_location(field_name)

            new_ts = FieldTimeSeries{LX,LY,LZ}(CPU(), grid, times, boundary_conditions)

            # Loop over time steps to re-grid each constituent field in `field_time_series`
            for n = 1:length(times)
                regrid!(new_ts[n], ts[n])
            end

            new_field_time_serieses[field_name] = new_ts
        end

        field_time_serieses = NamedTuple(new_field_time_serieses)
    end

    # validate_data(fields, grid, times) # might be a good idea to validate the data...
    file = jldopen(path)
    metadata = NamedTuple(Symbol(group) => read_group(file[group]) for group in filter(n -> n ∉ not_metadata_names, keys(file)))
    close(file)

    normalization = Dict(name => normalize(field_time_serieses[name]) for name in keys(field_time_serieses))

    return SyntheticObservations(field_time_serieses, grid, times, path, metadata, normalization)
end

observation_times(observation) = observation.times

function observation_times(obs::Vector)
    @assert all([o.times ≈ obs[1].times for o in obs]) "Observations must have the same times."

    return observation_times(first(obs))
end

#####
##### set! for simulation models and observations
#####

function set!(model, ts::SyntheticObservations, index=1)
    # Set initial condition
    for name in keys(fields(model))

        model_field = fields(model)[name]

        if name ∈ keys(ts.field_time_serieses)
            ts_field = ts.field_time_serieses[name][index]
            set!(model_field, ts_field)
        else
            set!(model_field, 0)
        end
    end

    return nothing
end

"""
    column_ensemble_interior(observations::Vector{<:SyntheticObservations}, field_name, time_indices::Vector, N_ens)

Returns an `N_cases × N_ens × Nz` array of the interior of a field `field_name` defined on a
`OneDimensionalEnsembleGrid` of size `N_cases × N_ens × Nz`, given a list of `SyntheticObservations` objects
containing the `N_cases` single-column fields at time index in `time_index`.
"""
function column_ensemble_interior(observations::Vector{<:SyntheticObservations}, field_name, time_index, ensemble_size)
    zeros_column = zeros(size(observations[1].field_time_serieses[1].grid))
    Nt = length(observation_times(observations))

    batch = []
    for observation in observations
        fts = observation.field_time_serieses
        if field_name in keys(fts) && time_index <= Nt
            push!(batch, interior(fts[field_name][time_index]))
        else
            push!(batch, zeros_column)
        end
    end

    batch = cat(batch..., dims = 2) # (n_batch, n_z)
    ensemble_interior = cat([batch for i = 1:ensemble_size]..., dims = 1) # (ensemble_size, n_batch, n_z)

    return ensemble_interior
end

function set!(model, observations::Vector{<:SyntheticObservations}, index = 1)

    for name in keys(fields(model))

        model_field = fields(model)[name]

        field_ts_data = column_ensemble_interior(observations, name, index, model.grid.Nx)

        arch = architecture(model_field)

        # Reshape `field_ts_data` to the size of `model_field`'s interior
        reshaped_data = arch_array(arch, reshape(field_ts_data, size(model_field)))

        # Sets the interior of field `model_field` to values of `reshaped_data`
        model_field .= reshaped_data
    end

    return nothing
end

struct FieldTimeSeriesCollector{G, D, F, T}
    grid :: G
    field_time_serieses :: D
    collected_fields :: F
    times :: T
end

"""
    FieldTimeSeriesCollector(collected_fields, times; architecture=CPU())

Returns a `FieldTimeSeriesCollector` for `fields` of `simulation`.
`fields` is a `NamedTuple` of `AbstractField`s that are to be collected.
"""
function FieldTimeSeriesCollector(collected_fields, times; architecture=CPU())

    grid = on_architecture(architecture, first(collected_fields).grid)
    field_time_serieses = Dict{Symbol, Any}()

    for name in keys(collected_fields)
        field = collected_fields[name]
        LX, LY, LZ = location(field)
        field_time_series = FieldTimeSeries{LX, LY, LZ}(grid, times)
        field_time_serieses[name] = field_time_series
    end

    # Convert to NamedTuple
    field_time_serieses = NamedTuple(name => field_time_serieses[name] for name in keys(collected_fields))

    return FieldTimeSeriesCollector(grid, field_time_serieses, collected_fields, times)
end

function (collector::FieldTimeSeriesCollector)(simulation)
    for field in collector.collected_fields
        compute!(field)
    end

    current_time = simulation.model.clock.time
    time_index = findfirst(t -> t >= current_time, collector.times)

    for name in keys(collector.collected_fields)
        field_time_series = collector.field_time_serieses[name]
        set!(field_time_series[time_index], collector.collected_fields[name])
    end

    return nothing
end

function initialize_simulation!(simulation, observations, time_series_collector, time_index = 1)
    set!(simulation.model, observations, time_index)

    times = observation_times(observations)

    initial_time = times[time_index]
    simulation.model.clock.time = initial_time
    simulation.model.clock.iteration = 0
    simulation.model.timestepper.previous_Δt = Inf

    # Zero out time series data
    for time_series in time_series_collector.field_time_serieses
        time_series.data .= 0
    end

    simulation.callbacks[:data_collector] = Callback(time_series_collector, SpecifiedTimes(times...))

    simulation.stop_time = times[end]

    return nothing
end

summarize_metadata(::Nothing) = ""
summarize_metadata(metadata) = keys(metadata)

Base.show(io::IO, ts::SyntheticObservations) =
    print(io, "SyntheticObservations with fields $(propertynames(ts.field_time_serieses))", '\n',
              "├── times: $(ts.times)", '\n',
              "├── grid: $(short_show(ts.grid))", '\n',
              "├── path: \"$(ts.path)\"", '\n',
              "├── metadata: ", summarize_metadata(ts.metadata), '\n',
              "└── normalization: $(short_show(ts.normalization))")

end # module


text:
function for_with_break()
    for i in (0:4-1)
        if i == 2
            break
        end
        println(join([i], " "))
    end
end

function for_with_continue()
    for i in (0:4-1)
        if i == 2
            continue
        end
        println(join([i], " "))
    end
end

function for_with_else()
    for i in (0:4-1)
        println(join([i], " "))
    end
end

function while_with_break()
    i = 0
    while true
        if i == 2
            break
        end
        println(join([i], " "))
        i += 1
    end
end

function while_with_continue()
    i = 0
    while i < 5
        i += 1
        if i == 2
            continue
        end
        println(join([i], " "))
    end
end

function main()
    for_with_break()
    for_with_continue()
    while_with_break()
    while_with_continue()
end

main()


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libsass"], :libsass_so),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/piever/SassBuilder/releases/download/v3.5.5-1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/libsass.v3.5.5.aarch64-linux-gnu.tar.gz", "6c04e15fe9485207cd7cc1327b87b883c85d9f0e0e5f8aafd81b5244f9db4a80"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/libsass.v3.5.5.aarch64-linux-musl.tar.gz", "60bc726ea57191c03acf8ac96c6f0f5630515b73ca0a7b86299909b328efa448"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/libsass.v3.5.5.arm-linux-gnueabihf.tar.gz", "024822bbf15fc82c18966c0e16d17fab8792757328ceb44518752e52c02f6807"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/libsass.v3.5.5.arm-linux-musleabihf.tar.gz", "7ead3aa331cd2082ef3e14f81101c4abebe2cd6d1ae4133813d1c01ab064d34c"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/libsass.v3.5.5.i686-linux-gnu.tar.gz", "0afce3dced23efa55a7ddc762cb3064fabcea7fcd37ef3c72f2a290f38d24307"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/libsass.v3.5.5.i686-linux-musl.tar.gz", "b46d2f1c92dfa3ff4fc028380da00ac1d0b79b019d5694d6f79fa378e7208e92"),
    Windows(:i686) => ("$bin_prefix/libsass.v3.5.5.i686-w64-mingw32.tar.gz", "bad6acf124d6649b8ec4b021cdfd88aecfc4151734f2c5ce1683fda07f6fb236"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/libsass.v3.5.5.powerpc64le-linux-gnu.tar.gz", "fcfc2d96cbead84577cdf2d8e9810612037d8a67fc4861bb7e92d2009edee588"),
    MacOS(:x86_64) => ("$bin_prefix/libsass.v3.5.5.x86_64-apple-darwin14.tar.gz", "a2f8768639f49a7444fd815b6a4e572c650495059272c70af73dec149fde7ee0"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/libsass.v3.5.5.x86_64-linux-gnu.tar.gz", "28316732cbe46c1047db3068824468b8a620b750b311bcef1b430c00de81092e"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/libsass.v3.5.5.x86_64-linux-musl.tar.gz", "25c7925c9b57fce19a646c040484bd5f6a21603a27561469f319b96c462f8241"),
    FreeBSD(:x86_64) => ("$bin_prefix/libsass.v3.5.5.x86_64-unknown-freebsd11.1.tar.gz", "12a2b5739b186a14189b369101fac2241790104697cc7480ded0783f757231a1"),
    Windows(:x86_64) => ("$bin_prefix/libsass.v3.5.5.x86_64-w64-mingw32-gcc7.tar.gz", "c24fb85e81b8e2329651592406b4d371bea1e24f9dc7157f1534660945582ec9"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
const isWin64 = Sys.iswindows() && Sys.WORD_SIZE == 64
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    @static if isWin64
        install(dl_info...; prefix=prefix, force=true, verbose=verbose, ignore_platform=true)
    else
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end
end
# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using DataStructures: Deque

mutable struct TierraModel
    time::UInt64

    next_key::UInt64

    organism_keys::Vector{UInt64}
    organisms::Dict{UInt64, TierrianOrganism}

    slice_index::UInt16

    memory::Vector{UInt8}
    memory_size::UInt32
    free_blocks::Vector{FreeMemoryBlock}
    used_memory::UInt32

    reaper_queue::Vector{UInt64}

    logger::Logger

    function TierraModel(ancestor_program::Vector{TierrianInstruction})
        TierraModel(ancestor_program, DoNothingLogger())
    end
    function TierraModel(ancestor_program::Vector{TierrianInstruction}, logger::Logger)
        TierraModel(convert_instructions(ancestor_program), logger)
    end

    function TierraModel(ancestor_program::Vector{UInt8})
        TierraModel(ancestor_program, DoNothingLogger())
    end
    function TierraModel(ancestor_program::Vector{UInt8}, logger::Logger)
        @assert 0 < length(ancestor_program)

        memory_size = typemax(UInt16) + 1
        memory = zeros(UInt8, memory_size)

        ancestor_length = convert(UInt16, length(ancestor_program))
        ancestor_key = 1
        ancestor = TierrianOrganism(UInt16(0), ancestor_length)
        ancestor.key = ancestor_key
        ancestor.hash = bytes2hex(sha256(ancestor_program))

        memory[1:ancestor_length] = ancestor_program
        remaining_memory_size = convert(UInt32, memory_size - ancestor_length)
        free_block = FreeMemoryBlock(ancestor_length, remaining_memory_size)

        new(0, 2, [ancestor_key], Dict(ancestor_key => ancestor),
            UInt16(1), memory, memory_size, [free_block],
            ancestor_length, [ancestor_key], logger)
    end
end

function add_organism!(model::TierraModel, program::Vector{TierrianInstruction})
    add_organism!(model, convert_instructions(program))
end

function add_organism!(model::TierraModel, program::Vector{UInt8})
    length = convert(UInt16, Base.length(program))
    start_address = allocate_free_memory!(model, length)

    if start_address == -1
        return
    end

    new_organism = TierrianOrganism(start_address, length)

    key = add_organism!(model, new_organism)

    write_memory(model, new_organism, new_organism.start_address, program)
    set_hash!(new_organism, model)

    return key
end

function add_organism!(model::TierraModel, organism::TierrianOrganism)
    key = model.next_key

    organism.key = key
    model.organisms[key] = organism
    push!(model.organism_keys, key)
    push!(model.reaper_queue, key)

    model.next_key += 1

    organism.time_birth = model.time

    return key
end

function remove_organism!(model::TierraModel, key; remove_from_reaper_queue = true)
    organism_to_be_removed = model.organisms[key]

    clear_memory!(model, organism_to_be_removed.start_address, organism_to_be_removed.length)
    if organism_to_be_removed.has_daughter
        clear_memory!(model, organism_to_be_removed.daughter_address, organism_to_be_removed.daughter_length)
    end

    delete!(model.organisms, key)
    delete!(model.organism_keys, key)

    if remove_from_reaper_queue
        delete!(model.reaper_queue, key)
    end

    log_death(model.logger, model, organism_to_be_removed)
end

function set_hash!(organism::TierrianOrganism, model::TierraModel)
    program = read_memory(model, organism.start_address, organism.length)
    organism.hash = bytes2hex(sha256(program))
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# This is an example on how to use MIMO with AbstractSDRs assuming that SDR support several channels.
# The following example has been built for a USRP e310 and should work for any multiple channel USRP.

# ----------------------------------------------------
# --- Package dependencies
# ----------------------------------------------------
using AbstractSDRs
using Plots
using FFTW

# ----------------------------------------------------
# --- Radio parameters
# ----------------------------------------------------
radioType   = :uhd      # We target UHDBindings here
carrierFreq = 2410e6     # Carrier frequency (ISM band)
samplingRate = 2e6      # Not high but issue with e310 stream :(
gain        = 50        # In dB

# ----------------------------------------------------
# --- MIMO and board specificities
# ----------------------------------------------------
# We need to specify the channels and the board as we use a USRP => Note that the way the parameters are set is based on how it is done at UHD level. You can have a look at the mimo example provided by UHD.
# This may vary depending on the chosen hardware
# /!\ We need to specify how many antenna we want. Be sure that the configuration is handled by the hardware you have otherwise UHD will raise an error (and maybe segfault ?)
nbAntennaRx = 2
nbAntennaTx = 0
# Antenna and board config
channels    = [0;1]     # First channel is the first path
subdev      = "A:0 A:1" # e310 board names
# For the antenna, we need to specify the dictionnary of the allocated antenna for both Tx and Rx. In our case we will only do full Rx mode so we only specify the :Rx key
antennas    = Dict(:Rx => ["TX/RX";"RX2"])

# ----------------------------------------------------
# --- Open radio
# ----------------------------------------------------
radio = openSDR(
                radioType,
                carrierFreq,
                samplingRate,
                gain ;
                nbAntennaTx, nbAntennaRx,
                channels,
                antennas,
                subdev)

# --- Receive a buffer
# We specify the size of each buffer
# We will have one buffer per channel so a Vector of Vector
# sigVect[1] and sigVect[2] will have the same size (of nbSamples) each for a channel
nbSamples = 4096
sigVect = recv(radio,nbSamples)

# ----------------------------------------------------
# --- Plot the spectrum
# ----------------------------------------------------
# Get the rate from the radio using the accessor
s = getSamplingRate(radio)
# X axis
xAx = ((0:nbSamples-1)./nbSamples .- 0.5) * s
# PSD
y = 10*log10.(abs2.(fftshift(fft(sigVect[1]))))
plt = plot(xAx,y,label="First Channel")
y = 10*log10.(abs2.(fftshift(fft(sigVect[2]))))
plot!(plt,xAx,y,label="Second Channel")
xlabel!("Frequency [Hz]")
ylabel!("Magnitude [dB]")
display(plt)

# ----------------------------------------------------
# --- Close radio
# ----------------------------------------------------
close(radio)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct InstrumentationLibrary <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function InstrumentationLibrary(; kwargs...)
        obj = new(meta(InstrumentationLibrary), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct InstrumentationLibrary
const __meta_InstrumentationLibrary = Ref{ProtoMeta}()
function meta(::Type{InstrumentationLibrary})
    ProtoBuf.metalock() do
        if !isassigned(__meta_InstrumentationLibrary)
            __meta_InstrumentationLibrary[] = target = ProtoMeta(InstrumentationLibrary)
            allflds = Pair{Symbol,Union{Type,String}}[
                :name=>AbstractString,
                :version=>AbstractString,
            ]
            meta(
                target,
                InstrumentationLibrary,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_InstrumentationLibrary[]
    end
end
function Base.getproperty(obj::InstrumentationLibrary, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :version
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct AnyValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AnyValue(; kwargs...)
        obj = new(meta(AnyValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct AnyValue (has cyclic type dependency)
const __meta_AnyValue = Ref{ProtoMeta}()
function meta(::Type{AnyValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AnyValue)
            __meta_AnyValue[] = target = ProtoMeta(AnyValue)
            allflds = Pair{Symbol,Union{Type,String}}[
                :string_value => AbstractString,
                :bool_value   => Bool,
                :int_value    => Int64,
                :double_value => Float64,
                :array_value  => "ArrayValue",
                :kvlist_value => "KeyValueList",
                :bytes_value  => Vector{UInt8},
            ]
            oneofs = Int[1, 1, 1, 1, 1, 1, 1]
            oneof_names = Symbol[Symbol("value")]
            meta(
                target,
                AnyValue,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                oneofs,
                oneof_names,
            )
        end
        __meta_AnyValue[]
    end
end
function Base.getproperty(obj::AnyValue, name::Symbol)
    if name === :string_value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :bool_value
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :int_value
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :double_value
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :array_value
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :kvlist_value
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :bytes_value
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct ArrayValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ArrayValue(; kwargs...)
        obj = new(meta(ArrayValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ArrayValue (has cyclic type dependency)
const __meta_ArrayValue = Ref{ProtoMeta}()
function meta(::Type{ArrayValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ArrayValue)
            __meta_ArrayValue[] = target = ProtoMeta(ArrayValue)
            allflds = Pair{Symbol,Union{Type,String}}[:values=>Base.Vector{AnyValue}]
            meta(
                target,
                ArrayValue,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ArrayValue[]
    end
end
function Base.getproperty(obj::ArrayValue, name::Symbol)
    if name === :values
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AnyValue}
    else
        getfield(obj, name)
    end
end

mutable struct KeyValueList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function KeyValueList(; kwargs...)
        obj = new(meta(KeyValueList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct KeyValueList (has cyclic type dependency)
const __meta_KeyValueList = Ref{ProtoMeta}()
function meta(::Type{KeyValueList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_KeyValueList)
            __meta_KeyValueList[] = target = ProtoMeta(KeyValueList)
            allflds = Pair{Symbol,Union{Type,String}}[:values=>"Base.Vector{KeyValue}"]
            meta(
                target,
                KeyValueList,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_KeyValueList[]
    end
end
function Base.getproperty(obj::KeyValueList, name::Symbol)
    if name === :values
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct KeyValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function KeyValue(; kwargs...)
        obj = new(meta(KeyValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct KeyValue (has cyclic type dependency)
const __meta_KeyValue = Ref{ProtoMeta}()
function meta(::Type{KeyValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_KeyValue)
            __meta_KeyValue[] = target = ProtoMeta(KeyValue)
            allflds =
                Pair{Symbol,Union{Type,String}}[:key=>AbstractString, :value=>AnyValue]
            meta(
                target,
                KeyValue,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_KeyValue[]
    end
end
function Base.getproperty(obj::KeyValue, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AnyValue
    else
        getfield(obj, name)
    end
end

export AnyValue,
    ArrayValue,
    KeyValueList,
    KeyValue,
    InstrumentationLibrary,
    AnyValue,
    ArrayValue,
    KeyValueList,
    KeyValue


text:
module MachineLearningWithJulia

# package code goes here

end # module


text:
# Dependencies
# Pkg.add("WignerSymbols")
# Pkg.add("LinearAlgebra")
# Pkg.add("JuliennedArrays")
# Pkg.add("Statistics")
# Pkg.add("GeometryTypes")
# Pkg.add("Makie")
# Pkg.add("Colors")
# Pkg.add("SpecialFunctions")

using WignerSymbols
using LinearAlgebra
using SparseArrays

export Spherical, harmonic, cleb, rotate, rotate!

# Our first step is to make spherical tensors semi-native. Note this only works with integer valued spherical harmonics
# This takes two steps.  First we define and indexer for the spherical harmonics (the MRange type)
# This is a special unit range type which defines the size of the tensor
# It signals to the compile that broadcasting must be initiated with another spherical tensor.
# Plus the methods for getting Base to understand how it works
# Then we define the spherical tensor type and methods for accessing it's index in the natural way

struct MRange <: AbstractUnitRange{Int}
    start::Int
    stop::Int
end

function MRange(L)
    MRange(-L,L)
end

# definition
struct Spherical{L} <: AbstractArray{Complex{Float64},1}
    x::Array{Complex{Float64},1}
end

# standard constructor
function Spherical(x)
    Spherical{Int((length(x)-1)/2)}(x)
end

# initialization: unit vector in m = 0
function Spherical{L}() where L
    Spherical(Complex{Float64}.( (-L:L) .== 0))
end

# interface with Base
Base.axes(y::Spherical{L}) where L = (MRange(L),)
Base.size(y::Spherical{L}) where L = (2*L+1,)
function Base.getindex(y::Spherical{L}, m::Int) where L
    y.x[L+m+1]
end
function Base.setindex!(y::Spherical{L}, v ,m::Int) where L
    y.x[L+m+1] = v
end
function Base.similar(A::AbstractArray, T::Type, shape::Tuple{MRange,Vararg{MRange}})
    Spherical(A)
end
function Base.similar(f::Any, shape::Tuple{MRange,Vararg{MRange}})
    #something didn't work here from the recipe book, might be a bug
    Spherical(first(shape))
end

# Basic Operations on Spherical Tensors

# clebsch product into L3 irrep
function cleb(y1::Spherical{L1},y2::Spherical{L2},L3::Int) where {L1, L2}
    if abs(L1 - L2) <= L3 <= L1 + L2
        return Spherical{L3}(
            [ sum( clebschgordan(L1, m1, L2, m3 - m1, L3, m3) * y1[m1] * y2[m3 - m1]
                for m1 in -L1:L1 if abs(m3 - m1) <= L2)
            for m3 in -L3:L3 ] # comprehension and index work makes this easy!
         )
    else
        nothing
    end
end

# Generator matrices and rotations
# lowerJ |l,m> = sqrt(l(l+1) - m(m-1)) |l,m-1>
# raiseJ |l,m> = sqrt(l(l+1) - m(m+1)) |l,m+1>
function raiseJ(L)
    rows = L .+ 1 .- ((- L):(L-1))
    cols = L .+ 1 .- ((1-L):L)
    vals = sqrt.( L*(L+1) .- (-L:(L-1)) .* ((-L+1):L))
    sparse(rows,cols, vals, 2L+1, 2L+1)
end

function lowerJ(L)
    rows = L .+ 1 .- ((1-L):L)
    cols = L .+ 1 .- ((- L):(L-1))
    vals = sqrt.( L*(L+1) .- (-L:(L-1)) .* ((-L+1):L))
    sparse(rows,cols, vals, 2L+1, 2L+1)
end

function Jz(L)
    rows = L .+ 1 .+ (-L:L)
    cols = L .+ 1 .+ (-L:L)
    vals = (-L:L)
    sparse(rows,cols, vals, 2L+1, 2L+1)
end

# the representation of jx and jy in irrep L
Jx(L) = (raiseJ(L) + lowerJ(L))/2
Jy(L) = (raiseJ(L) - lowerJ(L))/(2*im)

function rotate!(y::Spherical{L}, nvec::Array) where L
    # rotates the spherical tensor around the axis defined by nvec
    # angle is given by the magnitude of nvec.
    mat = exp( Matrix(  im .* (
        Jx(L) .* nvec[1] .+
        Jy(L) .* nvec[2] .+
        Jz(L) .* nvec[3]) ))
    y = Spherical{L}(mat * y.x); #must be made dense to get this method to work
end

function rotate(y::Spherical{L}, nvec::Array) where L
    # rotates the spherical tensor around the axis defined by nvec
    # angle is given by the magnitude of nvec.
    mat = exp( Matrix( im .* (
        Jx(L) .* nvec[1] .+
        Jy(L) .* nvec[2] .+
        Jz(L) .* nvec[3]) ))
    return Spherical{L}(mat * y.x);
    # must be made dense to use Base.exp method
    # a sparse matrix exponentiator from a package could speed this up,
    # not sure how much we will need this function (it's already fast)
end

function zAxisTo(new::Array)
    # return the smalles nvec that takes the z axis to point new
    b = normalize(new)
    if abs(b[3]) < 1 # have to deal with the degenerate case
        return normalize([-b[2], b[1], 0])*acos(b[3])
    else
        return [1, 0, 0]*acos(b[3])
    end
end

function harmonic(y::Spherical{L}, point::Array) where L
    # returns the value of the spherical harmonic at point projected onto the unit sphere
    # val::Complex = Σ_m Ylm(point) y[m] = f(point) c.f. W.K. Tung 8.3 - 14
    rotate(y::Spherical{L}, zAxisTo(point))[0]*sqrt((2L + 1)/(4π))
    # The m=0 component of the spherical tensor with point defining the z-axis
    # Is normalization over sphere ∫dΩ conj(y(Ω)))* y(Ω) = dot(y,y)
end

function harmonic(L::Int, point::Array)
    # returns a spherical tensor composed of spherical harmonics
    # Ylm(phat) = y[m] c.f. W.K. Tung 8.3 - 15
    conj(rotate(Spherical{L}(), -zAxisTo(point))*sqrt((2L + 1)/(4π)))
    # orthogonality relations ∫dΩ conj(y(Ω)[m])* y(Ω)[m'] = δ(m',m)
end


text:
module IndexTransforms

using ArrayInterface
using ArrayInterface: offsets, offset1, known_length, known_offsets, known_size, buffer,
    known_step, ndims_index, static_first, static_step, static_last, device
using ArrayInterface: OptionallyStaticUnitRange, OptionallyStaticStepRange, OptionallyStaticRange,
    CanonicalInt, AbstractDevice, AbstractCPU, CPUPointer
using LinearAlgebra
using SparseArrays
using SparseArrays: getcolptr, AbstractSparseMatrixCSC
using Static

using Base: @propagate_inbounds, AbstractCartesianIndex, Fix2, isdone, ReshapedArray, Generator

abstract type CoordinateTransform{I,O} end

abstract type IndexTransform{I,O} <: CoordinateTransform{I,O} end

"""
    UnkownTransform{N}

Represents an unkown transform that throws an error when called.
"""
struct UnkownTransform{N} <: CoordinateTransform{N,N} end

(::UnkownTransform)(x) = error("Attempt to execute an unkown transform.")

"""
    IdentityTransform

Returns the input unchanged.
"""
struct IdentityTransform{N} <: CoordinateTransform{N,N} end

(::IdentityTransform)(x) = x

include("utils.jl")
include("Size.jl")
include("TransformedIndex.jl")
include("Permute.jl")
include("StrideTransform.jl")
include("LinearView.jl")
include("View.jl")
include("Reshape.jl")
include("sparse.jl")
include("ComposedTransform.jl")

@inline function _to_linear(x)
    StrideTransform(ArrayInterface.size_to_strides(ArrayInterface.size(x), static(1)), offsets(x))
end

@inline _to_cartesian(x) = reshape(Size(x))

## constructors
CoordinateTransform{N}(x) where {N} = UnkownTransform{N}()
CoordinateTransform{1}(::DenseArray{<:Any,1}) = IdentityTransform{1}()
CoordinateTransform{1}(::DenseArray{<:Any,N}) where {N} = IdentityTransform{1}()
CoordinateTransform{N}(A::DenseArray{<:Any,N}) where {N} = StrideTransform(A)

CoordinateTransform{2}(::Diagonal) = Permute{(1,1),(1,)}()

CoordinateTransform{2}(x::SparseArrays.AbstractSparseMatrixCSC) = CompactSparseColumn(x)
CoordinateTransform{2}(x::AbstractSparseVector) = SparseLinear(x)

CoordinateTransform{1}(x::ReshapedArray) = IdentityTransform{1}()
CoordinateTransform{N}(x::ReshapedArray) where {N} = _to_linear(x)

# TODO should we only define index constructors for explicit types?
CoordinateTransform{1}(x::AbstractRange) = OffsetTransform(offset1(x) - static(1))

## SubArray
CoordinateTransform{N}(x::SubArray{<:Any,N}) where {N} = View(getfield(x, :indices))
@inline function CoordinateTransform{1}(x::SubArray{T,N}) where {T,N}
    if N === 1
        return View(getfield(x, :indices))
    else
        return ComposedTransform(View(getfield(x, :indices)), _to_cartesian(x))
    end
end
CoordinateTransform{1}(x::Base.FastContiguousSubArray) = LinearView(getfield(x, :offset1), static(1))
CoordinateTransform{1}(x::Base.FastSubArray) = LinearView(getfield(x, :offset1), getfield(x, :stride1))

## PermutedDimsArray
CoordinateTransform{1}(::PermutedDimsArray{<:Any,1,(1,),(1,)}) = IdentityTransform{1}()
@inline CoordinateTransform{N}(::PermutedDimsArray{<:Any,N,perm,iperm}) where {N,perm,iperm} = Permute(static(perm), static(iperm))
function CoordinateTransform{1}(x::PermutedDimsArray{<:Any,N,perm,iperm}) where {N,perm,iperm}
    ComposedTransform(CoordinateTransform{N}(x), _to_cartesian(x))
end

## Transpose/Adjoint{Real}
CoordinateTransform{1}(::Transpose{<:Any,<:AbstractVector{<:Any}}) = IdentityTransform{1}()
CoordinateTransform{2}(::Transpose{<:Any,<:AbstractVector{<:Any}}) = Permute((static(2), static(1)), (static(2),))
CoordinateTransform{2}(::Transpose{<:Any,<:AbstractMatrix{<:Any}}) = Permute((static(2), static(1)), (static(2), static(1)))
function CoordinateTransform{1}(x::Transpose{<:Any,<:AbstractMatrix{<:Any}})
    ComposedTransform(Permute((static(2), static(1)), (static(2), static(1))), _to_cartesian(x))
end
CoordinateTransform{1}(::Adjoint{<:Real,<:AbstractVector{<:Real}}) = IdentityTransform{1}()
CoordinateTransform{2}(::Adjoint{<:Real,<:AbstractVector{<:Real}}) = Permute((static(2), static(1)), (static(2),))
CoordinateTransform{2}(::Adjoint{<:Real,<:AbstractMatrix{<:Real}}) = Permute((static(2), static(1)), (static(2), static(1)))
function CoordinateTransform{1}(x::Adjoint{<:Real,<:AbstractMatrix{<:Real}})
    ComposedTransform(Permute((static(2), static(1)), (static(2), static(1))), _to_cartesian(x))
end

include("layout.jl")

end


text:
using CancerImagingArchive, DataFrames
using Test

#######
# SETUP
#######

# Use global variable for filenames because we want to delete them if they already exist
zip_file = "test.zip"
dicom_file = "test.dcm"
csv_file = "test.csv"
json_file = "test.json"

for file in [zip_file, dicom_file, csv_file, json_file]
    rm(file, force = true)
end

# Helper function for comparing CSV/DataFrames vs JSON/DictionaryArrays
function compare_csv_vs_json(csv, json; max_names = Inf)
    names_in_csv = find_names_in_csv(csv)
    names_in_json = find_names_in_json(json)
    @test length(json) == length(csv[!, names_in_csv[1]])
    @test length(names_in_json) == length(names_in_csv)
    if length(names_in_json) > max_names
        # Some cases have too many rows and they are not directly comparable because of types
        names_in_json = names_in_json[1:max_names]
        names_in_csv = names_in_csv[1:max_names]
    end
    for (idx, json_element) in enumerate(json)
        for name in names_in_csv
            if ismissing(csv[!, name][idx])
                # Missing entries lead to this messy situation. Just check if both json/csv are empty/missing
                @test isempty(json_element) || isempty(json_element[string(name)])
                @test ismissing(csv[!, name][idx])
            else
                @test csv[!, name][idx] == json_element[string(name)]
            end
        end
    end
    return nothing
end

function find_names_in_csv(csv)
    names_in_csv = names(csv)
    # AnnotationsFlag field exists in csv but not json, so remove it for comparisons
    filter!(name -> name != "AnnotationsFlag", names_in_csv)
    return names_in_csv
end

function find_names_in_json(json_array)
    # JSON names get ignored if the entry is missing so we can't just do collect(keys(json_array[1]))
    found_names = []
    num_names = 0
    for json in json_array
        cur_names = keys(json)
        if length(cur_names) > num_names
            found_names = collect(cur_names)
            num_names = length(found_names)
        end
    end
    return found_names
end

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

@testset "Queries - Collection" begin
    @test_throws ErrorException tcia_collections(format = "unknown")
    collections_csv = tcia_collections()
    collections_json = tcia_collections(format = "json")
    @test length(collections_json) > 90
    compare_csv_vs_json(collections_csv, collections_json)
end

@testset "Queries - Modalities" begin
    @test length( tcia_modalities(collection = "TCGA-GBM", format = "json") ) > 2
    @test length( tcia_modalities(bodypart = "BREAST", format = "json") ) > 5
    compare_csv_vs_json(
        tcia_modalities(collection = "TCGA-GBM", bodypart = "BRAIN"),
        tcia_modalities(collection = "TCGA-GBM", bodypart = "BRAIN", format = "json"))
end

@testset "Queries - BodyParts" begin
    @test "BRAIN" in tcia_bodyparts(modality = "MR").BodyPartExamined
    compare_csv_vs_json(
        tcia_bodyparts(collection = "CPTAC-HNSCC"),
        tcia_bodyparts(collection = "CPTAC-HNSCC", format = "json"))
end

@testset "Queries - Manufacturers" begin
    compare_csv_vs_json(
        tcia_manufacturers(collection = "TCGA-KICH", modality = "MR"),
        tcia_manufacturers(collection = "TCGA-KICH", modality = "MR", format = "json"))
    compare_csv_vs_json(
        tcia_manufacturers(bodypart = "BREAST"),
        tcia_manufacturers(bodypart = "BREAST", format = "json"))
end

@testset "Queries - Patients" begin
    compare_csv_vs_json(
        tcia_patients(collection = "TCGA-THCA"),
        tcia_patients(collection = "TCGA-THCA", format = "json"))

    # Following criteria should only find one patient
    found_patient = tcia_patients_by_modality(collection = "ACRIN-FLT-Breast", modality = "OT")
    @test length(found_patient.PatientID) == 1
    @test found_patient.PatientID[1] == "ACRIN-FLT-Breast_066"

    # Following criteria should find at least two patients
    new_gbm_patients = tcia_newpatients(collection = "TCGA-GBM", date = "2015-01-01", format = "json")
    @test length(new_gbm_patients) > 1
end

@testset "Queries - Studies" begin
    # The CSV version requires a few manual changes, so we do them first
    studies_csv = tcia_studies(collection = "TCGA-SARC")
    # 1. Convert the date to plain strings so that they can be compared with the json version
    studies_csv.StudyDate = string.(studies_csv.StudyDate)
    # 2. Remove the escape characters in the string. These occur in the study description
    for (idx, description) in enumerate(studies_csv.StudyDescription)
        studies_csv.StudyDescription[idx] = replace(description, "\\" => "")
    end

    compare_csv_vs_json(
        studies_csv,
        tcia_studies(collection = "TCGA-SARC", format = "json"))

    # Following criteria should find at least three series
    @test length(tcia_newstudies(collection="TCGA-GBM", date="2015-01-01", format="json")) > 2
end

@testset "Queries - Series" begin
    compare_csv_vs_json(
        tcia_series(collection = "TCGA-THCA"),
        tcia_series(collection = "TCGA-THCA", format = "json"), max_names = 3)
    compare_csv_vs_json(
        tcia_series(study = "1.3.6.1.4.1.14519.5.2.1.3023.4024.298690116465423805879206377806"),
        tcia_series(study = "1.3.6.1.4.1.14519.5.2.1.3023.4024.298690116465423805879206377806", format = "json"), max_names = 3)
    compare_csv_vs_json(
        tcia_series(bodypart = "CHEST", modality = "CT", manufacturer = "TOSHIBA"),
        tcia_series(bodypart = "CHEST", modality = "CT", manufacturer = "TOSHIBA", format = "json"), max_names = 3)

    # Can not use compare_csv_vs_json() on tcia_series_size() because TotalSizeInBytes has different types
    dce_series_json = tcia_series_size(series = "1.3.6.1.4.1.14519.5.2.1.4591.4001.241972527061347495484079664948", format="json")[1]
    @test dce_series_json["TotalSizeInBytes"] == "149149266.000000"
    dce_series_csv = tcia_series_size(series = "1.3.6.1.4.1.14519.5.2.1.4591.4001.241972527061347495484079664948")
    @test dce_series_csv.TotalSizeInBytes[1] ≈ 149149266
    @test dce_series_csv.ObjectCount[1] == dce_series_json["ObjectCount"] == 1120
end

@testset "Queries - SOP" begin
    compare_csv_vs_json(
        tcia_sop(series = "1.3.6.1.4.1.14519.5.2.1.4591.4001.241972527061347495484079664948"),
        tcia_sop(series = "1.3.6.1.4.1.14519.5.2.1.4591.4001.241972527061347495484079664948", format = "json"))
end

@testset "Data Download" begin
    patient_studies = tcia_studies(collection = "TCGA-THCA")
    chosen_study = patient_studies.StudyInstanceUID[1]
    imaging_series = tcia_series(study = chosen_study)
    chosen_series = imaging_series.SeriesInstanceUID[1]
    series_sops = tcia_sop(series = chosen_series)
    chosen_sop = series_sops.SOPInstanceUID[1]

    tcia_images(series = chosen_series, file = zip_file)
    @test isfile(zip_file)
    @test filesize(zip_file) == 945849

    tcia_single_image(series = chosen_series, sop = chosen_sop, file = dicom_file)
    @test isfile(dicom_file)
    @test filesize(dicom_file) == 980794
end

@testset "Download series" begin
  series = tcia_series(collection = "AAPM-RT-MAC", patient = "RTMAC-LIVE-001")
  seriesjs = tcia_series(collection = "AAPM-RT-MAC", patient = "RTMAC-LIVE-001", format="json")
  download_series(series, "./testdf")
  download_series(seriesjs, "./testjs")
  download_series(series, "./testdf"; overwrite = false)
end

@testset "Utilities - remove_empty!()" begin
    dict_potentialy_empty_values = Dict(1 => "", 2 => "hello", 3 => "b", 4 => "", 5 => "ye")
    CancerImagingArchive.remove_empty!(dict_potentialy_empty_values)
    nonempty_keys = []
    for (key, value) in dict_potentialy_empty_values
        @test !isempty(value)
        push!(nonempty_keys, key)
    end
    @test sort(nonempty_keys) == [2, 3, 5]
end

@testset "Utilities - Data writer" begin
    tabular_data = tcia_collections()
    dataframe_to_csv(dataframe = tabular_data, file = csv_file)
    @test isfile(csv_file)
    println("Size of csv file: $(filesize(csv_file))")
    @test filesize(csv_file) >= 1346

    dict_array = tcia_collections(format = "json")
    dictionary_to_json(dictionary = dict_array, file = json_file)
    @test isfile(json_file)
    println("Size of json file: $(filesize(json_file))")
    @test filesize(json_file) >= 4816
end


text:
############################################################################################
#Import
include("memory.jl")
include("diagnostics.jl")
include("tune.jl")

############################################################################################
#export
#export


text:
"""
    catmullrom_by_arclength( points )
    catmullrom_by_arclength( points, (min_between_points, max_between_points))

Calculates a Catmull-Rom fit through given `points`, where the number of interpoint knots
varys inversely with the approximate arclength of that arc.

`arcpoints_min` and `arcpoints_max` refer any one arc between two points.
""" catmullrom_by_arclength

function catmullrom_by_arclength(points::P, arcpoints::Tuple{I,I}=(ArcpointsMin, ArcpointsMax); extend::Bool=true) where {P, I<:Integer}
    catmullrom_requirement(npoints(points))

    crpoints = deepcopy(points)
    if extend
        extend_seq(crpoints)
    end

    return catmullrombyarclength(crpoints, minmax(arcpoints...))
end

function catmullrombyarclength(points::P, arcpoints::Tuple{I,I}=(ArcpointsMin, ArcpointsMax)) where {P, I<:Integer}
    n_points = npoints(points)
    pointsperarc = arclength_interpolants(points, arcpoints)
    total_points = sum(pointsperarc) + n_points
    n_coords = ncoords(points)
    T = coordtype(points)
    result = []
    for i=1:n_coords
        push!(result, T[])
    end
    point_seqs = [points[i:i+3] for i=1:(length(points)-3)]

    n_points_3 = n_points - 3
    for idx=1:(n_points-3)
        fourpoints = point_seqs[idx]
        n_between  = pointsperarc[idx]
        fitted = catmullrom_splines(fourpoints, n_between)
        for j=1:n_coords
            append!(result[j], fitted[j][1:(end-1+(idx==n_points_3))])
        end
    end

    return result
end

"""
    arclength_interpolants(points; arcpoints_min=2, arcpoints_max=64)
Convert a sequence of points and a count of additional points to interpolate
to a sequence of arc-length specific point counts where each arc has at least
`arcpoints_min` and at most `arcpoints_max` points.
"""
function arclength_interpolants(points::P, arcpoints::Tuple{I,I}=(ArcpointsMin, ArcpointsMax)) where {P, I<:Integer}
    normalized_arclengths = normalized_catmullrom_arclengths(points)

    arcpoints_min, arcpoints_max = arcpoints
    arclength_min = minimum(normalized_arclengths)
    multiplier = nextfloat(arcpoints_min / arclength_min)
    arcpoints = floor.(Int, normalized_arclengths .* multiplier)
    arcpoints = arcpoints .+ isodd.(arcpoints)

    most_arcpoints = maximum(arcpoints)
    if most_arcpoints > arcpoints_max
        multiplier = prevfloat(arcpoints_max / (most_arcpoints + arcpoints_min))
        arcpoints  = floor.(Int, muladd.(arcpoints .- arcpoints_min, multiplier, arcpoints_min))
        arcpoints = arcpoints .+ isodd.(arcpoints)
    end

    return arcpoints
end

"""
    normalized_catmullrom_arclengths(points)
Convert a sequence of CatmullRom points, those given
as data and those interpolated through the data,
to a corresponding sequence of normalized arclengths.
The sum of the normalized arclengths is 1.0.  Each
arc's approximated length is a portion of unity.
The most probable arclength, and the mean of these
normalized arclengths are roughly equal; given by:
`1//number_of_arcs`.
The way to use normalized arclengths is to multiply
the fractional lengths by the number of "plottable"
points (the given data + the interpoint interpolants).
```
smallest_arc = minimum(normalized_arclengths)
rational_arc = rationalize(Float16(smallest_arc))
scalefactor = points_to_realize / denominator(rational_arc)
# scaled_arclengths_are_point_counts counting points on the arc
points_per_arc = round.(Int, normalized_arclengths .* scalefactor)
```
"""
function normalized_catmullrom_arclengths(points::P) where P
    arclengths = approx_catmullrom_arclengths(points)
    sum_of_arcs = sum(arclengths)
    arclengths[:] = arclengths ./ sum_of_arcs
    return arclengths
end

"""
    approx_catmullrom_arclengths(points)
Convert a sequence of CatmullRom points, those given
as data and those interpolated through the data,
to a corresponding sequence of arclength approximations.
"""
function approx_catmullrom_arclengths(points::P) where P
    L = float(coordtype(P))
    n_points = npoints(points)
    # the first two points and the last two points are boundary + anchor
    # boundary anchor first_internal_point ... final_internal_point anchor boundary
    #       noarc   arc                    arcs                   arc   noarc
    #    pt1     pt2     pt3                            ptN-2       ptN-1     ptN
    #                                   pt4 .. ptN-3
    #               1 arc               ((N-3)-4 arcs)          1 arc
    #  N-3-4+1+1 = N-1-4 = Npts-5 arcs
    n_arcs = n_points - 3
    result = Array{L, 1}(undef, n_arcs)

    # points[1:4], points[2:5], .. points[idx:idx+3] .., points[N-3:N]
    for idx = 1:n_points-3
        arc = approx_catmullrom_arclength(points[idx:idx+3]...,)
        result[idx] = arc
    end
    return result
end

"""
    approx_catmullrom_arclength(pt0, pt1, pt2, pt3)
Given 4 ND points along a centripetal Catmull-Rom span,
roughly approximate the arclength of the curvilinear segment
that would be determined by the two bounding points
and the tangents they determine [the arc between p1 and p2].
This well-behaved approximation was developed by Jens Gravesen
```
    (2*corddist + (n-1)*bezdist)/(n+1), n is degree of the curve
    deg=2 --> (2*corddist + bezdist)/(3)
    deg=3 --> (2*corddist + 2*bezdist)/(4) --> (corddist + bezdist)/2
    deg=4 --> (2*corddsit + 3*bezdist)/(5)
```
"""
function approx_catmullrom_arclength(p0::T, p1::T, p2::T, p3::T) where T
    b0, b1, b2, b3 = catmullrom_as_bezier(p0, p1, p2, p3)
    cordal_dist = norm(b3 .- b0)
    bezier_dist = norm(b1 .- b0) + norm(b2 .- b1) + norm(b3 .- b2)

    return (cordal_dist + bezier_dist) * 0.5
end

function catmullrom_as_bezier(p0::T, p1::T, p2::T, p3::T) where T
    b0 = p1
    b3 = p2

    d1 = norm(p1 .- p0); d1a = sqrt(d1)
    d2 = norm(p2 .- p1); d2a = sqrt(d2)
    d3 = norm(p3 .- p2); d3a = sqrt(d3)

    b1n = @. (d1 * p2) - (d2 * p0) + ((2*d1 + 3*d1a*d2a+d2) * p1)
    b1d = 3*d1a*(d1a+d2a)
    b1 = b1n ./ b1d
    b2n = @. (d3 * p1) - (d2 * p3) + ((2*d3 + 3*d3a*d2a+d2) * p2)
    b2d = 3*d3a*(d3a+d2a)
    b2 = b2n ./ b2d

    return b0, b1, b2, b3
end


text:
"""
    COFFSections

COFF Section header table type, containing information about the number of
sections within the COFF object, the location of the section headers, etc...
"""
struct COFFSections{H <: COFFHandle} <: Sections{COFFHandle}
    handle::H
end

Sections(oh::COFFHandle) = COFFSections(oh)

# Implement Sections API
handle(sections::COFFSections) = sections.handle
endof(sections::COFFSections) = num_sections(header(handle(sections)))

"""
    COFFSection

COFF Section header type, containing information about a `Section` of the COFF
file such as its name, its size, location within memory, etc...
"""
@io struct COFFSection{H <: COFFHandle} <: Section{H}
    Name::fixed_string{UInt64}
    VirtualSize::UInt32
    VirtualAddress::UInt32
    SizeOfRawData::UInt32
    PointerToRawData::UInt32
    PointerToRelocations::UInt32
    PointerToLinenumbers::UInt32
    NumberOfRelocations::UInt16
    NumberOfLinenumbers::UInt16
    Characteristics::UInt32
end

function section_name(s::COFFSection)
    name = unsafe_string(s.Name)
    if name[1] == '/'
        # Wow, COFF files are weird
        return string("strtab@", parse(Int, name[2:end]))
    end
    return name
end

section_size(s::COFFSection) = s.SizeOfRawData
section_address(s::COFFSection) = s.VirtualAddress
section_offset(s::COFFSection) = s.PointerToRawData

"""
    COFFSectionRef

COFF `SectionRef` type, containing a `COFFSection` and important metadata..
"""
struct COFFSectionRef{H <: COFFHandle} <: SectionRef{H}
    sections::COFFSections{H}
    section::COFFSection{H}
    idx::UInt32
end

function SectionRef(sections::COFFSections, s::COFFSection, idx)
    return COFFSectionRef(sections, s, UInt32(idx))
end
deref(s::COFFSectionRef) = s.section
sections(s::COFFSectionRef) = s.sections
handle(s::COFFSectionRef) = handle(sections(s))
section_number(s::COFFSectionRef) = s.idx
function section_name(s::COFFSectionRef)
    name = unsafe_string(deref(s).Name)
    if !isempty(name) && name[1] == '/'
        # Wow, COFF files are weird.
        strtab = StrTab(handle(s))
        return strtab_lookup(strtab, parse(Int, name[2:end]))
    end
    return name
end


text:
# Optimization tests

@testset "Optimization tests" begin
@testset "fcc 7" begin
	fcc7 = MkCell.cellopt(cellfcc, 7, 2)
    # best cell
    bestcell = fcc7[1]
    @test bestcell ≈ [-0.5 -1.0 -0.5; -0.5 0.5 -1.0; 1.0 -0.5 -0.5]
    @test abs(det(bestcell)) ≈ 7.0*abs(det(cellfcc))
    # best matrix
    @test fcc7[2] ≈ [1.0 -2.0 0.0; -1.0 -1.0 2.0; -1.0 0.0 -1.0]
    score_tuple = fcc7[3]
    @test MkCell.score_sphere(bestcell) ≈ 1.118033988749895  # Rmax
    @test MkCell.score_image_distance(bestcell) ≈ 1.224744871391589  # dmin
    @test MkCell.score_angles_cell(bestcell) ≈ 0.08411657770309003  # angle deviation
    @test MkCell.score_lengths_cell(bestcell) ≈ 0.0  # length deviation
end;
@testset "cubic 13" begin
	cubic13 = MkCell.cellopt(cellcubic, 13, 2)
    # best cell
    bestcell = cubic13[1]
    @test bestcell ≈ [2.0 -1.0 -1.0; -1.0 -2.0 -1.0; 0.0 1.0 -2.0]
    @test abs(det(bestcell)) ≈ 13.0*abs(det(cellcubic))
    # best matrix
    @test cubic13[2] ≈ [2.0 -1.0 -1.0; -1.0 -2.0 -1.0; 0.0 1.0 -2.0]
    score_tuple = cubic13[3]
    @test MkCell.score_sphere(bestcell) ≈ 2.29128784747792  # Rmax
    @test MkCell.score_image_distance(bestcell) ≈ 2.23606797749979  # dmin
    @test MkCell.score_angles_cell(bestcell) ≈ 0.06174929182486454  # angle deviation
    @test MkCell.score_lengths_cell(bestcell) ≈ 0.09109769979335518  # length deviation
end;
end;


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ImarisWriter"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JaneliaSciComp/ImarisWriter.git", "da070d3fce8a1da4ca739c9e91054261d2df79a7")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ImarisWriter/
mkdir build && cd build
if [[ "${target}" == *-mingw* ]]; then
    FLAGS="-DHDF5_FALLBACK_LIBRARIES=${libdir}/libhdf5-0.dll"
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    "${FLAGS}" \
    ..
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
#
# Mac build fails due to lack of features.h. See H5_HAVE_FEATURES_H in
# /opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root/usr/local/include/H5pubconf.h
# /opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root/usr/local/include/H5public.h
#
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libbpImarisWriter96", :libbpImarisWriter96)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="HDF5_jll", uuid="0234f1f7-429e-5d53-9886-15a909be8d59"))
    Dependency(PackageSpec(name="Lz4_jll", uuid="5ced341a-0733-55b8-9ab6-a4889d929147"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import MathOptInterface
const MOI = MathOptInterface

mutable struct Optimizer{T} <: MOI.AbstractOptimizer
    problem::SATProblem
    status::Symbol
    results::Vector{Int}
end
const NOT_CALLED = :notcalled
function Optimizer{T}() where {T}
    return Optimizer{T}(SATProblem(0, Vector{Int}[]), NOT_CALLED, Int[])
end
# JuMP will expects `VariablePrimal` to return `Float64`
Optimizer() = Optimizer{Float64}()
function MOI.empty!(model::Optimizer)
    model.problem.num_variables = 0
    empty!(model.problem.clauses)
    model.status = NOT_CALLED
    empty!(model.results)
end
function MOI.is_empty(model::Optimizer)
    return iszero(model.problem.num_variables) && isempty(model.problem.clauses)
end
MOI.supports_add_constrained_variable(::Optimizer, ::Type{MOI.ZeroOne}) = true
function MOI.add_constrained_variable(model::Optimizer, ::MOI.ZeroOne)
    model.problem.num_variables += 1
    id = model.problem.num_variables
    return MOI.VariableIndex(id), MOI.ConstraintIndex{MOI.SingleVariable,MOI.ZeroOne}(id)
end
struct CNF <: MOI.AbstractVectorSet
    not::Vector{Bool}
end
MOI.dimension(set::CNF) = length(set.not)
Base.copy(set::CNF) = CNF(Base.copy(set.not))
_not(id::Int, not::Bool) = not ? -id : id
MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{CNF}) = true
function MOI.add_constraint(model::Optimizer, func::MOI.VectorOfVariables, set::CNF)
    clause = Int[_not(func.variables[i].value, set.not[i]) for i in eachindex(set.not)]
    push!(model.problem, clause)
    return MOI.ConstraintIndex{typeof(func),typeof(set)}(length(model.problem.clauses))
end
function MOI.optimize!(model::Optimizer)
    model.status, model.results = solve(model.problem)
end
const TERM = Dict(NOT_CALLED => MOI.OPTIMIZE_NOT_CALLED, :sat => MOI.OPTIMAL, :unsat => MOI.INFEASIBLE, :unknown => MOI.TIME_LIMIT)
MOI.get(model::Optimizer, ::MOI.TerminationStatus) = TERM[model.status]
MOI.get(model::Optimizer, ::MOI.ResultCount) = model.status == :sat ? 1 : 0
function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if 1 <= attr.N <= MOI.get(model, MOI.ResultCount())
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end
MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION
MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex) = model.results[vi.value] > 0


text:
# Autogenerated wrapper script for ruby_jll for x86_64-linux-gnu
export RUBYLIB, bundle, bundler, erb, gem, irb, libruby, racc, racc2y, rake, rdoc, ri, ruby

using Libiconv_jll
using OpenSSL_jll
using Readline_jll
using Zlib_jll
using Gdbm_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `RUBYLIB`
const RUBYLIB_splitpath = ["RUBYLIB"]

# This will be filled out by __init__() for all products, as it must be done at runtime
RUBYLIB_path = ""

# RUBYLIB-specific global declaration
# This will be filled out by __init__()
RUBYLIB = ""

# Relative path to `bundle`
const bundle_splitpath = ["bin", "bundle"]

# This will be filled out by __init__() for all products, as it must be done at runtime
bundle_path = ""

# bundle-specific global declaration
function bundle(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(bundle_path)
    end
end

# Relative path to `bundler`
const bundler_splitpath = ["bin", "bundler"]

# This will be filled out by __init__() for all products, as it must be done at runtime
bundler_path = ""

# bundler-specific global declaration
function bundler(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(bundler_path)
    end
end

# Relative path to `erb`
const erb_splitpath = ["bin", "erb"]

# This will be filled out by __init__() for all products, as it must be done at runtime
erb_path = ""

# erb-specific global declaration
function erb(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(erb_path)
    end
end

# Relative path to `gem`
const gem_splitpath = ["bin", "gem"]

# This will be filled out by __init__() for all products, as it must be done at runtime
gem_path = ""

# gem-specific global declaration
function gem(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(gem_path)
    end
end

# Relative path to `irb`
const irb_splitpath = ["bin", "irb"]

# This will be filled out by __init__() for all products, as it must be done at runtime
irb_path = ""

# irb-specific global declaration
function irb(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(irb_path)
    end
end

# Relative path to `libruby`
const libruby_splitpath = ["lib", "libruby.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libruby_path = ""

# libruby-specific global declaration
# This will be filled out by __init__()
libruby_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libruby = "libruby.so.2.7"

# Relative path to `racc`
const racc_splitpath = ["bin", "racc"]

# This will be filled out by __init__() for all products, as it must be done at runtime
racc_path = ""

# racc-specific global declaration
function racc(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(racc_path)
    end
end

# Relative path to `racc2y`
const racc2y_splitpath = ["bin", "racc2y"]

# This will be filled out by __init__() for all products, as it must be done at runtime
racc2y_path = ""

# racc2y-specific global declaration
function racc2y(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(racc2y_path)
    end
end

# Relative path to `rake`
const rake_splitpath = ["bin", "rake"]

# This will be filled out by __init__() for all products, as it must be done at runtime
rake_path = ""

# rake-specific global declaration
function rake(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(rake_path)
    end
end

# Relative path to `rdoc`
const rdoc_splitpath = ["bin", "rdoc"]

# This will be filled out by __init__() for all products, as it must be done at runtime
rdoc_path = ""

# rdoc-specific global declaration
function rdoc(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(rdoc_path)
    end
end

# Relative path to `ri`
const ri_splitpath = ["bin", "ri"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ri_path = ""

# ri-specific global declaration
function ri(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ri_path)
    end
end

# Relative path to `ruby`
const ruby_splitpath = ["bin", "ruby"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ruby_path = ""

# ruby-specific global declaration
function ruby(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ruby_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"ruby")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Libiconv_jll.PATH_list, OpenSSL_jll.PATH_list, Readline_jll.PATH_list, Zlib_jll.PATH_list, Gdbm_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Libiconv_jll.LIBPATH_list, OpenSSL_jll.LIBPATH_list, Readline_jll.LIBPATH_list, Zlib_jll.LIBPATH_list, Gdbm_jll.LIBPATH_list,))

    global RUBYLIB_path = normpath(joinpath(artifact_dir, RUBYLIB_splitpath...))

    global RUBYLIB = RUBYLIB_path
    global bundle_path = normpath(joinpath(artifact_dir, bundle_splitpath...))

    push!(PATH_list, dirname(bundle_path))
    global bundler_path = normpath(joinpath(artifact_dir, bundler_splitpath...))

    push!(PATH_list, dirname(bundler_path))
    global erb_path = normpath(joinpath(artifact_dir, erb_splitpath...))

    push!(PATH_list, dirname(erb_path))
    global gem_path = normpath(joinpath(artifact_dir, gem_splitpath...))

    push!(PATH_list, dirname(gem_path))
    global irb_path = normpath(joinpath(artifact_dir, irb_splitpath...))

    push!(PATH_list, dirname(irb_path))
    global libruby_path = normpath(joinpath(artifact_dir, libruby_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libruby_handle = dlopen(libruby_path)
    push!(LIBPATH_list, dirname(libruby_path))

    global racc_path = normpath(joinpath(artifact_dir, racc_splitpath...))

    push!(PATH_list, dirname(racc_path))
    global racc2y_path = normpath(joinpath(artifact_dir, racc2y_splitpath...))

    push!(PATH_list, dirname(racc2y_path))
    global rake_path = normpath(joinpath(artifact_dir, rake_splitpath...))

    push!(PATH_list, dirname(rake_path))
    global rdoc_path = normpath(joinpath(artifact_dir, rdoc_splitpath...))

    push!(PATH_list, dirname(rdoc_path))
    global ri_path = normpath(joinpath(artifact_dir, ri_splitpath...))

    push!(PATH_list, dirname(ri_path))
    global ruby_path = normpath(joinpath(artifact_dir, ruby_splitpath...))

    push!(PATH_list, dirname(ruby_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
module Tomoko
##
using Random
using Distributions
using SpecialFunctions
using JuliennedArrays
using DataFrames
using Optim
using StatsBase
using HypergeometricFunctions

include("pop_state.jl")
include("run_sim.jl")
include("pop_stats.jl")
include("env_var.jl")
include("df_analysis.jl")
include("escape_prob.jl")
include("sensitivity_analysis.jl")
include("WrightSampler.jl")

end # module

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
using GCoptimization

width = 10
height = 5
pixelNum = width * height
labelNum = 7
result = zeros(pixelNum)

gco = GCoptimizationGridGraph(width, height, labelNum)

function dataFn(pixel::GCOSiteID, label::GCOLabelID)::GCOEnergyTermType
    if pixel < 25
        return label == 0 ? 0 : 10
    else
        return label == 5 ? 0 : 10
    end
end
setDataCost(gco, dataFn)

smoothFn(p1::GCOSiteID, p2::GCOSiteID, l1::GCOLabelID, l2::GCOLabelID)::GCOEnergyTermType = (l1-l2)*(l1-l2) <= 4 ? (l1-l2)*(l1-l2) : 4
setSmoothCost(gco, smoothFn)

println("Before optimization energy is ", compute_energy(gco))
expansion(gco, 2)
println("After optimization energy is ", compute_energy(gco))

result = [whatLabel(gco, i-1) for i = 1:pixelNum]

showall(result)


text:
@polly function kernel_heat_3d(tsteps, A, B)
    n = size(A,1)

    for t = 1:tsteps
        for i = 2:(n-1), j = 2:(n-1), k = 2:(n-1)
            B[i,j,k] =   0.125 * (A[i+1,j,k] - 2.0 * A[i,j,k] + A[i-1,j,k]) +
                         0.125 * (A[i,j+1,k] - 2.0 * A[i,j,k] + A[i,j-1,k]) +
                         0.125 * (A[i,j,k+1] - 2.0 * A[i,j,k] + A[i,j,k-1]) +
                         A[i,j,k]
        end
        for i = 2:(n-1), j = 2:(n-1), k = 2:(n-1)
            A[i,j,k] =   0.125 * (B[i+1,j,k] - 2.0 * B[i,j,k] + B[i-1,j,k]) +
                         0.125 * (B[i,j+1,k] - 2.0 * B[i,j,k] + B[i,j-1,k]) +
                         0.125 * (B[i,j,k+1] - 2.0 * B[i,j,k] + B[i,j,k-1]) +
                         B[i,j,k]
        end
    end
end

let
    tsteps = 500
    n = 120

    A = zeros(Float32, n, n, n)
    B = zeros(Float32, n, n, n)

    for i = 1:n, j = 1:n, k = 1:n
        A[i,j,k] = B[i,j,k] = (i + j + (n-k))* 10 / (n)
    end

    SUITE["heat-3d"] = @benchmarkable kernel_heat_3d($tsteps, A, B) setup = (A = copy($A); B = copy($B))
end