text:
module AircraftGeometry

## Package imports
#==========================================================================================#

using Base.Math
using Base.Iterators
using DelimitedFiles
using StaticArrays
using CoordinateTransformations
using Rotations
using LinearAlgebra
using SplitApplyCombine
using Interpolations

# Math tools
import ..MathTools: uniform_spacing, linear_spacing, sine_spacing, cosine_spacing, cosine_interp, splitat, adj3, slope, columns, forward_sum, forward_division, forward_difference, weighted_vector, vectarray, extend_yz

# Panel geometry
import ..PanelGeometry: Panel2D, Panel3D, panel_area, panel_normal, transform, make_panels

import ..AeroMDAO: properties

## Types
#==========================================================================================#

abstract type AbstractAircraft end

abstract type AbstractWing <: AbstractAircraft end

## Foil geometry
#==========================================================================================#

include("Foils/foil.jl")
include("Foils/class_shape_transformation.jl")
include("Foils/naca_airfoils.jl")

## Fuselage geometry
#==========================================================================================#

include("fuselage.jl")

## Wing geometry
#==========================================================================================#

aspect_ratio(span, area) = span^2 / area
mean_geometric_chord(span, area) = area / span
taper_ratio(root_chord, tip_chord) = tip_chord / root_chord
area(span, chord) = span * chord
mean_aerodynamic_chord(root_chord, taper_ratio) = (2/3) * root_chord * (1 + taper_ratio + taper_ratio^2)/(1 + taper_ratio)
y_mac(y, b, λ) = y + b / 2 * (1 + 2λ) / 3(1 + λ)
quarter_chord(chord) = 0.25 * chord

include("Wings/halfwing.jl")
include("Wings/wing.jl")
include("Wings/mesh_tools.jl")
include("Wings/mesh_wing.jl")
include("Wings/controls.jl")

"""
    span(wing :: AbstractWing)

Compute the planform span of an `AbstractWing`.
"""
span(wing :: AbstractWing) = span(wing)

"""
    projected_area(wing :: AbstractWing)

Compute the projected area of an `AbstractWing`` by summing the trapezoidal areas.
"""
projected_area(wing :: AbstractWing) = projected_area(wing)

"""
    mean_aerodynamic_chord(wing :: AbstractWing)

Compute the mean aerodynamic chord of an `AbstractWing`.
"""
mean_aerodynamic_chord(wing :: AbstractWing) = mean_aerodynamic_chord(wing)

"""
    mean_aerodynamic_chord(wing :: AbstractWing)

Compute the coordinates of the mean aerodynamic center of an `AbstractWing`.
"""
mean_aerodynamic_center(wing :: AbstractWing) = mean_aerodynamic_center(wing)

"""
    aspect_ratio(wing :: AbstractWing)

Compute the aspect ratio of an `AbstractWing`.
"""
aspect_ratio(wing) = aspect_ratio(span(wing), projected_area(wing))

"""
    taper_ratio(wing :: AbstractWing)

Compute the taper ratio of an `AbstractWing`, defined as the tip chord length divided by the root chord length.
"""
taper_ratio(wing :: AbstractWing) = taper_ratio(wing)

"""
    properties(wing :: AbstractWing)

Compute the generic properties of interest (span, area, etc.) of an `AbstractWing`.
"""
properties(wing :: AbstractWing) = [ aspect_ratio(wing), span(wing), projected_area(wing), mean_aerodynamic_chord(wing), mean_aerodynamic_center(wing) ]

function Base.show(io :: IO, wing :: AbstractWing)
    println(io, supertype(typeof(wing)),  " with ", length(spans(wing)), " spanwise section(s).")
    println(io, "Aspect Ratio: ", aspect_ratio(wing))
    println(io, "Span (m): ", span(wing))
    println(io, "Projected Area (m): ", projected_area(wing))
    println(io, "Mean Aerodynamic Chord (m): ", mean_aerodynamic_chord(wing))
    println(io, "Mean Aerodynamic Center (m): ", mean_aerodynamic_center(wing))

    nothing
end

end


text:
function vec2param!(Θvec,Θparam::TrackedArray)
    Θparam.data .= reshape(Θvec, size(Θparam))
    return Θparam
end
function vec2param!(Θvec,Θparam::AbstractArray)
    Θparam .= reshape(Θvec, size(Θparam))
    return Θparam
end

function vec2param!(Θvec,Θparm::Tuple)
    cnt = 0
    for k=1:length(Θparm)
        nk = lengthvec(Θparm[k])
        vec2param!(Θvec[cnt+1:cnt+nk], Θparm[k])
        cnt+=nk
    end
    return Θparm
end

lengthvec(Θparm::AbstractArray) = length(Θparm)

function lengthvec(Θparm::Tuple)
    cnt = 0;
    for k=1:length(Θparm)
        cnt+= lengthvec(Θparm[k])
    end
    return cnt
end

getParmsType(Θ::TrackedArray) = typeof(Θ[1].data)
getParmsType(Θ::AbstractArray) = typeof(Θ[1])
getParmsType(Θ::Tuple)        = getParmsType(Θ[1])

function param2vec(Θparm::Tuple)
    cnt = lengthvec(Θparm)
    R = getParmsType(Θparm)
    Θvec = zeros(R, cnt)
    return param2vec!(Θparm,Θvec)
end

param2vec(Θparm::TrackedArray) = vec(Θparm.data)
param2vec(Θparm::AbstractArray) = vec(Θparm)
function param2vec!(Θparm,Θvec)
    cnt = 0
    for k=1:length(Θparm)
        nk = lengthvec(Θparm[k])
        Θvec[cnt+1:cnt+nk] .= param2vec(Θparm[k])
        cnt+=nk
    end
    return Θvec
end

ext: jl


text:
function load(path_to_file::String)::DataFrame
    return CSV.read(path_to_file, DataFrame)
end

function load_tga_training_data()

    # load the data -
    path_data_file_TF = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF.csv")
    full_data_table_TF = load(path_data_file_TF)
    path_data_file_TF_TM = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF-TM.csv")
    full_data_table_TF_TM = load(path_data_file_TF_TM)

    # return a complete set -
    return vcat(full_data_table_TF, full_data_table_TF_TM)
end

function load_fibrinolysis_training_data()
    # load the data -
    path_data_file_TF = joinpath(_PATH_TO_DATA, "Training-Fibrinolysis-TF.csv")
    full_data_table_TF = load(path_data_file_TF)
    return full_data_table_TF
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module Graficos

using Plots
import ..FEA: Estrutura, forcas_internas
import ..Otm: Otimizador

export plotar_estrutura_inicial,
       plotar_estrutura_otimizada

"""
    maior_dimensao(estrutura::Estrutura)::Float64

Retorna o tamanho da maior dimensão (horizontal ou vertical) da estrutura.
"""
function maior_dimensao(estrutura::Estrutura)::Float64
    min_x = Inf
    max_x = -Inf
    min_y = Inf
    max_y = -Inf

    for no in estrutura.nos
        min_x = min(min_x, no.x)
        max_x = max(max_x, no.x)
        min_y = min(min_y, no.y)
        max_y = max(max_y, no.y)
    end

    return max(max_x - min_x, max_y - min_y)
end

"""
    plotar_forcas_apoios!(p, estrutura::Estrutura, estrutura_inicial=true, ponderador=0.05)

Plota no gráfico `p` as forças e os apoios da estrutura.

# Argumentos
- `p`: Gráfico (plot).
- `estrutura::Estrutura`: Estrutura analisada.
- `estrutura_inicial::Bool=true`: Parâmetro que altera as cores dos indicadores de forças e apoios
                                  caso a estrutura seja a inicial ou a otimizada.
- `ponderador::Real`: Parâmetro que indica o tamanho dos indicadores de forças e apoios. É proporcional
                      à maior dimensão da estrutura.
"""
function plotar_forcas_apoios!(p, estrutura::Estrutura, estrutura_inicial::Bool=true, ponderador::Real=0.05)
    lmax = maior_dimensao(estrutura)
    cor_forcas = estrutura_inicial ? :red : :green
    cor_apoios = estrutura_inicial ? :blue : :black
    # Plotagem dos apoios e das forças.
    for no in estrutura.nos
        # Apoios
        if no.rx
            x = [NaN, no.x - ponderador/2*lmax, no.x]
            y = [NaN, no.y, no.y]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_apoios)
        end
        if no.ry
            x = [NaN, no.x, no.x]
            y = [NaN, no.y - ponderador/2*lmax, no.y]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_apoios)
        end

        # Forças
        if no.fx > 0
            x = [NaN, no.x, no.x + ponderador*lmax]
            y = [NaN, no.y, no.y]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_forcas, arrow=true, arrowsize=ponderador*10)
        elseif no.fx < 0
            x = [NaN, no.x, no.x - ponderador*lmax]
            y = [NaN, no.y, no.y]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_forcas, arrow=true, arrowsize=ponderador*10)
        end
        if no.fy > 0
            x = [NaN, no.x, no.x]
            y = [NaN, no.y, no.y + ponderador*lmax]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_forcas, arrow=true, arrowsize=ponderador*10)
        elseif no.fy < 0
            x = [NaN, no.x, no.x]
            y = [NaN, no.y, no.y - ponderador*lmax]
            plot!(p, x, y, lw=ponderador*50, label="", color=cor_forcas, arrow=true, arrowsize=ponderador*10)
        end
    end
end

"""
    plotar_estrutura_inicial(estrutura::Estrutura)

Plota a estrutura antes da otimização.
"""
function plotar_estrutura_inicial(estrutura::Estrutura)
    p = plot()
    # Plotagem dos elementos.
    x = []
    y = []
    for el in estrutura.elementos
        push!(x, [NaN, el.no₁.x, el.no₂.x]...)
        push!(y, [NaN, el.no₁.y, el.no₂.y]...)
    end
    plot!(p, x, y, label="", color=:black, aspect_ratio=:equal)

    plotar_forcas_apoios!(p, estrutura, true)

    display(p)
    sleep(10000)
end

"""
    plotar_estrutura_otimizada(otimizador::Otimizador)

Plota a estrutura após a otimização.
"""
function plotar_estrutura_otimizada(otimizador::Otimizador)
    x_norm = otimizador.x / maximum(otimizador.x)
    forcas = forcas_internas(otimizador.estrutura)

    p = plot()
    for i in eachindex(otimizador.estrutura.elementos)
        area = x_norm[i]
        if isapprox(area, 0; atol=eps())
            continue
        end

        forca = forcas[i]
        el = otimizador.estrutura.elementos[i]

        x = [NaN, el.no₁.x, el.no₂.x]
        y = [NaN, el.no₁.y, el.no₂.y]

        plot!(p, x, y, label="", lw=area*5, color= forca ≤ 0 ? :blue : :red, aspect_ratio=:equal)
    end

    plotar_forcas_apoios!(p, otimizador.estrutura, false)

    display(p)
    sleep(10000)
end

end

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
# Use Julia's `jl_type_morespecific` function to emulate Julia's multiple dispatch across
# generic functions.
#
# Origin:
# https://github.com/JuliaLang/julia/blob/master/doc/src/devdocs/types.md#subtyping-and-method-sorting
type_morespecific(a, b) = ccall(:jl_type_morespecific, Bool, (Any, Any), a, b)

"""
    anonymous_signature(m::Method) -> Type{<:Tuple}

Construct a Tuple of the methods signature with the function type removed

# Example
```jldoctest
julia> m = first(methods(first, (String,)));

julia> m.sig
Tuple{typeof(first),Any}

julia> anonymous_signature(m)
Tuple{Any}
```
"""
anonymous_signature(m::Method) = anonymous_signature(m.sig)
anonymous_signature(sig::DataType) = Tuple{sig.parameters[2:end]...}
anonymous_signature(sig::UnionAll) = UnionAll(sig.var, anonymous_signature(sig.body))

"""
    anon_morespecific(a::Method, b::Method) -> Bool

Determine which method is more specific for multiple dispatch without considering the
function type. By not considering the function type we can determine which method is more
specific as if they are a part of the same generic function.
"""
function anon_morespecific(a::Method, b::Method)
    # Drop the function type from the parameter
    a_sig = anonymous_signature(a)
    b_sig = anonymous_signature(b)

    return type_morespecific(a_sig, b_sig)
end

"""
    dispatch(funcs::AbstractVector, args...) -> Tuple{Method, Any}

Choose which method to execute based upon the provided arguments (values not types).
Emulates Julia's multiple dispatch system but allows for dispatching between methods of
multiple generic functions instead of just methods of a single generic function. Returns a
tuple of the selected method and the generic function of the method.

When the function to dispatch to is ambiguous last ambiguous function in the vector is used.
"""
function dispatch(funcs::AbstractVector, args...)
    arg_types = map(Core.Typeof, args)

    best_method = nothing
    best_function = nothing
    for f in reverse(funcs)
        # Since arguments will be using concrete types `methods` should only return up to
        # one method.
        for m in methods(f, arg_types)
            if best_method === nothing || anon_morespecific(m, best_method)
                best_method = m
                best_function = f
            end
        end
    end

    return best_method, best_function
end


text:
# ---
# title: 1359. Count All Valid Pickup and Delivery Options
# id: problem1359
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/count-all-valid-pickup-and-delivery-options/description/>
# hidden: true
# ---
#
# Given `n` orders, each order consist in pickup and delivery services.
#
# Count all valid pickup/delivery possible sequences such that delivery(i) is
# always after of pickup(i).
#
# Since the answer may be too large, return it modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: 1
#     Explanation: Unique order (P1, D1), Delivery 1 always is after of Pickup 1.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: 6
#     Explanation: All possible orders:
#     (P1,P2,D1,D2), (P1,P2,D2,D1), (P1,D1,P2,D2), (P2,P1,D1,D2), (P2,P1,D2,D1) and (P2,D2,P1,D1).
#     This is an invalid order (P1,D2,P2,D1) because Pickup 2 is after of Delivery 2.
#
#
# **Example 3:**
#
#
#
#     Input: n = 3
#     Output: 90
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 500`
#
#
## @lc code=start
using LeetCode

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


text:
include("stop_conditions.jl")
include("hooks.jl")
include("core.jl")
include("iteration_methods.jl")


text:
#-------PHREDS AND PROBS--------
#copypasta from rifraf

const Phred = Int8
const Prob = Float64
const LogProb = Float64

const MIN_PHRED = Phred(1)
const MAX_PHRED = Phred(Int('~') - 33)

"""
    phred_to_log_p(x)

Conversion from phred value to log probability value.
"""
@generated function phred_to_log_p(x)
    return quote
        return x / (-10.0)
    end
end

"""
    phred_to_p(q::Phred)

Conversion from phred value to real probability value.
"""
function phred_to_p(q::Phred)
    return exp10(phred_to_log_p(q))
end

"""
    phred_to_p(x::Vector{Phred})

Conversion from phred values to real probability values in vector.
"""
function phred_to_p(x::Vector{Phred})
    return exp10.(phred_to_log_p(x))
end

"""
    p_to_phred(p::Prob)

Conversion from real probability value to phred value.
"""
function p_to_phred(p::Prob)
    return Phred(min(round(-10.0 * log10(p)), MAX_PHRED))
end

"""
    p_to_phred(p::Prob)

Conversion from real probability values to phred values in vector.
"""
function p_to_phred(x::Vector{LogProb})
    return Phred[p_to_phred(p) for p in x]
end

"""
    error_probs_to_phreds(eps::Vector{Float64})

Conversion from error probability values to phred values in vector.
"""
function error_probs_to_phreds(eps::Vector{Float64})
    return [Phred(min(round(-10.0 * log10(p)),99)) for p in eps]
end

#-------ALTERNATE FILTER FUNCTION----------

"""
    quality_filter(infile, outfile=join(split(infile, ".")[1:end-1], ".") * ".filt.fasta";
                   errorRate = 0.01, minLength = 100, labelPrefix = "seq", errorOut = true)

Writes file with sequences from input file that have all sites within error rate margin.
"""
function quality_filter(infile, outfile=join(split(infile, ".")[1:end-1], ".") * ".filt.fasta" ;
                        errorRate = 0.01, minLength = 100, labelPrefix = "seq", errorOut = true)
    seqs, scores, names = read_fastq(infile)

	p_vals = [phred_to_p(score) for score in scores]

    inds = quality_filter_inds(seqs, scores, errorRate=errorRate, minLength=minLength)

    #names = names[inds]
    if errorOut == true
        names = ["$labelPrefix$(i)|ee=$(mean(p_vals[i]))" for i in inds]
    else
        names = ["$labelPrefix$(i)" for i in 1:length(inds)]
    end

    if last(outfile) == 'a'
        write_fasta(outfile, seqs[inds], names=names)
    else
        write_fastq(outfile, seqs[inds], scores[inds], names=names)
    end
end

"""
    quality_filter(seqs::Array{String, 1}, scores::Array, names::Array{String, 1}; errorRate=0.01, minLength=0)

Returns sub arrays with sequences that have site-wise error all within given margin.
"""
function quality_filter(seqs::Array{String, 1}, scores::Array, names::Array{String, 1}; errorRate=0.01,minLength=0)
    inds = quality_filter_inds(seqs, scores, errorRate=errorRate, minLength=minLength)
    return seqs[inds], scores[inds], names[inds]
end

"""
    quality_filter_inds(seqs::Array{String, 1}, scores::Array; errorRate=0.01,minLength=100)

Returns indices of sequences with site-wise error within given margin.
"""
function quality_filter_inds(seqs::Array{String, 1}, scores::Array; errorRate=0.01,minLength=100)
    p_vals = [phred_to_p(score) for score in scores]
    inds = filter!(x -> (mean(p_vals[x]) < errorRate && (length(seqs[x]) >= minLength)) , collect(1:length(seqs)))
    return inds
end

# temporary: for compatibility
seq_filter = quality_filter

#--------QUALITY INTERROGATING FUNCTIONS, GIVEN EITHER EE NAME TAGGING OR FASTQ FORMAT -------

"""
    length_vs_qual(fasta_path; plot_title = "Length Vs Errors")

Creates plot of lengths of sequences vs. mean error rates of sequences.
"""
function length_vs_qual(fasta_path; plot_title = "Length Vs Errors", alpha=0.3)

	if last(fasta_path) == 'a'
		names, seqs = read_fasta_with_names(fasta_path)
    	error_rates = [parse(Float64, split(i,['|', '='])[3]) for i in names]
	else
		seqs, errs, names = read_fastq(fasta_path)
		error_rates = [mean(phred_to_p(e)) for e in errs]
    end

	# To remove that annoying rounding that happens in usearch
    # dithered_rates = [(rand()*10-5)*((i%1)==0.0&&i>90.0)+i for i in error_rates]
	# Following line commented since usearch is no longer used
	# dithered_rates = [(1+rand()*0.1-0.05)*i for i in error_rates]
    string_lengths = [length(i) for i in seqs]
    plot(string_lengths, error_rates, ".", alpha=alpha)
    title(plot_title)
    xlabel("Length")
    ylabel("Predicted Error Rate")
    # return dithered_rates
end

"""
    qual_hist(fasta_path; plot_title = "Quality Hist")

Creates histogram of mean quality scores per sequence.
"""
function qual_hist(fasta_path; plot_title = "Quality Hist")
	if last(fasta_path) == 'a'
		names, seqs = read_fasta_with_names(fasta_path)
    	error_rates = [parse(Float64, split(i,['|', '='])[3]) for i in names]
	else
		seqs, errs, names = read_fastq(fasta_path)
		error_rates = [mean(phred_to_p(e)) for e in errs]
    end
    # to remove that annoying rounding that happens in usearch
    # Following line commented since usearch is no longer used
	#dithered_rates = [(rand()*10-5)*((i%1)==0.0&&i>90.0)+i for i in error_rates]
    string_lengths = [length(i) for i in seqs]
    plt[:hist](error_rates, 40)
    title(plot_title)
    xlabel("Error Rate")
    ylabel("Frequency")
    # return dithered_rates
end


text:
module Layered

using StaticArrays

using Reexport
@reexport using Colors
using Match
import Cairo

export Transform, Layer, Shape, upward_transform, solve!, →, Angle, rad, deg, needed_attributes

include("types.jl")
include("angles.jl")
include("geometry.jl")
include("gradient.jl")
include("attributes.jl")
include("transform.jl")
include("layer.jl")
include("shape.jl")
include("canvas.jl")
include("drawing.jl")
include("video.jl")
include("helpers.jl")

Base.Broadcast.broadcastable(g::GeometricObject) = Ref(g)

function add_shape!(shape, layer::Layer, args...)

end

for geom in geoms

    lowerc = Symbol(lowercase(String(geom)))
    lowerc_plural = Symbol(lowercase(String(geom)) * "s")
    lowerc_mutating = Symbol(lowercase(String(geom)) * "!")
    lowerc_plural_mutating = Symbol(lowercase(String(geom)) * "s!")

    @eval begin
        $geom(g::$geom) = $geom((getfield(g, k) for k ∈ fieldnames($geom))...)

        """
        $($lowerc)(args...)

        Creates a shape containing a `GeometricObject` of type `$($geom)`,
        passing any trailing arguments to the constructor `$($geom)()`."""
        $lowerc(args...) = Shape($geom, args...)

        """
        $($lowerc)(f::Function, args...)

        Creates a shape containing a `GeometricObject` of type `$($geom)`,
        storing any trailing arguments as dependencies, later to be passed as arguments
        to the given `Function` `f` that should return a `GeometricObject` of type `$($geom)`
        and will be evaluated when `solve!` is called on the shape during
        the drawing process."""
        $lowerc(f::Function, args...) = Shape(f, $geom, args...)

        """
        $($lowerc_mutating)(layer::Layer, args...)

        Creates a shape containing a `GeometricObject` of type `$($geom)`,
        passing any trailing arguments to the constructor `$($geom)()`. Then
        appends the created shape to the given `Layer` `layer`."""
        function $lowerc_mutating(layer::Layer, args...)
            x = $lowerc(args...)
            push!(layer, x)
            x
        end

        """
        $($lowerc_mutating)(f::Function, layer::Layer, args...)

        Creates a shape containing a `GeometricObject` of type `$($geom)`,
        storing any trailing arguments as dependencies, later to be passed as arguments
        to the given `Function` `f` that should return a `GeometricObject` of type `$($geom)`
        and will be evaluated when `solve!` is called on the shape during
        the drawing process. This function then appends the created shape to the given `Layer` `layer`."""
        function $lowerc_mutating(f::Function, layer::Layer, args...)
            x = $lowerc(f, args...)
            push!(layer, x)
            x
        end

        """
        $($lowerc_plural)(args...)

        Creates `Shapes` containing `GeometricObject`s of type `$($geom)`,
        passing any trailing arguments to the constructor `$($geom).()`."""
        $lowerc_plural(args...) = Shapes($geom, args...)

        """
        $($lowerc_plural)(f::Function, args...)

        Creates shapes containing `GeometricObject`s of type `$($geom)`,
        storing any trailing arguments as dependencies, later to be passed as arguments
        to the given `Function` `f` that should return `GeometricObject`s of type `$($geom)`
        and will be evaluated when `solve!` is called on the shapes during
        the drawing process."""
        $lowerc_plural(f::Function, args...) = Shapes(f, $geom, args...)

        """
        $($lowerc_plural_mutating)(layer::Layer, args...)

        Creates `Shapes` containing `GeometricObject`s of type `$($geom)`,
        passing any trailing arguments to the constructor `$($geom).()`. Then
        appends the created `Shapes` to the given `Layer` `layer`."""
        function $lowerc_plural_mutating(layer::Layer, args...)
            x = $lowerc_plural(args...)
            push!(layer, x)
            x
        end

        """
        $($lowerc_plural_mutating)(f::Function, layer::Layer, args...)

        Creates multiple shapes containing `GeometricObject`s of type `$($geom)`,
        storing any trailing arguments as dependencies, later to be passed as arguments
        to the given `Function` `f` that should return the `GeometricObject`s of type `$($geom)`
        and will be evaluated when `solve!` is called on the shapes during
        the drawing process. This function then appends the created shapes to the given `Layer` `layer`."""
        function $lowerc_plural_mutating(f::Function, layer::Layer, args...)
            x = $lowerc_plural(f, args...)
            push!(layer, x)
            x
        end

        export $geom, $lowerc, $lowerc_mutating, $lowerc_plural, $lowerc_plural_mutating

    end
end

end # module


text:
TrainingPages = [
    "System" => [
   # "index" ="TrainingPages.md",
   "2008" => [
      "Training/HarmonieSystemTraining2008.md",
      "Training/HarmonieSystemTraining2008/Lecture/CommonMonitoring.md",
      "Training/HarmonieSystemTraining2008/Lecture/DAdataflow.md",
      "Training/HarmonieSystemTraining2008/Lecture/DomainAndCoupling.md",
      "Training/HarmonieSystemTraining2008/Lecture/Forecast.md",
      "Training/HarmonieSystemTraining2008/Lecture/General.md",
      "Training/HarmonieSystemTraining2008/Lecture/Installation.md",
      "Training/HarmonieSystemTraining2008/Lecture/JbStatistics.md",
      "Training/HarmonieSystemTraining2008/Lecture/JobControl.md",
      "Training/HarmonieSystemTraining2008/Lecture/LocalInstallation.md",
      "Training/HarmonieSystemTraining2008/Lecture/ObsHandling.md",
      "Training/HarmonieSystemTraining2008/Lecture/ObsMonitoringDiagnosis.md",
      "Training/HarmonieSystemTraining2008/Lecture/Opening.md",
      "Training/HarmonieSystemTraining2008/Lecture/PostppVerification.md",
      "Training/HarmonieSystemTraining2008/Lecture/SourceCode.md",
      "Training/HarmonieSystemTraining2008/Lecture/SurfaceAssimilation.md",
      "Training/HarmonieSystemTraining2008/Lecture/WikiAuthoring.md",
      "Training/HarmonieSystemTraining2008/Tasklist.md",
      "Training/HarmonieSystemTraining2008/Trainin/DataAssimilation.md",
      "Training/HarmonieSystemTraining2008/Training/DataAssimilation.md",
      "Training/HarmonieSystemTraining2008/Training/DomainAndCoupling.md",
      "Training/HarmonieSystemTraining2008/Training/Forecast.md",
      "Training/HarmonieSystemTraining2008/Training/Installation.md",
      "Training/HarmonieSystemTraining2008/Training/JobControl.md",
      "Training/HarmonieSystemTraining2008/Training/LocalInstallation.md",
      "Training/HarmonieSystemTraining2008/Training/PostppVerification.md",
      "Training/HarmonieSystemTraining2008/Training/SourceCode.md",
      "Training/HarmonieSystemTraining2008/Training/WikiAuthoring.md"
    ],
    "2011" => [
       "Training/HarmonieSystemTraining2011.md",
       "Training/HarmonieSystemTraining2011/Agenda.md",
       "Training/HarmonieSystemTraining2011/Exercises.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/Archiving.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/Building.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/ConfFiles.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/ConfFilesExp.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/ConstFiles.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/GridDepFiles.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/InputFiles.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/LocalChanges.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/LogFiles.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/ModelOutput.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/OnError.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/PreqHard.md",
       "Training/HarmonieSystemTraining2011/Lecture/Installation/PreqSoft.md"
    ],

    "2012" => "Training/HarmonieSystemTraining2012.md",
    "2014" => "Training/HarmonieSystemTraining2014.md",
    "2015" => "Training/HarmonieSystemTraining2015.md",
    "2016" => [
        "Training/HarmonieSystemTraining2016/EPyGrAM.md",
        "Training/HarmonieSystemTraining2016/Exercises.md"
    ]
   ]
]


text:
# Quick example

using BinaryTraits
using BinaryTraits.Prefix: Can

# Let's define two ability traits
abstract type Ability end
@trait Fly as Ability
@trait Swim as Ability

# Define interface contracts for the type.
@implement Can{Fly} by liftoff(_)
@implement Can{Fly} by fly(_, direction::Float64, altitude::Float64)
@implement Can{Fly} by speed(_)::Float64

# Define a data type and assign it traits.
struct Crane end
@assign Crane with Can{Fly},Can{Swim}

# Check compliance.
@check(Crane)
#=
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ liftoff(🔹)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits.jl/src/interface.jl:59
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ fly(🔹, ::Float64, ::Float64)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits.jl/src/interface.jl:59
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ speed(🔹)::Float64
└ @ BinaryTraits ~/.julia/dev/BinaryTraits.jl/src/interface.jl:59
❌ Crane is missing these implementations:
1. BinaryTrait{Fly}: Positive{Fly} ⇢ liftoff(🔹)::Any
2. BinaryTrait{Fly}: Positive{Fly} ⇢ fly(🔹, ::Float64, ::Float64)::Any
3. BinaryTrait{Fly}: Positive{Fly} ⇢ speed(🔹)::Float64
=#

# What about composite traits?
@trait FlySwim with Can{Fly},Can{Swim}

# Define a new type and assign with composite trait
struct Swan end
@assign Swan with Can{Fly}, Can{Swim}

# Check compliance. It should automatically drill down to figure out
# the required interface contracts.
@check(Swan)
#=
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ liftoff(🔹)::Any
└ @ BinaryTraits ~/Library/Mobile Documents/com~apple~CloudDocs/Programming/Julia/BinaryTraits.jl/src/interface.jl:59
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ fly(🔹, ::Float64, ::Float64)::Any
└ @ BinaryTraits ~/Library/Mobile Documents/com~apple~CloudDocs/Programming/Julia/BinaryTraits.jl/src/interface.jl:59
┌ Warning: Missing implementation
│   contract = BinaryTrait{Fly}: Positive{Fly} ⇢ speed(🔹)::Float64
└ @ BinaryTraits ~/Library/Mobile Documents/com~apple~CloudDocs/Programming/Julia/BinaryTraits.jl/src/interface.jl:59
❌ Swan is missing these implementations:
1. BinaryTrait{Fly}: Positive{Fly} ⇢ liftoff(🔹)::Any
2. BinaryTrait{Fly}: Positive{Fly} ⇢ fly(🔹, ::Float64, ::Float64)::Any
3. BinaryTrait{Fly}: Positive{Fly} ⇢ speed(🔹)::Float64
=#


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BallSampler(radius; [options])

A method for sampling isolated points from spatial objects using
a ball neighborhood of given `radius`.

## Options

* `metric`  - Metric for the ball (default to `Euclidean()`)
* `maxsize` - Maximum size of the resulting sample (default to none)
"""
struct BallSampler{T,M} <: AbstractSampler
  radius::T
  metric::M
  maxsize::Union{Int,Nothing}
end

BallSampler(radius; metric=Euclidean(), maxsize=nothing) =
  BallSampler(radius, metric, maxsize)

function sample(object, sampler::BallSampler)
  N = ncoords(object)
  T = coordtype(object)
  radius = sampler.radius
  metric = sampler.metric
  msize  = sampler.maxsize ≠ nothing ? sampler.maxsize : Inf

  # neighborhood search method with ball
  ball = BallNeighborhood(radius, metric)
  searcher = NeighborhoodSearcher(object, ball)

  # pre-allocate memory for coordinates
  coords = MVector{N,T}(undef)

  locations = Vector{Int}()
  notviewed = trues(nelms(object))
  while length(locations) < msize && any(notviewed)
    location = rand(findall(notviewed))
    coordinates!(coords, object, location)

    # neighbors (including the location)
    neighbors = search(coords, searcher)

    push!(locations, location)
    notviewed[neighbors] .= false
  end

  view(object, locations)
end


text:
const DEFINED_DIFFRULES = Tuple{Symbol,Int}[]

struct DiffRule{F} end

(::Type{DiffRule{F}})(args...) where {F} = error("no derivative rule defined for $F with arguments $args")

"""
    @define_diffrule f(x) = :(df_dx(\$x))
    @define_diffrule f(x, y) = :(df_dx(\$x, \$y)), :(df_dy(\$x, \$y))
    ⋮

Define a new differentiation rule for the function `f` and the given arguments, which should
be treated as bindings to Julia expressions.

The RHS should be a function call with a non-splatted argument list, and the LHS should be
the derivative expression, or in the `n`-ary case, an `n`-tuple of expressions where the
`i`th expression is the derivative of `f` w.r.t the `i`th argument. Arguments should be
interpolated wherever they are used on the RHS.

Note that differentiation rules are purely symbolic, so no type annotations should be used.

Examples:

    @define_diffrule cos(x)          = :(-sin(\$x))
    @define_diffrule /(x, y)         = :(inv(\$y)), :(-\$x / (\$y^2))
    @define_diffrule polygamma(m, x) = :NaN,       :(polygamma(\$m + 1, \$x))

"""
macro define_diffrule(def)
    @assert isa(def, Expr) && def.head == :(=)
    rhs = def.args[1]
    lhs = def.args[2]
    @assert isa(rhs, Expr) && rhs.head == :call
    f = rhs.args[1]
    args = rhs.args[2:end]
    rhs.args[1] = :(::Type{DiffRule{$(Expr(:quote, f))}})
    key = (f, length(args))
    in(DEFINED_DIFFRULES, key) || push!(DEFINED_DIFFRULES, key)
    return esc(def)
end

"""
    diffrule(f::Symbol, args...)

Return the derivative expression for `f` at the given argument(s), with the argument(s)
interpolated into the returned expression.

In the `n`-ary case, an `n`-tuple of expressions will be returned where the `i`th expression
is the derivative of `f` w.r.t the `i`th argument.

Examples:

    julia> DiffBase.diffrule(:sin, 1)
    :(cos(1))

    julia> DiffBase.diffrule(:sin, :x)
    :(cos(x))

    julia> DiffBase.diffrule(:sin, :(x * y^2))
    :(cos(x * y ^ 2))

    julia> DiffBase.diffrule(:^, :(x + 2), :c)
    (:(c * (x + 2) ^ (c - 1)), :((x + 2) ^ c * log(x + 2)))
"""
diffrule(f::Symbol, args...) = DiffRule{f}(args...)

"""
    hasdiffrule(f::Symbol, arity::Int)

Return `true` if a differentiation rule is defined for `f` and `arity`, or return `false`
otherwise.

Here, `arity` refers to the number of arguments accepted by `f`.

Examples:

    julia> DiffBase.hasdiffrule(:sin, 1)
    true

    julia> DiffBase.hasdiffrule(:sin, 2)
    false

    julia> DiffBase.hasdiffrule(:-, 1)
    true

    julia> DiffBase.hasdiffrule(:-, 2)
    true

    julia> DiffBase.hasdiffrule(:-, 3)
    false
"""
hasdiffrule(f::Symbol, arity::Int) = in((f, arity), DEFINED_DIFFRULES)

################
# General Math #
################

# unary #
#-------#

@define_diffrule +(x)       = :(  1                                  )
@define_diffrule -(x)       = :( -1                                  )
@define_diffrule sqrt(x)    = :(  inv(2 * sqrt($x))                  )
@define_diffrule cbrt(x)    = :(  inv(3 * cbrt($x)^2)                )
@define_diffrule abs2(x)    = :(  $x + $x                            )
@define_diffrule inv(x)     = :( -abs2(inv($x))                      )
@define_diffrule log(x)     = :(  inv($x)                            )
@define_diffrule log10(x)   = :(  inv($x) / log(10)                  )
@define_diffrule log2(x)    = :(  inv($x) / log(2)                   )
@define_diffrule log1p(x)   = :(  inv($x + 1)                        )
@define_diffrule exp(x)     = :(  exp($x)                            )
@define_diffrule exp2(x)    = :(  exp2($x) * log(2)                  )
@define_diffrule expm1(x)   = :(  exp($x)                            )
@define_diffrule sin(x)     = :(  cos($x)                            )
@define_diffrule cos(x)     = :( -sin($x)                            )
@define_diffrule tan(x)     = :(  1 + tan($x)^2                      )
@define_diffrule sec(x)     = :(  sec($x) * tan($x)                  )
@define_diffrule csc(x)     = :( -csc($x) * cot($x)                  )
@define_diffrule cot(x)     = :( -(1 + cot($x)^2)                    )
@define_diffrule sind(x)    = :(  (π / 180) * cosd($x)               )
@define_diffrule cosd(x)    = :( -(π / 180) * sind($x)               )
@define_diffrule tand(x)    = :(  (π / 180) * (1 + tand($x)^2)       )
@define_diffrule secd(x)    = :(  (π / 180) * secd($x) * tand($x)    )
@define_diffrule cscd(x)    = :( -(π / 180) * cscd($x) * cotd($x)    )
@define_diffrule cotd(x)    = :( -(π / 180) * (1 + cotd($x)^2)       )
@define_diffrule asin(x)    = :(  inv(sqrt(1 - $x^2))                )
@define_diffrule acos(x)    = :( -inv(sqrt(1 - $x^2))                )
@define_diffrule atan(x)    = :(  inv(1 + $x^2)                      )
@define_diffrule asec(x)    = :(  inv(abs($x) * sqrt($x^2 - 1))      )
@define_diffrule acsc(x)    = :( -inv(abs($x) * sqrt($x^2 - 1))      )
@define_diffrule acot(x)    = :( -inv(1 + $x^2)                      )
@define_diffrule asind(x)   = :(  180 / π / sqrt(1 - $x^2)           )
@define_diffrule acosd(x)   = :( -180 / π / sqrt(1 - $x^2)           )
@define_diffrule atand(x)   = :(  180 / π / (1 + $x^2)               )
@define_diffrule asecd(x)   = :(  180 / π / abs($x) / sqrt($x^2 - 1) )
@define_diffrule acscd(x)   = :( -180 / π / abs($x) / sqrt($x^2 - 1) )
@define_diffrule acotd(x)   = :( -180 / π / (1 + $x^2)               )
@define_diffrule sinh(x)    = :(  cosh($x)                           )
@define_diffrule cosh(x)    = :(  sinh($x)                           )
@define_diffrule tanh(x)    = :(  sech($x)^2                         )
@define_diffrule sech(x)    = :( -tanh($x) * sech($x)                )
@define_diffrule csch(x)    = :( -coth($x) * csch($x)                )
@define_diffrule coth(x)    = :( -(csch($x)^2)                       )
@define_diffrule asinh(x)   = :(  inv(sqrt($x^2 + 1))                )
@define_diffrule acosh(x)   = :(  inv(sqrt($x^2 - 1))                )
@define_diffrule atanh(x)   = :(  inv(1 - $x^2)                      )
@define_diffrule asech(x)   = :( -inv($x * sqrt(1 - $x^2))           )
@define_diffrule acsch(x)   = :( -inv(abs($x) * sqrt(1 + $x^2))      )
@define_diffrule acoth(x)   = :(  inv(1 - $x^2)                      )
@define_diffrule deg2rad(x) = :(  π / 180                            )
@define_diffrule rad2deg(x) = :(  180 / π                            )
@define_diffrule gamma(x)   = :(  digamma($x) * gamma($x)            )
@define_diffrule lgamma(x)  = :(  digamma($x)                        )

# binary #
#--------#

@define_diffrule +(x, y) = :(1),                  :(1)
@define_diffrule -(x, y) = :(1),                  :(-1)
@define_diffrule *(x, y) = :($y),                 :($x)
@define_diffrule /(x, y) = :(inv($y)),            :(-$x / ($y^2))
@define_diffrule ^(x, y) = :($y * ($x^($y - 1))), :(($x^$y) * log($x))

@define_diffrule atan2(x, y) = :($y / ($x^2 + $y^2)), :(-$x / ($x^2 + $y^2))
@define_diffrule hypot(x, y) = :($x / hypot($x, $y)), :($y / hypot($x, $y))
@define_diffrule mod(x, y)   = :(first(promote(ifelse(isinteger($x / $y), NaN, 1), NaN))), :(z = $x / $y; first(promote(ifelse(isinteger(z), NaN, -floor(z)), NaN)))
@define_diffrule rem(x, y)   = :(first(promote(ifelse(isinteger($x / $y), NaN, 1), NaN))), :(z = $x / $y; first(promote(ifelse(isinteger(z), NaN, -trunc(z)), NaN)))

####################
# SpecialFunctions #
####################

# unary #
#-------#

@define_diffrule erf(x)         = :(  (2 / sqrt(π)) * exp(-$x * $x)       )
@define_diffrule erfinv(x)      = :(  (sqrt(π) / 2) * exp(erfinv($x)^2)   )
@define_diffrule erfc(x)        = :( -(2 / sqrt(π)) * exp(-$x * $x)       )
@define_diffrule erfcinv(x)     = :( -(sqrt(π) / 2) * exp(erfinv($x)^2)   )
@define_diffrule erfi(x)        = :(  (2 / sqrt(π)) * exp($x * $x)        )
@define_diffrule erfcx(x)       = :(  (2 * x * erfcx($x)) - (2 / sqrt(π)) )
@define_diffrule dawson(x)      = :(  1 - (2 * x * dawson($x))            )
@define_diffrule digamma(x)     = :(  trigamma($x)                        )
@define_diffrule invdigamma(x)  = :(  inv(trigamma(invdigamma($x)))       )
@define_diffrule trigamma(x)    = :(  polygamma(2, $x)                    )
@define_diffrule airyai(x)      = :(  airyaiprime($x)                     )
@define_diffrule airyaiprime(x) = :(  $x * airyai($x)                     )
@define_diffrule airybi(x)      = :(  airybiprime($x)                     )
@define_diffrule airybiprime(x) = :(  $x * airybi($x)                     )
@define_diffrule besselj0(x)    = :( -besselj1($x)                        )
@define_diffrule besselj1(x)    = :(  (besselj0($x) - besselj(2, $x)) / 2 )
@define_diffrule bessely0(x)    = :( -bessely1($x)                        )
@define_diffrule bessely1(x)    = :(  (bessely0($x) - bessely(2, $x)) / 2 )

# TODO:
#
# eta
# zeta
# airyaix
# airyaiprimex
# airybix
# airybiprimex

# binary #
#--------#

@define_diffrule besselj(ν, x)   = :NaN, :(  (besselj($ν - 1, $x) - besselj($ν + 1, $x)) / 2   )
@define_diffrule besseli(ν, x)   = :NaN, :(  (besseli($ν - 1, $x) + besseli($ν + 1, $x)) / 2   )
@define_diffrule bessely(ν, x)   = :NaN, :(  (bessely($ν - 1, $x) - bessely($ν + 1, $x)) / 2   )
@define_diffrule besselk(ν, x)   = :NaN, :( -(besselk($ν - 1, $x) + besselk($ν + 1, $x)) / 2   )
@define_diffrule hankelh1(ν, x)  = :NaN, :(  (hankelh1($ν - 1, $x) - hankelh1($ν + 1, $x)) / 2 )
@define_diffrule hankelh2(ν, x)  = :NaN, :(  (hankelh2($ν - 1, $x) - hankelh2($ν + 1, $x)) / 2 )
@define_diffrule polygamma(m, x) = :NaN, :(  polygamma($m + 1, $x)                             )
@define_diffrule beta(a, b)      = :(beta($a, $b)*(digamma($a) - digamma($a + $b))),
                                            :(beta($a, $b)*(digamma($b) - digamma($a + $b)))
@define_diffrule lbeta(a, b)     = :(digamma($a) - digamma($a + $b)), :(digamma($b) - digamma($a + $b))

# TODO:
#
# zeta
# besseljx
# besselyx
# besselix
# besselkx
# besselh
# besselhx
# hankelh1x
# hankelh2
# hankelh2x

# ternary #
#---------#

# TODO:
#
# besselh
# besselhx


text:
#---------------------------------------
# Imports
#---------------------------------------

import JSON
import Requests

#---------------------------------------
# Constants
#---------------------------------------

COINMARKETCAP_API_URL = "https://api.coinmarketcap.com"

# Indentation (in number of spaces) to use when writing the json wallet file.
JSON_INDENT = 2

# Name of the wallet file.
WALLET_FILENAME = "wallet.json"

#---------------------------------------
# Types
#---------------------------------------

mutable struct ApiCoin
  name::AbstractString
  ticker::AbstractString
  priceusd::Float64
  pricesek::Float64
  percentchange24h::Float64
end

struct WalletCoin
  name::AbstractString
  amount::Float64
end

struct Wallet
  coins::Array{WalletCoin}
end

mutable struct Model
  coindata::Dict{AbstractString, ApiCoin}
  wallet::Wallet
end

struct BuyCmd
  coinname::AbstractString
  amount::Float64
end

struct SellCmd
  coinname::AbstractString
  amount::Float64
end

struct HelpCmd
end

struct NoOpCmd
end

struct DisplayWalletCmd
end

struct SpeculateCmd
  coinname::AbstractString
  priceusd::Float64
end

struct FilterWalletCmd
  coinnames::Array{String}
end

struct TableCell
  text::AbstractString
  align::Symbol
  minsize::Int
  padding::Int
end

struct TableRow
  cells::Array{TableCell}
end

struct Table
  titleleft::AbstractString
  titlecenter::AbstractString
  titleright::AbstractString
  header::TableRow
  rows::Array{TableRow}
  footer::TableRow
end

#---------------------------------------
# Functions
#---------------------------------------

function coinmarketcap(endpoint; params=nothing, version="v1")
  url = "$COINMARKETCAP_API_URL/$version/$endpoint"

  if params != nothing
    url *= "?" * join(["$(p.first)=$(p.second)" for p in params], "&")
  end

  resp = Requests.get(url)
  JSON.parse(IOBuffer(resp.data))
end

function getcoindata(name)
  const FIELD_MAP = Dict([
    :name             => "name"
    :percentchange24h => "percent_change_24h"
    :priceusd         => "price_usd"
    :pricesek         => "price_sek"
    :ticker           => "id"
  ])

  data = coinmarketcap("ticker/$name", params=Dict("convert" => "SEK"))[1]

  values = []
  for (s, t) in zip(fieldnames(ApiCoin), ApiCoin.types)
    value = data[FIELD_MAP[s]]
    if t != AbstractString
      value = parse(t, data[FIELD_MAP[s]])
    end

    push!(values, value)
  end

  ApiCoin(values...)
end

function coinexists(name)
  data = coinmarketcap("ticker/$name")
  !(isa(data, Dict) && haskey(data, "error"))
end

function loadwallet(filename)
  # See https://github.com/JuliaIO/JSON.jl/issues/155 regarding use_mmap=false.
  walletdata = JSON.parsefile(filename, use_mmap=false)

  fields  = fieldnames(WalletCoin) .|> string
  tocoin = t -> WalletCoin([t[s] for s in fields]...)
  coins  = walletdata["coins"] .|> tocoin

  Wallet(coins)
end

function savewallet(wallet, filename)
  open(filename, "w") do f
    JSON.print(f, wallet, JSON_INDENT)
  end
end

function printtable(table)
  cellsize = cell -> max(length(cell.text), cell.minsize) + 2cell.padding
  colwidths = table.header.cells .|> cellsize

  for row in table.rows
    @assert length(row.cells) == length(colwidths)

    for i in 1:length(row.cells)
      cell = row.cells[i]
      colwidths[i] = max(colwidths[i], cell |> cellsize)
    end
  end

  @assert length(table.footer.cells) == length(colwidths)
  for i in 1:length(table.footer.cells)
    cell = table.footer.cells[i]
    colwidths[i] = max(colwidths[i], cell |> cellsize)
  end

  printedge = n -> (print("+"); print("-" ^ n); println("+"))
  printdelim = (n, c="-") -> (print("| "); print(c ^ (n - 2)); println(" |"))

  function printtitle(n, l, c, r)
    n -= length(l) + length(c) + length(r)
    print("|")
    print(l)
    print(" " ^ floor(Int, 0.5n))
    print(c)
    print(" " ^ ceil(Int, 0.5n))
    print(r)
    println("|")
   end

  function printaligned(n, s, a)
    print("|")
    n = max(n - length(s), 0)
    if a == :left
      print(s)
      print(" " ^ n)
    elseif a == :right
      print(" " ^ n)
      print(s)
    else
      print(" " ^ floor(Int, 0.5n))
      print(s)
      print(" " ^ ceil(Int, 0.5n))
    end
  end

  function printrow(row)
    for i in 1:length(row.cells)
      cell = row.cells[i]
      text = (" " ^ cell.padding) * cell.text * (" " ^ cell.padding)
      printaligned(colwidths[i], text, cell.align)
    end

    println("|")
  end

  tablewidth = sum(colwidths) + length(table.header.cells) - 1

  indent = () -> print("  ")

  indent(); printedge(tablewidth)
  indent(); printtitle(tablewidth, table.titleleft, table.titlecenter, table.titleright)
  indent(); printedge(tablewidth)
  indent(); printrow(table.header)
  indent(); printdelim(tablewidth)
  for row in table.rows
    indent(); printrow(row)
  end
  indent(); printdelim(tablewidth, "=")
  indent(); printrow(table.footer)
  indent(); printedge(tablewidth)
end

function parseargs(args)
  s = length(args) >= 1 && args[1] |> lowercase

  if s == false
    DisplayWalletCmd()
  elseif s == "-h" || s == "--help"
    HelpCmd()
  elseif s == "buy" && length(args) == 3 # <amount> <coin>
    try
      amount = args[2] |> a -> parse(Float64, a)
      coin   = args[3] |> lowercase

      BuyCmd(coin, amount)
    catch
      HelpCmd()
    end
  elseif s == "sell" && length(args) == 3 # <amount> <coin>
    try
      amount = args[2] |> a -> parse(Float64, a)
      coin   = args[3] |> lowercase

      SellCmd(coin, amount)
    catch
      HelpCmd()
    end
  elseif s == "speculate"
    try
      coin     = args[2] |> lowercase
      priceusd = args[3] |> a -> parse(Float64, a)

      SpeculateCmd(coin, priceusd)
    catch
      HelpCmd()
    end
  elseif s == "only"
    try
      coins = args[2:end] .|> lowercase

      FilterWalletCmd(coins)
    catch
      HelpCmd()
    end
  elseif s == "wallet"
    DisplayWalletCmd()
  else
    HelpCmd()
  end
end

function updatecoindata!(model)
  for coin in model.wallet.coins
    if !haskey(model.coindata, coin.name)
      model.coindata[coin.name] = ApiCoin(coin.name, "", -1.0, 0.0, 0.0)
    end
  end

  function updatecoin!(coindata)
    if coindata.priceusd < 0.0
      tmp = getcoindata(coindata.name)

      coindata.name             = tmp.name
      coindata.ticker           = tmp.ticker
      coindata.priceusd         = tmp.priceusd
      coindata.pricesek         = tmp.pricesek
      coindata.percentchange24h = tmp.percentchange24h
    end
  end

  @sync begin
    for entry in model.coindata
      @async updatecoin!(entry.second)
    end
  end
end

function execute!(buy::BuyCmd, model)
  !coinexists(buy.coinname) && fail("No such coin exists: '$(buy.coinname)'")

  coins = WalletCoin[]

  isnew = true
  for coin in model.wallet.coins
    if coin.name == buy.coinname
      isnew = false
      coin = WalletCoin(coin.name, coin.amount + buy.amount)
    end

    push!(coins, coin)
  end

  if isnew
    coin = WalletCoin(buy.coinname, buy.amount)
    push!(coins, coin)
  end

  model.wallet = Wallet(coins)
  savewallet(model.wallet, WALLET_FILENAME)

  execute!(DisplayWalletCmd(), model)
end

function execute!(sell::SellCmd, model)
  coins = WalletCoin[]

  for coin in model.wallet.coins
    if coin.name == sell.coinname
      coin = WalletCoin(coin.name, coin.amount - sell.amount)
    end

    if coin.amount > 0.0
      push!(coins, coin)
    end
  end

  model.wallet = Wallet(coins)
  savewallet(model.wallet, WALLET_FILENAME)

  execute!(DisplayWalletCmd(), model)
end

function execute!(speculate::SpeculateCmd, model)
  coindata = getcoindata(speculate.coinname)

  op = coindata.priceusd
  np = speculate.priceusd
  ch = coindata.percentchange24h
  coindata.name *= "*"
  coindata.percentchange24h = 100.0(np / ((1.0 - 0.01ch) * op) - 1.0)
  coindata.pricesek = (coindata.pricesek / coindata.priceusd) * speculate.priceusd
  coindata.priceusd = speculate.priceusd

  model.coindata[speculate.coinname] = coindata

  execute!(DisplayWalletCmd(), model)
end

function execute!(filter::FilterWalletCmd, model)
  model.wallet = Wallet([coin for coin in model.wallet.coins
                         if contains((y, x) -> x.name == y, filter.coinnames, coin)])

  execute!(DisplayWalletCmd(), model)
end

function execute!(::HelpCmd, model)
  s = basename(Base.source_path())
  println("Crypto Wallet v1.0")
  println()
  println("Usage:")
  println("  $s")
  println("  $s -h | --help")
  println()
  println("Options:")
  println("  <none>                     Display wallet conents.")
  println("  buy <amount> <coin>        Add investment to wallet.")
  println("  sell <amount> <coin>       Remove investment from wallet.")
  println("  speculate <coin> <price>   Assign speculative price and display wallet.")
  println("  only <coin>[, <coin> ...]  Only display specified coins.")
end

function execute!(::NoOpCmd, model)
  println("nothing to do")
end

function execute!(::DisplayWalletCmd, model)
  titleleft = ""
  titlecenter = "         W A L L E T  D A T A"
  titleright = Dates.format(now(), "d u yyyy HH:MM ")

  header = TableRow([
    TableCell("Altcoin"     , :left , 0, 1)
    TableCell("Price (USD)" , :right, 0, 1)
    #TableCell("Price (SEK)" , :right, 0, 1)
    TableCell("Change (24h)", :right, 0, 1)
    TableCell("Owned"       , :right, 0, 1)
    TableCell("USD"         , :right, 0, 1)
    TableCell("SEK"         , :right, 0, 1)
    TableCell("%"           , :right, 0, 1)
  ])

  insertsep = s -> length(s) < 3 ? s : insertsep(s[1:end-3]) * " " * s[end-2:end]
  fc  = f -> strip(insertsep(string(floor(Int, f))) * "." * @sprintf("%02d", floor(100(f - floor(f)))))
  fp1 = f -> (f < 0.0 ? "-" : "+") * @sprintf("%10.2f%%", abs(f))
  fp2 = f -> @sprintf("%.2f%%", f)

  usdtotal = 0.0
  sektotal = 0.0

  updatecoindata!(model)

  for coin in model.wallet.coins
    coindata = model.coindata[coin.name]

    usdvalue = coin.amount * coindata.priceusd
    sekvalue = coin.amount * coindata.pricesek

    usdtotal += usdvalue
    sektotal += sekvalue
  end

  rows = []
  for coin in sort(model.wallet.coins, by=coin -> coin.name)
    coindata = model.coindata[coin.name]

    usdvalue = coin.amount * coindata.priceusd
    sekvalue = coin.amount * coindata.pricesek

    row = TableRow([
      TableCell(coindata.name |> lowercase    , :left , 0, 1)
      TableCell(fc(coindata.priceusd)         , :right, 0, 1)
      #TableCell(fc(coindata.pricesek)         , :right, 0, 1)
      TableCell(fp1(coindata.percentchange24h), :right, 0, 1)
      TableCell(fc(coin.amount)               , :right, 0, 1)
      TableCell(fc(usdvalue)                  , :right, 0, 1)
      TableCell(fc(sekvalue)                  , :right, 0, 1)
      TableCell(fp2(100usdvalue/usdtotal)     , :right, 0, 1)
    ])

    push!(rows, row)
  end

  footer = TableRow([
    TableCell("Total"     , :left , 0, 1)
    TableCell("-"         , :right, 0, 1)
    #TableCell("-"         , :right, 0, 1)
    TableCell("-"         , :right, 0, 1)
    TableCell("-"         , :right, 0, 1)
    TableCell(fc(usdtotal), :right, 0, 1)
    TableCell(fc(sektotal), :right, 0, 1)
    TableCell("100.00%"   , :right, 0, 1)
  ])

  printtable(Table(titleleft, titlecenter, titleright, header, rows, footer))
end

function fail(reason, exitcode=1)
  println("[ERROR] " * reason)
  exit(exitcode)
end

function main()
  cd(dirname(Base.source_path()))

  cmd = parseargs(ARGS)

  wallet = Wallet([])
  if isfile(WALLET_FILENAME)
    wallet = loadwallet(WALLET_FILENAME)
  end

  model = Model(Dict(), wallet)

  execute!(cmd, model)
end

#---------------------------------------
# Entry Point
#---------------------------------------

!isinteractive() && main()


text:
mutable struct Momentum <: Optimizer
    xparams::Vector{XVariable}
    v::Vector
    lr::AbstractFloat
    inertia::AbstractFloat
    name::String
    function Momentum(xparams::Vector{XVariable}; lr=1e-4, inertia=0.9)
        num = length(xparams)
        vel = Vector(undef,num)
        for i = 1:num
            c , θ = xparams[i]
            vel[i] = Zeros(typeof(θ.value), θ.shape)
        end
        new(xparams, vel, lr, inertia, "Momentum")
    end
end

function Base.show(io::IO, O::Momentum)
    print("Momentum(lr=$(O.lr), inertia=$(O.inertia))")
end

function update!(O::Momentum; clipfn::Function=LPInfNormClip, clipvalue=10.0)
    vel = O.v
    μ = - O.lr
    ρ = O.inertia

    for i = 1:length(O.xparams)
        c , θ = O.xparams[i]
        ∇ = clipfn(setNanInfZero(θ.delta), clipvalue)
        @. vel[i] = ρ * vel[i] + ∇
        @. θ.value += μ * vel[i]
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
export connectNet

function connectNet(net,connectGraph,layersDims)
  topRangePointer=ones(Int64,size(connectGraph,1))
  for i=1:size(connectGraph,1)
    for j=1:size(connectGraph,2)
      if connectGraph[i,j]==1
        net[i].top=[net[i].top;net[j]]
        net[j].bottom=[net[j].bottom;net[i]]

        newTopRange=[topRangePointer[j]:topRangePointer[j]+size(net[i].rₒ,1)-1;]
        push!(net[i].topRange,newTopRange)
        topRangePointer[j]+=size(net[i].rₒ,1)
      end
    end
  end
end


text:
using BSplineKit

using BandedMatrices
using LinearAlgebra
using Random
using SparseArrays
using Test

import BSplineKit:
    AbstractDifferentialOp,
    DifferentialOpSum,
    LeftNormal,
    RightNormal

import BSplineKit.Recombinations:
    NoUniqueSolutionError,
    num_constraints,
    num_recombined

# Chebyshev (Gauss-Lobatto) points.
gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]

include("diffops.jl")
include("knots.jl")
include("bsplines.jl")
include("splines.jl")
include("natural.jl")
include("interpolation.jl")
include("approximation.jl")
include("recombination.jl")
include("collocation.jl")
include("galerkin.jl")
include("banded_tensors.jl")
include("airy.jl")

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
__precompile__()

module AnimalBreedingTools

using Printf
using Statistics
using LinearAlgebra
using Base.Threads
using StatsBase
using Jacobi

struct AnimalBreedingDataSetRaw
   pedlist::Matrix{Int}
   genid::Vector{Int}
   G::Matrix{Float64}
   y::Vector{Float64}
   X::Matrix{Float64}
end

export get_design_matrix, get_interaction_matrix, directsum, solve_pcg, approximate_trace_of_inverse
export v2r, c2r, r2v, r2c, vech, bending, bending2
export AnimalBreedingDataSetRaw, load_data_set
export chol_symtrid, ldlt_symtrid, takahashi_ldlt_symtrid!
export normalized_legendre, normalized_legendre_matrix, hrv_class
export disp

include("designmatrix.jl")
include("solver.jl")
include("tools.jl")
include("data.jl")
include("polynomials.jl")
include("io.jl")

end


text:
include(joinpath(dirname(@__DIR__), "src", "LinearShiftRegisters.jl"))
using Documenter, .LinearShiftRegisters

Documenter.makedocs(
    clean = true,
    doctest = true,
    modules = Module[LinearShiftRegisters],
    repo = "",
    highlightsig = true,
    sitename = "LinearShiftRegisters Documentation",
    expandfirst = [],
    pages = [
        "Index" => "index.md",
    ]
)

deploydocs(;
    repo  =  "github.com/jakewilliami/LinearShiftRegisters.jl.git",
)


text:
using Test, AnalogMTK

@info("No real test at this time. Simply `using AnalogMTK`")

#Last line


text:
"Provides `YearMonth` type."
module YearMonths

using Dates

export YearMonth

"""
`YearMonth` is an immutable value that represents a tuple (year, month).

Example:
```julia
julia> ym = YearMonth(2018, 12) # represents December of year 2018.
```
"""
struct YearMonth
    y::Int
    m::Int

    function YearMonth(y::Int, m::Int)
        @assert 1 <= m && m <= 12 "Month should be between 1 and 12."
        return new(y, m)
    end
end

YearMonth(dt::Date) = YearMonth(Dates.yearmonth(dt)...)

const RGX_YYYY_MM = r"^[0-9]+-[0-9][0-9]$" # yyyy-mm
const RGX_YYYYMM = r"^[0-9]+[0-9][0-9]$" # yyyymm

"""
    YearMonth(str::AbstractString)

Accepts two formats: "yyyy-mm" or "yyyymm".
"""
function YearMonth(str::AbstractString)
    @assert !isempty(str) "Cannot convert empty string to YearMonth."
    if occursin(RGX_YYYY_MM, str)
        y = parse(Int, str[1:(end-3)])
        m = parse(Int, str[(end-1):end])
        return YearMonth(y, m)
    elseif occursin(RGX_YYYYMM, str)
        y = parse(Int, str[1:(end-2)])
        m = parse(Int, str[(end-1):end])
        return YearMonth(y, m)
    else
        error("Invalid format to create a YearMonth: $str.")
    end
end

YearMonth(ym::Integer) = YearMonth(string(ym))
YearMonth(y::Integer, m::Integer) = YearMonth(Int(y), Int(m))

Dates.year(ym::YearMonth) = ym.y
Dates.month(ym::YearMonth) = ym.m
Dates.yearmonth(ym::YearMonth) = (ym.y, ym.m)

"""
    Date(ym::YearMonth)

Creates a date based on the first day of the month.
"""
Dates.Date(ym::YearMonth) = Date(yearmonth(ym)...)

Base.:+(ym::YearMonth, p::Year) = YearMonth(ym.y + Dates.value(p), ym.m)
Base.:-(ym::YearMonth, p::Year) = YearMonth(ym.y - Dates.value(p), ym.m)

Base.:+(ym::YearMonth, p::Month) = YearMonth(yearmonth(Date(ym) + p)...)
Base.:-(ym::YearMonth, p::Month) = YearMonth(yearmonth(Date(ym) - p)...)

Base.:+(p::P, ym::YearMonth) where {P<:Union{Year, Month}} = +(ym, p)

Base.:<(ym1::YearMonth, ym2::YearMonth) = Date(ym1) < Date(ym2)
Base.:(<=)(ym1::YearMonth, ym2::YearMonth) = ym1 == ym2 || ym1 < ym2
Base.isequal(ym1::YearMonth, ym2::YearMonth) = ym1.y == ym2.y && ym1.m == ym2.m
Base.isless(ym1::YearMonth, ym2::YearMonth) = isless(Date(ym1), Date(ym2))

Dates.firstdayofmonth(ym::YearMonth) = Date(ym)
Dates.lastdayofmonth(ym::YearMonth) = lastdayofmonth(Date(ym))

Base.string(ym::YearMonth) = "$(year(ym))-$(lpad(month(ym), 2, "0"))"
Base.show(io::IO, ym::YearMonth) = print(io, "YearMonth(\"", string(ym), "\")")

end # module YearMonths


text:
using Base.BinaryPlatforms

function augment_platform!(platform)
    # Can't use Preferences since we might be running this very early with a non-existing Manifest
    LAMMPS_UUID = Base.UUID("5b3ab26d-9607-527c-88ea-8fe5ba57cafe")
    abi = get(Base.get_preferences(LAMMPS_UUID), "abi", Sys.iswindows() ? "microsoftmpi" : "mpich")
    if !haskey(platform, "mpi")
        platform["mpi"] = abi
    end
    return platform
end

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
using Distributed
addprocs(8)

using Caesar
# @everywhere using Caesar
using DelimitedFiles, JLD
using IncrementalInference
using KernelDensityEstimatePlotting
using Gadfly, Cairo, Fontconfig

include(joinpath(@__DIR__,"slamUtils.jl"))
include(joinpath(@__DIR__,"expDataUtils.jl"))
include(joinpath(@__DIR__,"plotSASUtils.jl"))

expID = "dock"; datawindow = collect(410:450);
# expID = "drift"; datawindow = collect(1550:1750);

allpaths = getPaths(expID,"range", trialID = 1);
igt = loadGT(datawindow,expID,allpaths);
mfin = loadRanges(datawindow,allpaths)
posData = loadNav(datawindow,allpaths)
dposData = deepcopy(posData)

dposData = deepcopy(posData)
cumulativeDrift!(dposData,[0.0;0],[0.2,0.2])

window = 1:datawindow[end]-datawindow[1];
fg = initfg();
beacon = :l1
addVariable!(fg, beacon, Point2 )

for i in window
    sym = Symbol("x$i")
    addVariable!(fg, sym, Point2)
    initManual!(fg,sym,kde!(rand(MvNormal(dposData[i,:],Diagonal([0.5;0.5].^2)),100)))
end

rtkCov = Matrix(Diagonal([0.1;0.1].^2));
dpσ = Matrix(Diagonal([0.2;0.2].^2))

for i in window
    sym = Symbol("x$i")
    nextsymi = i+1;
    nextsym = Symbol("x$nextsymi")

    if i == window[1]
        pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
        addFactor!(fg, [sym;], pp, autoinit=false)

        dx = dposData[i+1,1] - dposData[i,1];
        dy = dposData[i+1,2] - dposData[i,2];
        p2p2 = Point2Point2(MvNormal([dx;dy],dpσ))
        addFactor!(fg, [sym;nextsym], p2p2, autoinit=false)
    elseif  i == window[end]
        pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
        addFactor!(fg, [sym;], pp, autoinit=false)
    else
        dx = dposData[i+1,1] - dposData[i,1];
        dy = dposData[i+1,2] - dposData[i,2];
        p2p2 = Point2Point2(MvNormal([dx;dy],dpσ))
        addFactor!(fg, [sym;nextsym], p2p2, autoinit=false)
    end
end

rangewindow = window[1]:5:window[end]
for i in rangewindow
    sym = Symbol("x$i")
    ppR = Point2Point2Range(AliasingScalarSampler(mfin[:,1,i],exp.(mfin[:,2,i]),SNRfloor=0.8))
    addFactor!(fg, [sym;beacon], ppR, autoinit=false)
end

writeGraphPdf(fg,viewerapp="", engine="neato", filepath = "/tmp/test.pdf")

getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = false
getSolverParams(fg).async = true
getSolverParams(fg).multiproc = true
getSolverParams(fg).downsolve = true
getSolverParams(fg).maxincidence = 100

tree, smt, hist = solveTree!(fg)
drawTree(tree,filepath = "/tmp/test.pdf")
# fg2 = deepcopy(fg)
# tree, smt, hist = solveTree!(fg,tree)

# plotSASDefault(fg,expID, posData,igt,datadir=allpaths[1],savedir=scriptHeader*"SASdefault.pdf")
plotSASDefault(fg,expID, posData,igt,dposData, datadir=allpaths[1], savedir="/tmp/caesar/test.pdf")
run(`evince /tmp/caesar/test.pdf`)

#PLOTTING ------------

plk= [];
push!(plk,plotBeaconGT(igt));
plotBeaconContours!(plk,fg);

for var in rangewindow
    mysym = Symbol("x$var")
    push!(plk, plotPoint(getVal(fg,mysym), colorIn = colorant"orange"))
end

# all forward convolves
LL = BallTreeDensity[];
IncrementalInference.proposalbeliefs!(fg, :l1, map(x->getFactor(fg,x),ls(fg, :l1)), LL, Dict{Int, Vector{BallTreeDensity}}())

# for i = 1:length(LL)
    # L1ac = rand(LL[i],100)
    # push!(plk,layer(x=L1ac[1,:],y=L1ac[2,:],Geom.histogram2d(xbincount=400, ybincount=400)))
# end

# product of new forward convolves
ll  = IncrementalInference.predictbelief(fg, :l1, ls(fg,:l1))
L1ac = ll[1];
push!(plk,layer(x=L1ac[1,:],y=L1ac[2,:],Geom.histogram2d(xbincount=200, ybincount=200)))
# L1est = manikde(ll, Point2().manifolds)

plotKDEMeans!(plk,fg);
push!(plk,plotPath(dposData,colorIn=colorant"blue"));
push!(plk,plotPath(posData));
if expID == "dock"
    push!(plk, Coord.cartesian(xmin=-40, xmax=140, ymin=-140, ymax=30,fixed=true))
elseif expID == "drift"
    push!(plk, Coord.cartesian(xmin=20, xmax=200, ymin=-220, ymax=0,fixed=true))
end
savefile = "/tmp/caesar/test.pdf"
Gadfly.plot(plk...) |> PDF(savefile)

plk= [];
push!(plk,plotBeaconGT(igt));
plotBeaconContours!(plk,fg);

for var in rangewindow
    mysym = Symbol("x$var")
    push!(plk, plotPoint(getVal(fg,mysym), colorIn = colorant"orange"))
end

L1ac = predictbelief(fg, :l1,ls(fg, :l1))
push!(plk,layer(x=L1ac[1,:],y=L1ac[2,:],Geom.histogram2d(xbincount=400, ybincount=400)))

plotKDEMeans!(plk,fg);
push!(plk,plotPath(posData));
push!(plk,plotPath(dposData,colorIn=colorant"blue"));
if expID == "dock"
    push!(plk, Coord.cartesian(xmin=-40, xmax=140, ymin=-140, ymax=30,fixed=true))
elseif expID == "drift"
    push!(plk, Coord.cartesian(xmin=20, xmax=200, ymin=-220, ymax=0,fixed=true))
end
savefile = scriptHeader*"predBel.pdf"
Gadfly.plot(plk...) |> PDF(savefile)


text:
remaps(mapping, x::Symbol) = get(mapping, x, x)
remaps(mapping, x::Expr) = Expr(x.head, remaps(mapping, x.args[1]), x.args[2:end]...)
remaps(mapping, x::Number) = x

type Op
    name::Symbol
    inputs::Vector
    outputs::Vector
    body::Vector
    info::Expr

    function Op(name, inputs, outputs, body, info, mapping = Dict{Symbol, Symbol}())
        if !isempty(body)
            remap(x) = remaps(mapping, x)
            map!(remap, inputs, inputs)
            uniques = Set{Symbol}()
            for op in body
                map!(remap, op.inputs, op.inputs)
                union!(uniques, filter(o->isa(o, Symbol), op.inputs))
                if isdefined(Symbol("δ$(op.name)_const")) || all(isconst, op.inputs)
                    [mapping[o] = Symbol("$(o)_const") for o in filter(isvar, op.outputs)]
                end
                [mapping[o] = gensym(o) for o in filter(o->o in inputs, op.outputs)]
                map!(remap, op.outputs, op.outputs)
                op.name = name_const(op.name, op.inputs, op.outputs)
            end
            map!(remap, outputs, outputs)
            name = name_const(name, inputs, outputs)
        end
        new(name, inputs, outputs, body, info)
    end
end

function name_const(name, inputs, outputs)
    if isdefined(Symbol("δ$(name)_const")) || all(isconst, inputs) || all(isvar, inputs) || all(isconst, outputs)
        return name
    end
    n = string(name)
    for k = eachindex(inputs)
        if !isvar(inputs[k])
            n *= "_const$k"
        end
    end
    Symbol(n)
end

function parse_function(expr, info, mapping = Dict{Symbol, Symbol}())
    @assert (expr.head == :function || expr.head == :(=)) && length(expr.args) == 2  "Only functions can be differentiated"
    header = expr.args[1]
    @assert header.head == :call "Only functions can be differentiated"
    name = header.args[1]
    inputs = header.args[2:end]
    body = expr.args[2]
    @assert body.head == :block "Body of the function is not found"

    local outputs
    ops = []
    [(info, outputs) = parse_line!(ops, info, line) for line in body.args]
    Op(name, inputs, outputs, ops, info, mapping)
end

parse_line!(ops, info, line::LineNumberNode) = info, [line]
parse_line!(ops, info, line::Symbol) = info, [line]
function parse_line!(ops, info, line::Expr)
    outputs = []
    if line.head == :(=)
        outputs = parse_assign!(ops, info, line.args...)
    elseif line.head == :call || line.head == :(.)
        outputs = [parse_arg!(ops, info, line)]
    elseif line.head == :return && isa(line.args[1], SlotNumber)
        outputs = [Symbol(string(line.args[1]))]
    elseif line.head == :return && isa(line.args[1], Symbol)
        outputs = [line.args[1]]
    elseif line.head == :tuple || line.head == :return && line.args[1].head != :tuple
        outputs = [parse_arg!(ops, info, arg) for arg in line.args]
    elseif line.head == :return && line.args[1].head == :tuple
        outputs = [parse_arg!(ops, info, arg) for arg in line.args[1].args]
    elseif line.head == :line
        info = line
    else error("In function $expr do not know how to handle $(line.head) on $line")
    end
    info, outputs
end

function parse_assign!(ops, info, vals, expr::Symbol)
    @assert vals.head == :tuple "In assignment $vals = $expr do not know how to handle $(vals.head)"
    outputs = [vals.args...]
    push!(ops, Op(:fanout, [expr], outputs, [], info))
    outputs
end

function parse_assign!(ops, info, vals, expr::Expr)
    func, inputs = parse_expr!(ops, info, expr)
    outputs = outs(vals)
    push!(ops, Op(func, inputs, outputs, [], info))
    outputs
end

outs(vals::Symbol) = [vals]
outs(vals::SlotNumber) = [Symbol(string(vals))]
outs(vals::Expr) = [vals.args...]

function parse_expr!(ops, info, expr::Expr)
    if expr.head == :tuple
        args = [parse_arg!(ops, info, arg) for arg in expr.args]
        :tuple, args
    elseif expr.head == :call && expr.args[1] == GlobalRef(Base, :broadcast)
        if isa(expr.args[2], GlobalRef)
            expr.args[2] = expr.args[2].name
        end
        ".$(expr.args[2])", [parse_arg!(ops, info, arg) for arg in expr.args[3:end]]
    elseif expr.head == :call
        if isa(expr.args[1], GlobalRef)
            expr.args[1] = expr.args[1].name
        end
        args = [parse_arg!(ops, info, arg) for arg in expr.args[2:end]]
        while length(args) > 2 && (expr.args[1] == :(+) || expr.args[1] == :(*))
            a = shift!(args)
            arg = Symbol("tmp$(length(ops)+1)")
            push!(ops, Op(opname(expr.args[1]), [a, args[1]], [arg], [], info))
            args[1] = arg
        end
        opname(expr.args[1]), args
    elseif expr.head == :(.)
        @assert expr.args[2].head == :tuple
        ".$(expr.args[1])", [parse_arg!(ops, info, arg) for arg in expr.args[2].args]
    elseif expr.head == :ref
        args = [parse_arg!(ops, info, arg) for arg in expr.args]
        :ref, args
    elseif expr.head == :(:)
        args = [parse_arg!(ops, info, arg) for arg in expr.args]
        :colon, args
    else error("In expr $expr do not know how to handle $(expr.head)")
    end
end

opname(name) = get(opnames, name, name)
const opnames = Dict(:(.*) => :dot_times, :(*) => :times, :(.+) => :dot_plus, :(+) => :plus,
                     :(./) => :dot_divide, :(/) => :divide, :(.-) => :dot_minus, :(-) => :minus,
                     :(.^) => :dot_power, :(^) => :power, :getindex => :ref)

function parse_arg!(ops, info, arg::Expr)
    func, inputs = parse_expr!(ops, info, arg)
    arg = Symbol("tmp$(length(ops)+1)")
    push!(ops, Op(func, inputs, [arg], [], info))
    arg
end

parse_arg!(ops, info, arg::Symbol) = arg
parse_arg!(ops, info, arg::Number) = arg
parse_arg!(ops, info, arg::SlotNumber) = Symbol(string(arg))


text:
@testset "coloring.jl" begin
    @testset "fill_with_generator!(::LegacyGraphColoringGenerator)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        nb_nodes = 10
        probability = 1.4

        generator = SeaPearl.LegacyGraphColoringGenerator(nb_nodes, probability)

        SeaPearl.fill_with_generator!(model, generator; seed=12)

        @test length(keys(model.variables)) == nb_nodes + 1
        @test length(model.constraints) == 24

    end

    @testset "fill_with_generator!(::HomogenousGraphColoringGenerator)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        nb_nodes = 10
        probability = 0.5

        generator = SeaPearl.HomogenousGraphColoringGenerator(nb_nodes, probability)

        SeaPearl.fill_with_generator!(model, generator; seed=12)

        @test length(keys(model.variables)) == nb_nodes + 1
        @test length(model.constraints) == 55

    end

    @testset "fill_with_generator!(::ClusterizedGraphColoringGenerator)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        nb_nodes = 10
        k = 3
        probability = 0.5

        generator = SeaPearl.ClusterizedGraphColoringGenerator(nb_nodes, k, probability)

        SeaPearl.fill_with_generator!(model, generator; seed=12)

        @test length(keys(model.variables)) == nb_nodes + 1

        # This condition is there because of the way random are generated can change from one version to another
        if VERSION >= v"1.6.0"
            @test length(model.constraints) == 38
        end

    end

end


text:
using Extern
using Base.Test
using StatsBase

@extern FakePackage function test1(x)
    return x + 1
end

@test !isdefined(Main, :test1)

@test @extern(FakePackage, 1 + 1) === nothing

# StatsBase is installed and callable from Main (via `using`)
@extern StatsBase type SomeType
    w::WeightVec
    x::Int
end

@test isdefined(Main, :SomeType)

# Package is installed (in test/REQUIRE) but not callable from Main
@test @extern(DataStructures, type OtherType <: OrderedDict end) === nothing
@test !isdefined(Main, :OtherType)


text:
## F updates

"""
    updateF!(M::HePPCATModel,Y,method)

Update the factor matrix in `M` with respect to data `Y` using `method`.
"""
function updateF! end

# Update method: Expectation Maximization
"""
    updateF!(M::HePPCATModel,Y,::ExpectationMaximization)

Expectation maximization update of `F`.
"""
function updateF!(M::HePPCATModel,Y,::ExpectationMaximization)
    n, L = size.(Y,2), length(Y)
    Λ = Diagonal(M.λ)
    Γ = [inv(Λ + M.v[l]*I) for l in 1:L]
    Z = [Γ[l]*sqrt(Λ)*M.U'*Y[l] for l in 1:L]
    num = sum(Y[l]*Z[l]'/M.v[l] for l in 1:L)
    den = sum(Z[l]*Z[l]'/M.v[l] + n[l]*Γ[l] for l in 1:L)

    F = svd((num / den) * M.Vt)
    M.U .= F.U
    M.λ .= F.S.^2
    M.Vt .= F.Vt
    return M
end


text:
@json_read mutable struct TextDocumentIdentifier
    uri::DocumentUri
end

@json_read mutable struct TextDocumentItem
    uri::DocumentUri
    languageId::String
    version::Int
    text::String
end

@json_read mutable struct VersionedTextDocumentIdentifier
    uri::DocumentUri
    version::Int
end

@json_read mutable struct TextDocumentPositionParams
    textDocument::TextDocumentIdentifier
    position::Position
end

@json_read mutable struct DocumentFilter
    language::Union{Nothing,String}
    scheme::Union{Nothing,String}
    pattern::Union{Nothing,String}
end

const DocumentSelector = Vector{DocumentFilter}

mutable struct TextDocumentRegistrationOptions
    documentSelector::DocumentSelector
end

mutable struct TextDocumentChangeRegistrationOptions
    documentSelector::DocumentSelector
    syncKind::Int
end

mutable struct TextDocumentEdit
    textDocument::VersionedTextDocumentIdentifier
    edits::Vector{TextEdit}
end

mutable struct WorkspaceEdit
    changes
    documentChanges::Union{Nothing,Vector{TextDocumentEdit}}
end

@json_read mutable struct DidOpenTextDocumentParams
    textDocument::TextDocumentItem
end

@json_read mutable struct TextDocumentContentChangeEvent
    range::Union{Nothing,Range}
    rangeLength::Union{Nothing,Int}
    text::String
end

@json_read mutable struct DidChangeTextDocumentParams
    textDocument::VersionedTextDocumentIdentifier
    contentChanges::Vector{TextDocumentContentChangeEvent}
end

@json_read mutable struct DidSaveTextDocumentParams
    textDocument::TextDocumentIdentifier
    text::Union{Nothing,String}
end

@json_read mutable struct DidCloseTextDocumentParams
    textDocument::TextDocumentIdentifier
end

mutable struct FileEvent
    uri::String
    _type::Int
end
FileEvent(d::Dict) = FileEvent(d["uri"], d["type"])

@json_read mutable struct DidChangeWatchedFilesParams
    changes::Vector{FileEvent}
end

mutable struct FileSystemWatcher
    globPattern::String
    kind::Union{Nothing,Int}
end

mutable struct DidChangeWatchedFilesRegistrationOptions
    watchers::Vector{FileSystemWatcher}
end

mutable struct WillSaveTextDocumentParams
    textDocument::TextDocumentIdentifier
    reason::Int
end

# TextDocumentIdentifier(d::Dict) = TextDocumentIdentifier(d["uri"])

# TextDocumentItem(d::Dict) = TextDocumentItem(d["uri"], d["languageId"], d["version"], d["text"])

# VersionedTextDocumentIdentifier(d::Dict) = VersionedTextDocumentIdentifier(d["uri"], d["version"])

# TextDocumentPositionParams(d::Dict) = TextDocumentPositionParams(TextDocumentIdentifier(d["textDocument"]), Position(d["position"]))

# function DocumentFilter(d::Dict)
#     language = haskeynotnull(d, "language") ? d["language"] : nothing
#     scheme = haskeynotnull(d, "scheme") ? d["scheme"] : nothing
#     pattern = haskeynotnull(d, "pattern") ? d["pattern"] : nothing
#     return DocumentFilter(language, scheme, pattern)
# end

# DidOpenTextDocumentParams(d::Dict) = DidOpenTextDocumentParams(TextDocumentItem(d["textDocument"]))

# function TextDocumentContentChangeEvent(d::Dict)
#       if length(d) == 1
#           text = d["text"]
#           rangeLength = length(text)
#           lines = split(text, "\n")
#           TextDocumentContentChangeEvent(Range(0, 0, length(lines) - 1, length(lines[end]) - 1), rangeLength, text)
#       else
#           TextDocumentContentChangeEvent(Range(d["range"]), d["rangeLength"], d["text"])
#       end
#   end

# DidChangeTextDocumentParams(d::Dict) = DidChangeTextDocumentParams(VersionedTextDocumentIdentifier(d["textDocument"]), TextDocumentContentChangeEvent.(d["contentChanges"]))

# DidSaveTextDocumentParams(d::Dict) = DidSaveTextDocumentParams(TextDocumentIdentifier(d["textDocument"]), get(d, "text", nothing))

# DidCloseTextDocumentParams(d::Dict) = DidCloseTextDocumentParams(TextDocumentIdentifier(d["textDocument"]))

# function DidChangeWatchedFilesParams(d::Dict)
#     DidChangeWatchedFilesParams(FileEvent.(d["changes"]))
# end

# WillSaveTextDocumentParams(d::Dict) = WillSaveTextDocumentParams(TextDocumentIdentifier(d["textDocument"]), d["reason"])


text:
"""
    DBAclustResult(centers,clustids,result)

Holds results of a DBAclust run.

"""
mutable struct DBAclustResult
    centers::AbstractArray{Sequence}
    clustids::Array{Int}
    converged::Bool
    iterations::Int
    dbaresult::DBAResult
end

"""
    avgseq, results = dba(sequences, [dist=SqEuclidean()]; kwargs...)

Perfoms DTW Barycenter Averaging (DBA) given a collection of `sequences`
and the current estimate of the average sequence.

Example usage:

    x = [1,2,2,3,3,4]
    y = [1,3,4]
    z = [1,2,2,4]
    avg,result = dba([x,y,z])
"""
function dbaclust(
        sequences::AbstractVector{Sequence{N,T}},
        nclust::Int,
        n_init::Int,
        _method::DTWMethod,
        _dist::SemiMetric = SqEuclidean();
        dbalen::Int = 0,
        iterations::Int = 100,
        inner_iterations::Int = 10,
        rtol::Float64 = 1e-4,
        rtol_inner::Float64 = rtol,
        n_jobs::Int =1,
        show_progress::Bool = true,
        store_trace::Bool = true,
        i2min::AbstractVector=[],
        i2max::AbstractVector=[]
    ) where {N,T}

    if n_jobs == 1
      best_result = []
      best_cost = []
      for i=1:n_init
      dbaclust_result = dbaclust_single(sequences,nclust,_method,_dist;dbalen=dbalen,iterations=iterations,inner_iterations=inner_iterations,rtol=rtol,rtol_inner=rtol_inner,show_progress=show_progress,store_trace=store_trace,i2min=i2min,i2max=i2max)
        if isempty(best_cost) || dbaclust_result.dbaresult.cost < best_cost
          best_result = deepcopy(dbaclust_result)
          best_cost = best_result.dbaresult.cost
        end
      end #1:n_init
    else
 #results = Array{DBAclustResult}(n)_
     error("parallelism for dbaclust not implemented yet")

    end # n_jobs
    return best_result
end

"""
    avgseq, results = dbaclust_single(sequences, [dist=SqEuclidean()]; kwargs...)

Perfoms a single DTW Barycenter Averaging (DBA) given a collection of `sequences`
and the current estimate of the average sequence.

Example usage:

    x = [1,2,2,3,3,4]
    y = [1,3,4]
    z = [1,2,2,4]
    avg,result = dba([x,y,z])
"""
function dbaclust_single(
        sequences::AbstractVector{Sequence{N,T}},
        nclust::Int,
        _method::DTWMethod,
        _dist::SemiMetric = SqEuclidean();
        init_centers::AbstractVector{Sequence{N,T}} = dbaclust_initial_centers(sequences, nclust, _method, _dist),
        dbalen::Int = 0,
        iterations::Int = 100,
        inner_iterations::Int = 10,
        rtol::Float64 = 1e-4,
        rtol_inner::Float64 = rtol,
        show_progress::Bool = true,
        store_trace::Bool = true,
        i2min::AbstractVector =[],
        i2max::AbstractVector = []
    ) where {N,T}

    # rename for convienence
    avgs = init_centers
    dbalen = length(avgs[1])

    # check initial centers have the same length
    if !all([ length(a) for a in avgs ] .== dbalen)
        throw(ArgumentError("all initial centers should be the same length"))
    end

    # dimensions
    nseq = length(sequences)
    maxseqlen = maximum([ length(s) for s in sequences ])

    # initialize procedure for computing DTW
    dtwdist = DTWDistance(_method, _dist)

    # TODO switch to ntuples?
    counts = [ zeros(Int,dbalen) for _ in 1:nclust ]
    sums = [ Sequence(Array{T}(dbalen)) for _ in 1:nclust ]

    # cluster assignments for each sequence
    clus_asgn = Array{Int}(nseq)
    c = 0

    # arrays storing path through dtw cost matrix
    i1,i2 = Int[],Int[]

    # variables storing optimization progress
    converged = false
    iter = 0
    inner_iter=0
    converged_inner=false
    last_cost = Inf
    total_cost = 0.0
    cost_trace = Float64[]
    costs = Array{Float64}(nseq)

    # main loop ##
    if show_progress
      prog = ProgressMeter.ProgressThresh(rtol)
      ProgressMeter.update!(prog, Inf; showvalues =[(:iteration,iter),
                                                 (Symbol("max iteration"),iterations),
                                                 (:cost,total_cost)])
    end#showprogress

    while !converged && iter < iterations
        show_progress && println("1")

        # first, update cluster assignments based on nearest
        # centroid (measured by dtw distance). Keep track of
        # total cost (sum of all distances to centers).
        total_cost = 0.0
        for s = 1:nseq
            show_progress && print("*")
            # process sequence s
            seq = sequences[s]

            # find cluster assignment for s
            costs[s] = Inf
            for c_ = 1:nclust
                # if one of the two is empty, use unconstrained window. If both are nonempty, but not the same lenght, distpath will throw error
                if isempty(i2min) && isempty(i2max)
                   cost, i1_, i2_ = distpath(dtwdist, avgs[c_], seq)
                else
                   cost, i1_, i2_ = distpath(dtwdist, avgs[c_], seq,i2min,i2max)
                end
                if cost < costs[s]
                    # store cluster, and match indices
                    c = c_
                    i1 = i1_
                    i2 = i2_
                    costs[s] = cost
                end
            end

            # s was assigned to cluster c
            clus_asgn[s] = c
            cnt = counts[c]
            sm = sums[c]
            avg = avgs[c]

            # update stats for barycentric average for
            # the assigned cluster
            for t=1:length(i2)
                cnt[i1[t]] += 1
                sm[i1[t]] += seq[i2[t]]
            end
        end

        show_progress && println("2")

        # if any centers are unused, and reassign them to the sequences
        # with the highest cost
        unused = setdiff(1:nclust, unique(clus_asgn))
        if ~isempty(unused)
            println("unused")
            # reinitialize centers
            for c in unused
                avgs[c] = deepcopy(sequences[indmax(costs)])
                for s = 1:nseq
                    if isempty(i2min) && isempty(i2max)
                      cost, = distpath(dtwdist, avgs[c], seq)
                    else
                      cost, = distpath(dtwdist, avgs[c], seq,i2min,i2max)
                    end
                    if costs[s] > cost
                        costs[s] = cost
                    end
                end
            end
            # we need to reassign clusters, start iteration over
            continue
        end

        # store history of cost while optimizing (optional)
        total_cost = sum(costs)
        store_trace && push!(cost_trace, total_cost)

        # check convergence
        Δ = (last_cost-total_cost)/total_cost
        if Δ < rtol
            converged = true
        else
            last_cost = total_cost
        end

        # update barycenter estimates
        for (a,s,c) in zip(avgs,sums,counts)
            for t = 1:dbalen
                c[t] == 0 && continue
                a[t] = s[t]/c[t]
            end
            # zero out sums and counts for next iteration
            scale!(s,0)
            scale!(c,0)
        end

        # add additional inner dba iterations

        for i = 1:nclust
            seqs = view(sequences, clus_asgn .== i)
            inner_iter = 0
            converged_inner = false
            oldcost = 1.0e100
            while !converged_inner  && inner_iter < inner_iterations
                show_progress && print("!")
                newcost = dba_iteration!(sums[i], avgs[i], counts[i], seqs, dtwdist; i2min=i2min,i2max=i2max)
                copy!(avgs[i], sums[i])
                inner_iter += 1
                δ = (oldcost-newcost)/oldcost
                if δ < rtol_inner
                  converged_inner = true
                else
                  oldcost=newcost
                end
            end
            show_progress && println(" ")
        end

        # update progress bar
        iter += 1
        show_progress && ProgressMeter.update!(prog, Δ; showvalues =[(:iteration,iter),
                                                 (Symbol("max iteration"),iterations),
                                                 (:cost,total_cost)])
    end

    return DBAclustResult(avgs, clus_asgn, converged, iter, DBAResult(total_cost,converged_inner,inner_iter,cost_trace))
end

# Wrapper for AbstractArray of one-dimensional time series.
function dbaclust( s::AbstractArray, args...; kwargs... )
    dbaclust(_sequentize(s), args...; kwargs...)
end

"""
   dbaclust_initial_centers(sequences, nclust, dist)

Uses kmeans++ (but with dtw distance) to initialize the centers
for dba clustering.
"""
function dbaclust_initial_centers(
        sequences::AbstractVector{Sequence{N,T}},
        nclust::Int,
        _method::DTWMethod,
        _dist::SemiMetric = SqEuclidean();
    ) where {N,T}

    # procedure for calculating dtw
    dtwdist = DTWDistance(_method,_dist)

    # number of sequences in dataset
    nseq = length(sequences)

    # distance of each datapoint to each center
    dists = zeros(nclust, nseq)

    # distances to closest center
    min_dists = zeros(1, nseq)

    # choose a center uniformly at random
    center_ids = zeros(Int,nclust)
    center_ids[1] = rand(1:nseq)

    # assign the rest of the centers
    for c = 1:(nclust-1)

        # first, compute distances for the previous center
        cent = sequences[center_ids[c]]
        for (i,seq) in enumerate(sequences)
            # this distance will be zero
            i == center_ids[c] && continue
            # else, compute dtw distance
            dists[c,i], = distpath(dtwdist, seq, cent)
        end

        # for each sequence, find distance to closest center
        minimum!(min_dists, dists[1:c,:])

        # square distances
        map!((x)->x^2, min_dists, min_dists)

        # sample the next center
        center_ids[c+1] = sample(1:nseq, Weights(view(min_dists,:)))
    end

    # return list of cluster centers
    return [ deepcopy(sequences[c]) for c in center_ids ]
end

# Wrapper for AbstractArray of one-dimensional time series.
function dbaclust_initial_centers(
        s::AbstractArray,
        nclust::Int,
        _method::DTWMethod,
        _dist::SemiMetric = SqEuclidean();
    )
    dbaclust_initial_centers(_sequentize(s), nclust, _method, _dist)
end


text:
using Hyperscript
using JSServe, Observables
using JSServe: Application, Session, evaljs, linkjs, update_dom!, div, active_sessions
using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session
using JSServe.DOM

function dom_handler(session, request)
    s1 = Slider(1:100)
    s2 = Slider(1:100)
    b = Button("hi")
    t = TextField("lol")
    s_value = s1.value
    linkjs(session, s1.value, s2.value)
    onjs(session, s1.value, js"(v)=> console.log(v)")
    on(t) do text
        println(text)
    end
    return DOM.div(s1, s2, b, t)
end;

# with_session() do session
#     dom_handler(session, nothing)
# end |> display
#
# pop!(JSServe.global_application[].routes.table)
# JSServe.route!(JSServe.global_application[], "/test") do ctx
#     JSServe.serve_dom(ctx, dom_handler)
# end

app = JSServe.Application(
    dom_handler,
    get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1"),
    parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081")),
    verbose = false
)
close(app)


text:
using SparseArrays

list_test_str = split(
"6,10
0,14
9,10
0,3
10,4
4,11
6,0
6,12
4,1
0,13
10,12
3,4
3,0
8,4
1,10
2,14
8,10
9,0

fold along y=7
fold along x=5",
'\n'
)

list_str = let input = @__DIR__() * "/day13_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list)
    n = findfirst(==(""), list)
    dots = fill(0, n-1, 2)
    folds = fill(0, length(list)-n, 2)
    for (i, s) in enumerate(@view(list[1:n-1]))
        dots[i, :] .= parse.(Int, split(s, ','))
    end
    for (i, s) in enumerate(@view(list[n+1:end]))
        folds[i, 1] = s[12] == 'x' ? 1 : 2
        folds[i, 2] = parse(Int, s[14:end])
    end
    fd = folds[1, 1]
    fl = folds[1, 2]
    for dt in eachrow(dots)
        if dt[fd] > fl
            dt[fd] = 2 * fl - dt[fd]
        end
    end
    return length(unique(eachrow(dots)))
end

testit("Test Part 1", part1(list_test_str), ==(17))
testit("Challenge Part 1", part1(list_str), ==(712))

# Part 2

function part2(list)
    n = findfirst(==(""), list)
    dots = fill(0, n-1, 2)
    folds = fill(0, length(list)-n, 2)
    for (i, s) in enumerate(@view(list[1:n-1]))
        dots[i, :] .= parse.(Int, split(s, ','))
    end
    for (i, s) in enumerate(@view(list[n+1:end]))
        folds[i, 1] = s[12] == 'x' ? 1 : 2
        folds[i, 2] = parse(Int, s[14:end])
    end
    m = maximum(dots, dims=1)
    for i in 1:size(folds, 1)
        fd = folds[i, 1]
        fl = folds[i, 2]
        for dt in eachrow(dots)
            if dt[fd] > fl
                dt[fd] = 2 * fl - dt[fd]
            end
        end
        m[fd] = fl
    end
    nx, ny = m
    v = spzeros(Int, ny, nx)
    for (j, i) in eachrow(dots)
        v[i+1, j+1] = 1
    end
    return v
end

v_test = part2(list_test_str)
display(v_test) # output O
v = part2(list_str)
display(v) # output BLHFJPJF


text:
"""
A 'Zipfian traveller population': interactors are located on a two-dimensional
lattice. Each interactor has a 'home' cell on the lattice, but also has the
option of travelling to another cell. The probability to travel to a cell
with coordinates ``x'`` and ``y'`` scales as the inverse of the distance
between that cell and the interactor's home (coordinates ``x`` and ``y``),
measured as the Manhattan distance ``d``:

```math
P \\propto \\frac{1}{(1 + d)^\\lambda}
```

where ``\\lambda`` is a scaling parameter. Interactors only rendezvous with
other interactors currently residing in the same cell.
"""
mutable struct ZipfTravellerPopulation{T<:AbstractInteractor} <: AbstractLatticePopulation
  # Census
  census::Array{T}

  # Spatial census
  spatial_census::Dict{Point{Int}, Array{Int}}

  # Interactors' homes
  coordinates::Array{Point{Int}}

  # Interactors' temporary (within-epoch) locations
  locations::Array{Point{Int}}

  # Lattice size
  size::Int

  # Decay exponent
  de::Float64

  # Probabilities
  probabilities::Array{Float64}

  # Distance function
  distance::Function

  # Range of possible distances
  possible_distances::Dict{Point{Int}, Array{Int}}

  # How many actions have occurred
  nactions::Int

  # Cells at a given distance from given cell
  lookup::Dict{Tuple{Point{Int}, Int}, Array{Point{Int}}}
end

"""
    ZipfTravellerPopulation(K::Int, lambda::Float64)

Construct a Zipf traveller population of size `K` by `K` with scaling
parameter `lambda`.
"""
function ZipfTravellerPopulation{T}(K::Int, lambda::Float64) where {T<:AbstractInteractor}
  # Possible distances on lattice
  possdist = collect(0:2*(K-1))

  # Possible distances for each cell
  possdistcell = Dict{Point{Int}, Array{Int}}()

  # Spatial census
  spatcens = Dict{Point{Int}, Array{Int}}()

  # Make lookup table
  lup = Dict{Tuple{Point{Int}, Int}, Array{Point{Int}}}()
  for x in 1:K
    for y in 1:K
      point = Point{Int}(x, y)
      for dist in possdist
        for i in 1:K
          for j in 1:K
            point2 = Point{Int}(i, j)
            if manhattan(point, point2) == dist
              if !haskey(possdistcell, point)
                possdistcell[point] = Array{Int}(undef, 0)
              end
              push!(possdistcell[point], dist)
              if !haskey(lup, (point, dist))
                lup[(point, dist)] = Array{Point{Int}}(undef, 0)
              end
              push!(lup[(point, dist)], point2)
            end
          end
        end
      end
      possdistcell[point] = unique(possdistcell[point])
      spatcens[point] = []
    end
  end

  # Construct
  ZipfTravellerPopulation(Array{T}(undef, 0), spatcens, Array{Point{Int}}(undef, 0), Array{Point{Int}}(undef, 0), K, lambda, zipf(2*(K-1) + 1, lambda), manhattan, possdistcell, 0, lup)
end

"""
    inject!(x::AbstractInteractor, y::ZipfTravellerPopulation)

Insert an interactor into a Zipf traveller population.
"""
function inject!(x::AbstractInteractor, y::ZipfTravellerPopulation)
  push!(y.census, x)
  home = Point{Int}(rand(1:y.size), rand(1:y.size))
  push!(y.coordinates, home)
  push!(y.locations, home)
  push!(y.spatial_census[home], length(y.census))
end

"""
    eject!(x::Int, y::ZipfTravellerPopulation)

Remove an interactor from a Zipf traveller population.
"""
function eject!(x::Int, y::ZipfTravellerPopulation)
  home = y.coordinates[x]
  loca = y.locations[x]
  splice!(y.census, x)
  splice!(y.coordinates, x)
  splice!(y.locations, x)
  y.spatial_census[loca] = setdiff(y.spatial_census[loca], x)
  for k in keys(y.spatial_census)
    for i in 1:length(y.spatial_census[k])
      if y.spatial_census[k][i] > x
        y.spatial_census[k][i] -= 1
      end
    end
  end
end

"""
    travel!(x::Int, y::ZipfTravellerPopulation)

Make interactor with index `x` travel in a Zipf traveller population.
"""
function travel!(x::Int, y::ZipfTravellerPopulation)
  home = y.coordinates[x]
  loca = y.locations[x]
  possdist = y.possible_distances[home]
  dist = Distributions.wsample(possdist, zipf(maximum(possdist) + 1, y.de))
  y.locations[x] = rand(y.lookup[(home, dist)])

  if y.locations[x] != loca
    y.spatial_census[loca] = setdiff(y.spatial_census[loca], x)
    push!(y.spatial_census[y.locations[x]], x)
  end
end

"""
    travel_all!(x::ZipfTravellerPopulation)

Make every interactor in a Zipf traveller population travel.
"""
function travel_all!(x::ZipfTravellerPopulation)
  for s in 1:length(x.census)
    travel!(s, x)
  end
end

"""
    travel!(x::ZipfTravellerPopulation)

Make a random interactor in a Zipf traveller population travel.
"""
function travel!(x::ZipfTravellerPopulation)
  travel!(rand(1:length(x.census)), x)
end

"""
    rendezvous!(x::ZipfTravellerPopulation)

Rendezvous in a Zipf traveller population (between two random interactors
currently in the same cell).
"""
function rendezvous!(x::ZipfTravellerPopulation)
  length(x.census) > 1 || return false
  int1 = rand(1:length(x.census))
  location = x.locations[int1]
  length(x.spatial_census[location]) > 1 || return false
  friends = setdiff(x.spatial_census[location], int1)
  int2 = rand(friends)
  rendezvous!(x, int1, int2)
end