text:
#=
Functions to estimate the blocking (factor) of scrap parts:

Here we assume a main line with length A_sum (which represents the cross
sectional surface of the scrap shaft) and n_p lines with length A_p which
represents the exposed surface of the individual scrap parts.
So A_sum/(A_p * n_p) represents the occupied (area) fraction (occ_frac),
determined by the radiation relevant area of earch scrap layer cell.

We distibute the n_p lines on the main line regulary and then randomly shift
the offsets of the line in the range between half the gap size, between the
original n_p lines. In this way we create a random distances, but avoid
random clustering effects and very big gaps, which  propably will not occur
in a real scrap bulk.

We then compare two randomly shiftet main lines, to each other. Under the
assumption that they are infinitely close to each other we can then calculate
the overlap factor. From this we calculate a "seen amout" representing the amout
of surface, which contributes to radiation heat exchange between a 1d scrap cell
i and its k neighboars (i-k:i-1,i+1:i+k).

As we see here due to symmetry between the lines (as each line has its own
"movement space") the number of lines does not influence the results.

Some results are shown under ./ex/blocking_eval.jl
=#
################################################################################
################################################################################

function oc_blocking_over_layers(oc_fracs, n_layers::Int = 10)
    nps = [1] # only 1 needed - since symmetry ...
    mean_overlap_oc_frac = get_meanoverlap_over_oc_frac(nps, oc_fracs)
    layers = collect(0:1:n_layers)
    layers_between = layers .-1

    blocking_oc_over_layers = zeros(length(oc_fracs), length(layers))

    for (i,oc_frac) in enumerate(oc_fracs)
        for j = 1:length(layers)
            if j == 1
                blocking_oc_over_layers[i,j] = 1
            elseif j == 2
                blocking_oc_over_layers[i,j] = mean_overlap_oc_frac[i,end]
            else
                blocking_oc_over_layers[i,j] = (1-mean_overlap_oc_frac[i,end]).^layers_between[j] * mean_overlap_oc_frac[i,end]
            end
        end
    end

    return blocking_oc_over_layers
end

function get_meanoverlap_over_oc_frac(nps, oc_frac)
    A_sum = 1000.0     # fixed values
    ### investigate occupied_frac
    progress = Progress(length(nps)*length(oc_frac), 1)

    update!(progress,0)
    mean_overlap_oc_frac = zeros(length(oc_frac),length(nps))

    for (j,n_p) in enumerate(nps)
        println("n_p = ", n_p)

        for i = 1:length(oc_frac)
            A_p_sum = A_sum * oc_frac[i]
            A_p = A_p_sum/n_p
            mean_overlap_oc_frac[i,j] = calc_mean_overlap(A_sum, A_p, n_p, 20_000)
        end

    end

    return mean_overlap_oc_frac
end

function plot_overlap_example(A_sum, A_p, n_p)
    A_p_rel = A_p/A_sum

    p1 = get_rand_pos(1, A_p_rel, n_p)
    p2 = get_rand_pos(1, A_p_rel, n_p)

    lines_p1 = line_from_pos.(p1,A_p_rel)
    y_lines_p1 = y_line_from_pos(p1, 1.0)

    lines_p2 = line_from_pos.(p2,A_p_rel)
    y_lines_p2 = y_line_from_pos(p2, 2.0)

    lines_overlap, overlap = calc_overlap_plot(p1,p2,A_p_rel)
    y_lines_overlap = y_overlapline(lines_overlap, 1.5)

    lines_overlap = filter_missing(lines_overlap)
    y_lines_overlap = filter_missing(y_lines_overlap)

    x1_line_overlap, x2_line_overlap, y_line_overlap =
            vertical_overlap_lines(lines_overlap, 1.0, 2.0)

    pl = plot(lines_p1, y_lines_p1, label ="", linecolor = :blue, linewidth = 3.0)
    plot!(pl, lines_p2, y_lines_p2,label ="", linecolor = :red, linewidth = 3.0)
    plot!(pl, lines_overlap, y_lines_overlap,label ="", linecolor = :green,
            linewidth = 3.0)
    plot!(pl, x1_line_overlap, y_line_overlap,label ="", linecolor = :lightgray,
            linestyle = :dash,linewidth = 1.0)
    plot!(pl, x2_line_overlap, y_line_overlap,label ="", linecolor = :lightgray,
            linestyle = :dash,linewidth = 1.0)
    plot!(pl, xlim=(0,1), ylim = (0.9,2.1))
    plot!(pl, yticks=([1, 1.5, 2], [L"A_{rad,\, layer \, 1}", L"overlap",
            L"A_{rad,\, layer \, 2}"]), bottom_margin = 5px, xlabel = "relative shaft cross section surface area")
    # display(pl)

    mean_mean_overlap = mean(overlap)

    println("Mean overlap: ", mean_mean_overlap)
    return pl
end

function filter_missing(lines::Vector{Vector{Union{Missing, Float64}}})

    lines_new = Vector{Vector{Float64}}(undef,0)

    for i = 1:length(lines)
        if !ismissing(lines[i][1]) && !ismissing(lines[i][2])
            push!(lines_new, [lines[i][1], lines[i][2]])
        end
    end

    return lines_new
end

function vertical_overlap_lines(lines_overlap, yp1, yp2)

    n = length(lines_overlap)
    x1_line_overlap = fill(Union{Float64,Missing}[0, 0],n)
    y_line_overlap = fill(Union{Float64,Missing}[0, 0],n)
    x2_line_overlap = fill(Union{Float64,Missing}[0, 0],n)

    for i = 1:length(lines_overlap)
        x1_line_overlap[i] = [lines_overlap[i][1], lines_overlap[i][1]]
        x2_line_overlap[i] = [lines_overlap[i][2], lines_overlap[i][2]]
        y_line_overlap[i] = [yp1, yp2]
    end

    return x1_line_overlap, x2_line_overlap, y_line_overlap
end

function line_from_pos(p_pos, A_p)
    return [p_pos - 0.5*A_p, p_pos + 0.5*A_p]
end

function y_line_from_pos(p_pos, yval)
    y_line_vec = ones(length(p_pos)) .* yval

    return y_val_vec.(y_line_vec)
end

function y_overlapline(overlap_lines, yval)
    n = length(overlap_lines)
    y_line_vec  = fill(Vector{Union{Missing,Float64}}(missing, 2),n)

    for i = 1:n
        if ismissing(overlap_lines[i][1])
            y_line_vec[i] = [missing, missing]
        else
            y_line_vec[i] = [yval, yval]
        end
    end

    return y_line_vec
end

function y_val_vec(yval)
    return [yval, yval]
end

function calc_overlap_plot(p1::Vector{T}, p2::Vector{T}, A_p::T) where T <: Real
    n = length(p1)
    lines_overlap = fill(Vector{Union{Missing,Float64}}(missing, 2),n)
    overlap = zeros(n)

    for i=1:n
        lines_overlap[i], overlap[i] = calc_overlap_plot(p1[i],p2[i],A_p)
    end

    return lines_overlap, overlap
end

function calc_overlap_plot(p1::Real, p2::Real, A_p::Real)
    overlap_tol = 1E-8

    plow, overlap = range_p_overlap(range_p(p1,A_p), range_p(p2,A_p), A_p)

    if !ismissing(plow)
        overlap_line = [plow, plow+overlap]
    else
        overlap_line = [missing, missing]
    end

    overlap = overlap / A_p  # in percent

    return overlap_line, overlap
end

function calc_mean_overlap(A_sum::Float64, A_p::Float64, n_p::Int64, n_runs = 20)
    mean_mean_overlap = zeros(n_runs)

    @inbounds for i = 1:n_runs
        p1 = get_rand_pos(A_sum, A_p, n_p)
        p2 = get_rand_pos(A_sum, A_p, n_p)
        overlap = calc_overlap(p1,p2,n_p,A_p)
        mean_mean_overlap[i] = mean(overlap)
    end

    mean_overlap = mean(mean_mean_overlap)
    return mean_overlap
end

function calc_overlap(p1, p2, n_p, A_p)
    overlap_tol = 1E-8
    overlap = zeros(Float64, n_p)

    @inbounds for i = 1:n_p
        local plow
        plow, overlap[i] = range_p_overlap(range_p(p1[i],A_p), range_p(p2[i],A_p), A_p)
        # in percent
        overlap[i] = overlap[i] / A_p
    end

    if ( (sum(overlap .< -overlap_tol) >= 1) ||  (sum(overlap .> 1+overlap_tol) >= 1) )
        println( overlap[(overlap .< 0) .| (overlap .> 1)] )
        error("Error in overlap calulation")
    end

    return overlap
end

function range_p(p_pos, A_p)
    return (p_pos - 0.5*A_p, p_pos + 0.5*A_p)
end

function range_p_overlap(p1_range, p2_range, A_p)
    if p1_range[1] > p2_range[1]
        # p1_min > p2_min
        if p1_range[1] > p2_range[2]
            # p1_min > p2_max
            return missing, 0
        end

        return p1_range[1], (A_p - (p1_range[1] - p2_range[1]))

    elseif p1_range[1] < p2_range[1]
        # p2_min > p1_min
        if p1_range[2] < p2_range[1]
            # p1_min > p2_max
            return missing, 0
        end

        return p2_range[1],(A_p - (p2_range[1] - p1_range[1]))
    else
        # equal
        return p1_range[1], A_p
    end
end

function get_rand_pos(A_sum, A_p, n_p)
    off_min = 0.5 * A_p
    off_gap =  (A_sum-A_p*n_p)/n_p
    off_regular = off_gap+2*off_min
    p_pos = collect(off_min+off_gap/2:off_regular:A_sum-off_gap/2)
    dist = zeros(Float64, n_p)

    @inbounds for i=1:n_p
        p_pos[i] = p_pos[i] + rand(-off_gap/2:(off_gap/50):off_gap/2)

        if i>1
            dist[i] = p_pos[i] - p_pos[i-1]
        end
    end

    error_check = falses(n_p)
    error_check = dist .< off_min-1E-16

    if sum(error_check[2:end]) != 0
        error("Something went wrong generating random pos")
    end

    return p_pos
end


text:
"""
Implementation of Genetic Algorithm

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `crossoverRate`: The fraction of the population at the next generation, not including elite children, that is created by the crossover function.
- `mutationRate`: Probability of chromosome to be mutated
- `ɛ`/`epsilon`: Positive integer specifies how many individuals in the current generation are guaranteed to survive to the next generation. Floating number specifies fraction of population.
- `selection`: [Selection](@ref) function
- `crossover`: [Crossover](@ref) function (default: `identity`)
- `mutation`: [Mutation](@ref) function (default: `identity`)
"""
mutable struct GA <: AbstractOptimizer
    populationSize::Int
    crossoverRate::Float64
    mutationRate::Float64
    ɛ::Real
    selection::Function
    crossover::Function
    mutation::Function

    GA(; populationSize::Int=50, crossoverRate::Float64=0.8, mutationRate::Float64=0.1,
        ɛ::Real=0, epsilon::Real=ɛ,
        selection::Function = ((x,n)->1:n),
        crossover::Function = identity, mutation::Function = identity) =
        new(populationSize, crossoverRate, mutationRate, epsilon, selection, crossover, mutation)
end
population_size(method::GA) = method.populationSize
default_options(method::GA) = (iterations=1000, abstol=1e-15)
summary(m::GA) = "GA[P=$(m.populationSize),x=$(m.crossoverRate),μ=$(m.mutationRate),ɛ=$(m.ɛ)]"
show(io::IO,m::GA) = print(io, summary(m))

mutable struct GAState{T,IT} <: AbstractOptimizerState
    N::Int
    eliteSize::Int
    fitness::T
    fitpop::Vector{T}
    fittest::IT
end
value(s::GAState) = s.fitness
minimizer(s::GAState) = s.fittest

"""Initialization of GA algorithm state"""
function initial_state(method::GA, options, objfun, population)
    T = typeof(value(objfun))
    N = length(first(population))
    fitness = zeros(T, method.populationSize)

    # setup state values
    eliteSize = isa(method.ɛ, Int) ? method.ɛ : round(Int, method.ɛ * method.populationSize)

    # Evaluate population fitness
    fitness = map(i -> value(objfun, i), population)
    minfit, fitidx = findmin(fitness)

    # setup initial state
    return GAState(N, eliteSize, minfit, fitness, copy(population[fitidx]))
end

function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::GA, options, itr) where {IT}
    @unpack populationSize,crossoverRate,mutationRate,ɛ,selection,crossover,mutation = method
    evaltype = options.parallelization

    offspring = similar(population)

    # Select offspring
    selected = selection(state.fitpop, populationSize)

    # Perform mating
    offidx = randperm(populationSize)
    offspringSize = populationSize - state.eliteSize
    for i in 1:2:offspringSize
        j = (i == offspringSize) ? i-1 : i+1
        if rand() < crossoverRate
            offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]])
        else
            offspring[i], offspring[j] = population[selected[i]], population[selected[j]]
        end
    end

    # Elitism (copy population individuals before they pass to the offspring & get mutated)
    fitidxs = sortperm(state.fitpop)
    for i in 1:state.eliteSize
        subs = offspringSize+i
        offspring[subs] = copy(population[fitidxs[i]])
    end

    # Perform mutation
    for i in 1:offspringSize
        if rand() < mutationRate
            mutation(offspring[i])
        end
    end

    # Create new generation & evaluate it
    for i in 1:populationSize
        population[i] = apply!(constraints, offspring[i])
    end
    # calculate fitness of the population
    value!(objfun, state.fitpop, population)
    # apply penalty to fitness
    penalty!(state.fitpop, constraints, population)

    # find the best individual
    minfit, fitidx = findmin(state.fitpop)
    state.fittest = population[fitidx]
    state.fitness = state.fitpop[fitidx]

    return false
end

ext: jl


text:
import Interpolations:
        LinearInterpolation,
        Line

"""
    linear_interpolation(x, y; extrapolation_bc = NaN, sort = true)

Create a linear interpolation object from two vectors `x` and `y` (to be interpolated).
"""
function linear_interpolation(x, y; extrapolation_bc = NaN, sort = true)
    # Knots must be sorted in ascending order.
    if sort
        sort_inds = sortperm(x)
        itp = LinearInterpolation(x[sort_inds], y[sort_inds], extrapolation_bc = extrapolation_bc)
    else
        itp = LinearInterpolation(x, y, extrapolation_bc = extrapolation_bc)
    end

    # Don't evaluate the interpolation object on an x-axis grid we provide, but return
    # the interpolation object that can be interpolated at a set of x values later
    # (i.e. itp(inpt_grid_x))
    return itp
end

export linear_interpolation


text:
include("save_result.jl");
include("open_externally.jl");
include("print_log.jl");
include("print_description.jl");
include("load_result.jl");
include("rename.jl");
include("print_all_descriptions.jl");
include("find_results.jl");
include("amend_description.jl");
include("res_dir_tree.jl");

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
"""
    Chain(rules...)
    Chain(rules::Tuple)

`Chain`s allow chaining rules together to be completed in a single processing step,
without intermediate reads or writes from grids.

They are potentially compiled together into a single function call, especially if you
use `@inline` on all `applyrule` methods. `Chain` can hold either all [`CellRule`](@ref)
or [`NeighborhoodRule`](@ref) followed by [`CellRule`](@ref).

[`SetRule`](@ref) can't be used in `Chain`, as it doesn't have a return value.

![Chain rule diagram](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/Chain.png)
"""
struct Chain{R,W,T<:Union{Tuple{},Tuple{Union{<:NeighborhoodRule,<:CellRule},Vararg{<:CellRule}}}} <: RuleWrapper{R,W}
    rules::T
end
Chain(rules...) = Chain(rules)
Chain(rules::Tuple) = begin
    rkeys = Tuple{union(map(k -> _asiterable(_readkeys(k)), rules)...)...}
    wkeys = Tuple{union(map(k -> _asiterable(_writekeys(k)), rules)...)...}
    Chain{rkeys,wkeys,typeof(rules)}(rules)
end

# Getter
rules(chain::Chain) = chain.rules
# Only the first rule in a chain can be a NeighborhoodRule
radius(chain::Chain) = radius(chain[1])
neighborhoodkey(chain::Chain) = neighborhoodkey(chain[1])
neighborhood(chain::Chain) = neighborhood(chain[1])
neighbors(chain::Chain) = neighbors(chain[1])
@inline function setwindow(chain::Chain{R,W}, win) where {R,W}
    rules = (setwindow(chain[1], win), tail(chain.rules)...)
    Chain{R,W,typeof(rules)}(rules)
end

function Base.tail(chain::Chain{R,W}) where {R,W}
    chaintail = tail(rules(chain))
    Chain{R,W,typeof(chaintail)}(chaintail)
end
Base.getindex(chain::Chain, i) = getindex(rules(chain), i)
Base.iterate(chain::Chain) = iterate(rules(chain))
Base.length(chain::Chain) = length(rules(chain))
Base.firstindex(chain::Chain) = firstindex(rules(chain))
Base.lastindex(chain::Chain) = lastindex(rules(chain))

ruletype(chain::Chain) = ruletype(first(chain))

@generated function applyrule(data::AbstractSimData, chain1::Chain{R,W,T}, state1, index) where {R,W,T}
    expr = Expr(:block)
    nrules = length(T.parameters)
    for i in 1:nrules
        # Variables are numbered to make debugging type stability easier
        state = Symbol("state$i")
        nextstate = Symbol("state$(i+1)")
        rule = Symbol("rule$i")
        read = Symbol("read$i")
        write = Symbol("write$i")
        chain = Symbol("chain$i")
        nextchain = Symbol("chain$(i+1)")
        rule_expr = quote
            $rule = $chain[1]
            # Get the state needed by this rule
            $read = _filter_readstate($rule, $state)
            # Run the rule
            $write = applyrule(data, $rule, $read, index)
            # Create new state with the result and state from other rules
            $nextstate = _update_chainstate($rule, $state, $write)
            $nextchain = tail($chain)
        end
        push!(expr.args, rule_expr)
    end
    laststate = Symbol("state$(nrules+1)")
    push!(expr.args, :(_filter_writestate(chain1, $laststate)))
    expr
end

function modifyrule(chain::Chain{R,W}, simdata) where {R,W}
    Chain{R,W}(_modifyrules(rules(chain), simdata))
end

# Get the state to pass to the specific rule as a `NamedTuple` or single value
@generated function _filter_readstate(::Rule{R,W}, state::NamedTuple) where {R<:Tuple,W}
    expr = Expr(:tuple)
    keys = Tuple(R.parameters)
    for k in keys
        push!(expr.args, :(state[$(QuoteNode(k))]))
    end
    :(NamedTuple{$keys}($expr))
end
@inline _filter_readstate(::Rule{R,W}, state::NamedTuple) where {R,W} = state[R]

# Get the state to write for the specific rule
@generated function _filter_writestate(::Rule{R,W}, state::NamedTuple) where {R<:Tuple,W<:Tuple}
    expr = Expr(:tuple)
    keys = Tuple(W.parameters)
    for k in keys
        push!(expr.args, :(state[$(QuoteNode(k))]))
    end
    expr
end

# Merge new state with previous state.
# Returns a new `NamedTuple` with all keys having the most recent state
@generated function _update_chainstate(rule::Rule{R,W}, state::NamedTuple{K,V}, writestate
                                     ) where {R,W,K,V}
    expr = Expr(:tuple)
    writekeys = W isa Symbol ? (W,) : W.parameters
    keys = (union(K, writekeys)...,)
    for (i, k) in enumerate(keys)
        if k in writekeys
            for (j, wkey) in enumerate(writekeys)
                if k == wkey
                    push!(expr.args, :(writestate[$j]))
                end
            end
        else
            push!(expr.args, :(state[$i]))
        end
    end
    :(NamedTuple{$keys}($expr))
end


text:
using PolyChaos, Test, QuadGK
import LinearAlgebra: norm

names = ["gaussian", "uniform01", "logistic", "legendre", "hermite", "laguerre"]
deg = 10
tol = 1e-7

@testset "Norms of orthgonal polynomials ($names)" begin
    for name in names
        # display(name)
        opq = OrthoPolyQ(name,deg)
        s = computeSP2(opq)
        t = [ quadgk(x->evaluate(d,x,opq)^2*opq.op.meas.w(x), opq.op.meas.dom...; order = max(10,10d))[1] for d=0:deg ]
        u = computeSP2(deg,opq.op)
        # display([s t])
        @test isapprox(norm(s-t,Inf)/maximum(s), 0; atol = tol)
        @test isapprox(norm(s-u,Inf)/maximum(s), 0; atol = tol)
        @test isapprox(computeSP([0, 0, 0], opq), t[1]; atol = tol)
        @test isapprox(computeSP([0, 0, 0, 1], opq), 0; atol = tol)
        @test isapprox(computeSP([0, 0, 0, deg-1, deg-1], opq), s[deg]; atol = tol)
    end
end

a, b = 1.1:0.5:4, 1.1:0.5:4
names = ["beta01", "jacobi", "genhermite", "genlaguerre"]

@testset "Norms of orthgonal polynomials ($names)" begin
    for α in a, β in b, name in names
        if name == "beta01" || name == "jacobi"
            pars = Dict(:shape_a=>α,:shape_b=>β)
        elseif name == "genhermite"
            pars = Dict(:mu=>α)
        elseif name == "genlaguerre"
            pars = Dict(:shape=>α)
        end
        opq = OrthoPolyQ(name,deg,pars)
        s = computeSP2(opq)
        t = [ quadgk(x->evaluate(d,x,opq)^2*opq.op.meas.w(x), opq.op.meas.dom...; order = max(10,10d))[1] for d=0:deg ]
        u = computeSP2(deg,opq.op)
        @test isapprox(norm(s-t,Inf)/maximum(s), 0; atol = tol)
        @test isapprox(norm(s-u,Inf)/maximum(s), 0; atol = tol)
        @test isapprox(computeSP([0, 0, 0], opq), t[1]; atol = tol)
        @test isapprox(computeSP([0, 0, 0, 1], opq), 0; atol = tol)
        @test isapprox(computeSP([0, 0, 0, deg-1, deg-1], opq), s[deg]; atol = tol)
    end
end


text:
"Movable ball in N dimensions"
mutable struct MovableBall{N,T}
    x::SVector{N,T}  # position of centre
    v::SVector{N,T}  # velocity
    r::T             # radius
    m::T             # mass
end

MovableBall(x, v, r) = MovableBall(x, v, r, one(r))

centre(b::MovableBall) = b.x

"Normal vector at point x on sphere"
normal(b::MovableBall, x) = normalize(x - centre(b))

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
# Temporary definitions until JuliaML stuff matures

abstract Algorithm

abstract Model
abstract BivariateModel <: Model  # LogisticRegression and SVMLike
# function Base.show(io::IO, m::Model)
#     s = string(typeof(m))
#     s = replace(s, "SparseRegression.", "")
#     print(s)
# end

# ====================================================================# Model methods
function lossvector!(m::Model, storage::Vector, y::Vector, η::Vector)
    for i in eachindex(y)
        @inbounds storage[i] = loss(m, y[i], η[i])
    end
end
function loss(m::Model, y::Vector, η::Vector)
    storage = zeros(length(y))
    lossvector!(m, storage, y, η)
    mean(storage)
end
function predict!(m::Model, storage::Vector, η::Vector)
    for i in eachindex(η)
        @inbounds storage[i] = predict(m, η[i])
    end
end
function predict(m::Model, η::Vector)
    storage = zeros(length(η))
    predict!(m, storage, η)
    storage
end
function classify!(m::BivariateModel, storage::Vector, η::Vector)
    for i in eachindex(η)
        @inbounds storage[i] = classify(m, η[i])
    end
end
function classify(m::BivariateModel, η::Vector)
    storage = zeros(length(η))
    classify!(m, storage, η)
    storage
end

#------------------------------------------------------------------# LinearRegression
immutable LinearRegression <: Model end
loss(m::LinearRegression, y::Real, η::Real) = 0.5 * (y - η) ^ 2
lossderiv(m::LinearRegression, y::Real, η::Real) = -(y - η)
predict(m::LinearRegression, η::Real) = η

#----------------------------------------------------------------------# L1Regression
immutable L1Regression <: Model end
loss(m::L1Regression, y::Real, η::Real) = abs(y - η)
lossderiv(m::L1Regression, y::Real, η::Real) = -sign(y - η)
predict(m::L1Regression, η::Real) = η

#----------------------------------------------------------------# LogisticRegression
"For data in {0, 1}"
immutable LogisticRegression <: BivariateModel end
loss(m::LogisticRegression, y::Real, η::Real) = -y * η + log(1.0 + exp(η))
lossderiv(m::LogisticRegression, y::Real, η::Real) = -(y - predict(m, η))
predict(m::LogisticRegression, η::Real) = 1.0 / (1.0 + exp(-η))
classify(m::LogisticRegression, η::Real) = Float64(η > 0.0)

#------------------------------------------------------------------# ProbitRegression
# TODO
# "For data in {0, 1}"
# immutable ProbitRegression <: BivariateModel end
# d = Ds.Normal()
# function loss(m::ProbitRegression, y::Real, η::Real)
#     -y * Ds.logcdf(d, η) - (1 - y) * Ds.logccdf(d, η)
# end
# loglikelihood(m::ProbitRegression, y::Real, η::Real) = -loss(m, y, η)
# function lossderiv(m::ProbitRegression, y::Real, η::Real)
#     (y - Ds.logcdf(d, η)) * Ds.pdf(η) / (Ds.cdf(d, η) * Ds.ccdf(d, η))
# end
# predict(m::ProbitRegression, η::Real) = Ds.cdf(d, η)
# classify(m::ProbitRegression, η::Real) = Float64(η > 0.0)

#----------------------------------------------------------------# PoissonRegression
immutable PoissonRegression <: Model end
loss(m::PoissonRegression, y::Real, η::Real) = -y * η + exp(η)
lossderiv(m::PoissonRegression, y::Real, η::Real) = -y + exp(η)
predict(m::PoissonRegression, η::Real) = exp(η)

#---------------------------------------------------------------------------# SVMLike
"For data in {-1, 1}"
immutable SVMLike <: BivariateModel end
loss(m::SVMLike, y::Real, η::Real) = max(0.0, 1.0 - y * η)
lossderiv(m::SVMLike, y::Real, η::Real) = 1.0 < y * η ? 0.0 : -y
predict(m::SVMLike, η::Real) = η
classify(m::SVMLike, η::Real) = sign(η)

#----------------------------------------------------------------# QuantileRegression
immutable QuantileRegression <: Model τ::Float64 end
function loss(m::QuantileRegression, y::Real, η::Real)
    r = y - η
    r * (m.τ - (r < 0.0))
end
lossderiv(m::QuantileRegression, y::Real, η::Real) = (y < η) - m.τ
predict(m::QuantileRegression, η::Real) = η

#-------------------------------------------------------------------# HuberRegression
immutable HuberRegression <: Model δ::Float64 end
function loss(m::HuberRegression, y::Real, η::Real)
    r = y - η
    abs(r) < m.δ ? 0.5 * r * r : m.δ * (abs(r) - 0.5 * m.δ)
end
function lossderiv(m::HuberRegression, y::Real, η::Real)
    r = y - η
    abs(r) <= m.δ ? -r : m.δ * sign(-r)
end
predict(m::HuberRegression, η::Real) = η

# ========================================================================= # Penalty
abstract Penalty
#------------------------------------------------------------------# abstract methods
function prox!(p::Penalty, β::VecF, step::Float64)
    for j in eachindex(β)
        @inbounds β[j] = prox(p, β[j], step)
    end
end
# if step different for each element
function prox!(p::Penalty, β::VecF, step::VecF)
    for j in eachindex(β)
        @inbounds β[j] = prox(p, β[j], step[j])
    end
end
function value(p::Penalty, β::VecF)
    v = zeros(β)
    for j in eachindex(v)
        @inbounds v[j] = value(p, β[j])
    end
    v
end
#-------------------------------------------------------------------------# NoPenalty
immutable NoPenalty <: Penalty end
value(p::NoPenalty, β::Float64) = 0.0
deriv(p::NoPenalty, βj::Float64) = 0.0
prox(p::NoPenalty, βj::Float64, s::Float64) = βj
#----------------------------------------------------------------------# RidgePenalty
immutable RidgePenalty <: Penalty
    λ::Float64
    RidgePenalty(λ::Real) = (@assert λ >= 0; new(λ))
end
value(p::RidgePenalty, β::VecF) = 0.5 * p.λ * sumabs2(β)
deriv(p::RidgePenalty, βj::Float64) = p.λ * βj
prox(p::RidgePenalty, βj::Float64, s::Float64) = βj / (1.0 + s * p.λ)
#----------------------------------------------------------------------# LassoPenalty
immutable LassoPenalty <: Penalty
    λ::Float64
    LassoPenalty(λ::Real) = (@assert λ >= 0; new(λ))
end
value(p::LassoPenalty, β::VecF) = p.λ * sumabs(β)
prox(p::LassoPenalty, βj::Float64, s::Float64) = sign(βj) * max(abs(βj) - s * p.λ, 0.0)
deriv(p::LassoPenalty, βj::Float64) = p.λ * sign(βj)
#-----------------------------------------------------------------# ElasticNetPenalty
immutable ElasticNetPenalty <: Penalty
    λ::Float64
    a::Float64  # proportion of Lasso
    function ElasticNetPenalty(λ::Real, a::Real = .9)
        @assert λ >= 0
        @assert 0 <= a <= 1
        new(λ, a)
    end
end
function value(p::ElasticNetPenalty, β::VecF)
    p.λ * (p.a * sumabs(β) + (1. - p.a) * 0.5 * sumabs2(β))
end
function prox(p::ElasticNetPenalty, βj::Float64, s::Float64)
    sign(βj) * max(abs(βj) - s * p.λ * p.a, 0.0) / (1.0 + s * p.λ * (1.0 - p.a))
end
function deriv(p::ElasticNetPenalty, βj::Float64)
    p.a * deriv(LassoPenalty(p.λ), βj) + (1 - p.a) * deriv(RidgePenalty(p.λ), βj)
end

#----------------------------------------------------------------------# SCADPenalty
# TODO
# immutable SCADPenalty <: Penalty
#     a::Float64
#     function SCADPenalty(a::Real = 3.7)
#         @assert a > 2
#         new(a)
#     end
# end
# Base.show(io::IO, p::SCADPenalty) = print(io, "SCADPenalty (a = $(p.a))")
# function penalty(p::SCADPenalty, β::VecF, λ::Float64)
#     val = 0.0
#     for j in eachindex(β)
#         βj = abs(β[j])
#         if βj < λ
#             val += λ * βj
#         elseif βj < λ * p.a
#             val -= 0.5 * (βj ^ 2 - 2.0 * p.a * λ * βj + λ ^ 2) / (p.a - 1.0)
#         else
#             val += 0.5 * (p.a + 1.0) * λ ^ 2
#         end
#     end
#     return val
# end
# function prox(p::SCADPenalty, βj::Float64, c::Float64)
#     if abs(βj) > p.a * c
#     elseif abs(βj) < 2.0 * c
#         βj = sign(βj) * max(abs(βj) - c, 0.0)
#     else
#         βj = (βj - c * sign(βj) * p.a / (p.a - 1.0)) / (1.0 - (1.0 / p.a - 1.0))
#     end
#     βj
# end

#--------------------------------------------------------------------# sweep! methods
"""
`sweep!(A, k, inv = false)`, `sweep!(A, k, v, inv = false)`

Symmetric sweep operator of the matrix `A` on element `k`.  `A` is overwritten.
`inv = true` will perform the inverse sweep.  Only the upper triangle is read and swept.

An optional vector `v` can be provided to avoid memory allocation.
This requires `length(v) == size(A, 1)`.  Both `A` and `v`
will be overwritten.

```julia
x = randn(100, 10)
xtx = x'x
sweep!(xtx, 1)
sweep!(xtx, 1, true)
```
"""
function sweep!{T<:Real}(A::AMat{T}, k::Integer, inv::Bool = false)
    n, p = size(A)
    # ensure @inbounds is safe
    @assert n == p "A must be square"
    @assert k <= p "pivot element not within range"
    @inbounds d = 1.0 / A[k, k]  # pivot
    # get column A[:, k] (hack because only upper triangle is available)
    akk = zeros(p)
    for j in 1:p
        if j <= k
            @inbounds akk[j] = A[j, k]
        else
            @inbounds akk[j] = A[k, j]
        end
    end
    BLAS.syrk!('U', 'N', -d, akk, 1.0, A)  # everything not in col/row k
    scale!(akk, d * (-1) ^ inv)
    for i in 1:k-1  # col k
        @inbounds A[i, k] = akk[i]
    end
    for j in k+1:p  # row k
        @inbounds A[k, j] = akk[j]
    end
    A[k, k] = -d  # pivot element
    A
end

function sweep!{T<:Real, I<:Integer}(A::AMat{T}, ks::AVec{I}, inv::Bool = false)
    for k in ks
        sweep!(A, k, inv)
    end
    A
end

function sweep!{T<:Real}(A::AMat{T}, k::Integer, v::AVecF, inv::Bool = false)
    n, p = size(A)
    # ensure that @inbounds is safe
    @assert n == p "A must be square"
    @assert length(v) == p "storage length ≠ size(A, 1)"
    @assert k <= p "pivot element not within range"
    @inbounds d = 1.0 / A[k, k]  # pivot
    for j in 1:p   # get column A[:, k]
        if j <= k
            @inbounds v[j] = A[j, k]
        else
            @inbounds v[j] = A[k, j]
        end
    end
    BLAS.syrk!('U', 'N', -d, v, 1.0, A)  # everything not in col/row k
    scale!(v, d * (-1) ^ inv)
    for i in 1:k-1  # col k
        @inbounds A[i, k] = v[i]
    end
    for j in k+1:p  # row k
        @inbounds A[k, j] = v[j]
    end
    @inbounds A[k, k] = -d  # pivot element
    A
end

function sweep!{T<:Real,I<:Integer}(A::AMat{T}, ks::AVec{I}, v::VecF, inv::Bool = false)
    for k in ks
        sweep!(A, k, v, inv)
    end
    A
end


text:
using LibExpat
using Base.Test

pd = xp_parse(open(readall, joinpath(splitdir(@__FILE__)[1],"t_s1.txt")))
@test isa(pd, ETree)
println("PASSED 1")

ret = find(pd, "/ListBucketResult")
@test isa(ret, Array)
@test length(ret) == 1
@test isa(ret[1], ETree)
println("PASSED 1.1")

ret = find(pd, "/ListBucketResult/Name")
@test isa(ret, Array)
println("PASSED 2")

ret = find(pd, "/ListBucketResult/Name#string")
@test ret == "bucket"
println("PASSED 3")

ret = find(pd, "/ListBucketResult/Contents")
@test isa(ret, Array)
@test length(ret) == 2
@test isa(ret[1], ETree)
@test isa(ret[2], ETree)
println("PASSED 4")

@test_throws ErrorException find(pd, "/ListBucketResult/Contents#string")
println("PASSED 5")

ret = split(strip(find(pd, "/ListBucketResult/Contents[1]#string")),'\n')
@test ret[1] == "C1C1C1"
println("PASSED 6")

@test (find(pd, "/ListBucketResult/Contents[1]#string") == find(pd, "Contents[1]#string"))
println("PASSED 6.1")

ret = split(strip(find(pd, "/ListBucketResult/Contents[2]#string")),'\n')
@test ret[1] == "C2C2C2"
println("PASSED 7")

ret = find(pd, "/ListBucketResult/Contents[1]/Owner/ID")
@test isa(ret, Array)
@test length(ret) == 1
@test isa(ret[1], ETree)
println("PASSED 8")

ret = find(pd, "/ListBucketResult/Contents[1]/Owner/ID#string")
@test ret == "11111111111111111111111111111111"
println("PASSED 9")

ret = find(pd, "/ListBucketResult/Contents[1]/Owner/ID{idk}")
@test ret == "IDKV1"
println("PASSED 10")

ret = find(pd, "/ListBucketResult/Contents[2]/Owner/ID{idk}")
@test ret == "IDKV2"
println("PASSED 11")

@test (find(pd, "/ListBucketResult/Contents[2]/Owner/ID{idk}") == find(pd, "Contents[2]/Owner/ID{idk}"))
println("PASSED 11.1")

@test (find(pd, "/I/Do/NOT/Exist") == [])
println("PASSED 12")

@test (find(pd, "/I/Do/NOT/Exist[1]") == nothing)
println("PASSED 12.1")

@test (find(pd, "/ListBucketResult/Contents[2]/Owner/JUNK#string") == nothing)
println("PASSED 12.2")

pd = xp_parse(open(readall, joinpath(splitdir(@__FILE__)[1],"utf8.xml")))
@test isa(pd, ParsedData)
println("PASSED 13")

pd = xp_parse(open(readall, joinpath(splitdir(@__FILE__)[1],"wiki.xml")))
@test isa(pd, ParsedData)
ret = find(pd, "/page/revision/id#string")
@test ret == "557462847"
println("PASSED 14")

# Check streaming parse functionality
found_start_element = false
found_end_element = false
start_element_name = ""
first_attr_start_element = ""
first_attrval_start_element = ""
end_element_name = ""

cbs = XPCallbacks()
cbs.start_element = function (h, name, attrs)
    global found_start_element = true
    global start_element_name = name
    global first_attr_start_element = collect(keys(attrs))[1]
    global first_attrval_start_element = attrs[first_attr_start_element]
end

cbs.end_element = function (h, name)
    global found_end_element = true
    global end_element_name = name
end

parse("<test id=\"someid\">somecontent</test>", cbs)

@test found_start_element
println("PASSED 15.1")
@test found_end_element
println("PASSED 15.2")
@test (start_element_name == "test")
println("PASSED 15.3")
@test (first_attr_start_element == "id")
println("PASSED 15.4")
@test (first_attrval_start_element == "someid")
println("PASSED 15.5")
@test (end_element_name == "test")
println("PASSED 15.6")

# Check streaming file functionality with more callbacks
cdata_started = false
cdata_ended = false
cdata_section = ""
comment = ""
in_cdata_section = false

cbs = XPCallbacks()

cbs.start_cdata = function (h)
    global in_cdata_section = true
    global cdata_started = true
end

cbs.end_cdata = function (h)
    global in_cdata_section = false
    global cdata_ended = true
end

cbs.comment = function (h, txt)
    global comment = strip(txt)
end

cbs.character_data = function (h, data)
    if in_cdata_section
        global cdata_section = string(cdata_section, data)
    end
end

parsefile(joinpath(splitdir(@__FILE__)[1], "graph.graphml"), cbs)

@test cdata_started
println("PASSED 16.1")
@test cdata_ended
println("PASSED 16.2")
@test (strip(cdata_section) == "<Some CDATA section>")
println("PASSED 16.3")
@test (comment == "This is just a comment.")
println("PASSED 16.4")


text:
module TestBroadcast
    using NullableArrays
    using Base.Test
    using Compat

    A1 = rand(Int, 10)
    M1 = rand(Bool, 10)
    n = rand(2:5)
    dims = [ rand(2:5) for i in 1:n]
    A2 = rand(10, dims...)
    M2 = rand(Bool, 10, dims...)
    C2 = Array{Float64}(10, dims...)
    i = rand(2:5)
    A3 = rand(10, [dims; i]...)
    M3 = rand(Bool, 10, [dims; i]...)
    C3 = Array{Float64}(10, [dims; i]...)

    m1 = broadcast((x,y)->x, M1, M2)
    m2 = broadcast((x,y)->x, M2, M3)
    R2 = reshape(Bool[ (x->(m1[x] | M2[x]))(i) for i in 1:length(M2) ], size(M2))
    R3 = reshape(Bool[ (x->(m2[x] | M3[x]))(i) for i in 1:length(M3) ], size(M3))
    r2 = broadcast((x,y)->x, R2, R3)
    Q3 = reshape(Bool[ (x->(r2[x] | R3[x]))(i) for i in 1:length(R3) ], size(M3))

    U1 = NullableArray(A1)
    U2 = NullableArray(A2)
    U3 = NullableArray(A3)
    V1 = NullableArray(A1, M1)
    V2 = NullableArray(A2, M2)
    V3 = NullableArray(A3, M3)
    Z2 = NullableArray(Float64, 10, dims...)
    Z3 = NullableArray(Float64, 10, [dims; i]...)

    f() = 5
    f(x::Real) = 5 * x
    f(x::Real, y::Real) = x * y
    f(x::Real, y::Real, z::Real) = x * y * z

    @inferred NullableArrays.lift(f, (1,))
    @inferred NullableArrays.lift(f, (1, 2.0))
    @inferred NullableArrays.lift(f, (1, 2.0, 3.0))

    for (dests, arrays, nullablearrays, mask) in
        ( ((C2, Z2), (A1, A2), (U1, U2), ()),
          ((C3, Z3), (A2, A3), (U2, U3), ()),
          ((C3, Z3), (A1, A2, A3), (U1, U2, U3), ()),

          ((C2, Z2), (A1, A2), (V1, V2), (R2,)),
          ((C3, Z3), (A2, A3), (V2, V3), (R3,)),
          ((C3, Z3), (A1, A2, A3), (V1, V2, V3), (Q3,)),
    )

        # Base.broadcast!(f, B::NullableArray, As::NullableArray...)
        broadcast!(f, dests[1], arrays...)
        broadcast!(f, dests[2], nullablearrays...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))

        # Base.broadcast(f, As::NullableArray...)
        D = broadcast(f, arrays...)
        X = broadcast(f, nullablearrays...)
        @test isequal(X, NullableArray(D, mask...))
    end

    # Base.broadcast!(f, X::NullableArray)
    for (array, nullablearray, mask) in
        ( (A1, U1, ()), (A2, U2, ()), (A3, U3, ()),
          (A1, V1, (M1,)), (A2, V2, (M2,)), (A3, V3, (M3,)),
    )
        broadcast!(f, array)
        broadcast!(f, nullablearray)
        @test isequal(nullablearray, NullableArray(array, mask...))
    end

    # test broadcasted arithmetic operators
    A = rand(10)
    X1 = NullableArray(A)
    n = rand(2:5)
    dims = rand(2:5, n)
    B = rand(Float64, 10, dims...)
    X2 = NullableArray(B)
    M = rand(Bool, 10, dims...)
    Y = NullableArray(B, M)

    if VERSION < v"0.6.0-dev.1632"
        for op in (
            (.+),
            (.-),
            (.*),
            (./),
            (.%),
            (.^),
            (.==),
            (.!=),
            (.<),
            (.>),
            (.<=),
            (.>=),
        )
            @test isequal(op(X1, X2), NullableArray(op(A, B)))
            @test isequal(op(X1, Y), NullableArray(op(A, B), M))
        end
    end

    A = rand(Bool, 100)
    B = rand(Bool, 100)
    M1 = rand(Bool, 100)
    M2 = rand(Bool, 100)
    X = NullableArray(A, M1)
    Y = NullableArray(B, M2)
    @test isequal(broadcast(&, X, Y), NullableArray(A .& B, M1 .| M2))
    @test isequal(broadcast(|, X, Y), NullableArray(A .| B, M1 .| M2))

    # Test broadcasting with constructor
    immutable SurvEvent
        time::Float64
        censored::Bool
    end
    t = NullableArray(rand(3))
    c = NullableArray(rand(Bool, 3))
    @test isequal(SurvEvent.(t, c), NullableArray([SurvEvent(get(t[i]), get(c[i])) for i in 1:3]))
    @test isa(SurvEvent.(t, c), NullableVector{SurvEvent})
    @inferred NullableArrays.lift(SurvEvent, (1, true))

end # module


text:
baremodule CoreTypes # Convenience
using ..SymbolServer
using Base: ==, @static

const Any = SymbolServer.stdlibs[:Core][:Any]
const DataType = SymbolServer.stdlibs[:Core][:DataType]
const Function = SymbolServer.stdlibs[:Core][:Function]
const Module = SymbolServer.stdlibs[:Core][:Module]
const String = SymbolServer.stdlibs[:Core][:String]
const Char = SymbolServer.stdlibs[:Core][:Char]
const Symbol = SymbolServer.stdlibs[:Core][:Symbol]
const Bool = SymbolServer.stdlibs[:Core][:Bool]
const Int = SymbolServer.stdlibs[:Core][:Int]
const UInt8 = SymbolServer.stdlibs[:Core][:UInt8]
const UInt16 = SymbolServer.stdlibs[:Core][:UInt16]
const UInt32 = SymbolServer.stdlibs[:Core][:UInt32]
const UInt64 = SymbolServer.stdlibs[:Core][:UInt64]
const Float64 = SymbolServer.stdlibs[:Core][:Float64]
const Vararg = SymbolServer.FakeTypeName(Core.Vararg)

iscoretype(x, name) = false
iscoretype(x::SymbolServer.VarRef, name) = x isa SymbolServer.DataTypeStore && x.name.name == name && x.name isa SymbolServer.VarRef && x.name.parent.name == :Core
iscoretype(x::SymbolServer.DataTypeStore, name) = x isa SymbolServer.DataTypeStore && x.name.name.name == name && x.name.name isa SymbolServer.VarRef && x.name.name.parent.name == :Core
isdatatype(x) = iscoretype(x, :DataType)
isfunction(x) = iscoretype(x, :Function)
ismodule(x) = iscoretype(x, :Module)
isstring(x) = iscoretype(x, :String)
ischar(x) = iscoretype(x, :Char)
issymbol(x) = iscoretype(x, :Symbol)
@static if Core.Int == Core.Int64
    isint(x) = iscoretype(x, :Int64)
else
    isint(x) = iscoretype(x, :Int32)
end
isfloat(x) = iscoretype(x, :Float64)
isvector(x) = iscoretype(x, :Vector)
isarray(x) = iscoretype(x, :Array)
isva(x::SymbolServer.FakeUnionAll) = isva(x.body)
@static if Core.Vararg isa Core.Type
    function isva(x)
        return (x isa SymbolServer.FakeTypeName && x.name.name == :Vararg &&
            x.name.parent isa SymbolServer.VarRef && x.name.parent.name == :Core)
    end
else
    isva(x) = x isa SymbolServer.FakeTypeofVararg
end
end


text:
module DynamicsModelInterface

export add_dynamics_constraints!, add_dynamics_jacobians!, next_x

function add_dynamics_constraints! end
function add_dynamics_jacobians! end
function next_x end

end


text:
#MDDatasets DataRS (Recursive Sweep) defninitions
#-------------------------------------------------------------------------------

#==Main types
===============================================================================#

#Linked-list representation of multi-dimensional datasets:
#-------------------------------------------------------------------------------
mutable struct DataRS{T} <: DataMD
	sweep::PSweep
	elem::Vector{T}

	function DataRS{T}(sweep::PSweep, elem::Vector{T}) where {T}
		if !elemallowed(DataRS, eltype(elem))
			msg = "Can only create DataRS{T} for T ∈ {DataRS, DataF1, DataFloat, DataInt, DataComplex}"
			throw(ArgumentError(msg))
		elseif length(sweep) != length(elem)
			throw(ArgumentError("sweep length does not match number of elem"))
		end
		return new(sweep, elem)
	end
end

#Shorthand (because default (non-parameterized) constructor was overwritten):
DataRS(sweep::PSweep, elem::Vector{T}) where T = DataRS{T}(sweep, elem)

elemallowed(::Type{DataRS}, t::Type{T}) where T = elemallowed(DataMD, t) #Allow basic types
elemallowed(::Type{DataRS}, ::Type{DataRS}) = true #Also allow recursive structures

#Generate empty DataRS structure:
(::Type{DataRS{T}})(sweep::PSweep) where T = DataRS{T}(sweep, Array{T}(undef, length(sweep)))

#==Type promotions
===============================================================================#
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:DataRS, T2<:Number} = DataRS

#==Accessor functions
===============================================================================#
Base.eltype(d::DataRS{T}) where T = T
Base.length(d::DataRS) = length(d.elem)

#==Help with construction
===============================================================================#

#Implement "fill(DataRS, ...) do sweepval" syntax:
function Base.fill!(fn::Function, d::DataRS)
	for i in 1:length(d.sweep)
		d.elem[i] = fn(d.sweep.v[i])
	end
	return d
end
Base.fill(fn::Function, ::Type{DataRS{T}}, sweep::PSweep) where T =
	fill!(fn, DataRS(sweep, Array{T}(undef, length(sweep))))
Base.fill(fn::Function, ::Type{DataRS}, sweep::PSweep) = fill(fn, DataRS{DataRS}, sweep)

#==Data generation
===============================================================================#
function _ensuresweepunique(d::DataRS, sweepid::String)
	if sweepid == d.sweep.id
		msg = "Sweep occurs multiple times in DataRS: $sweepid"
		throw(ArgumentError(msg))
	end
end

#Define "parameter".
#(Generates a DataRS object containing the value of a given swept parameter)
#-------------------------------------------------------------------------------

#Deal with non-leaf elements, once the sweep value is found:
function _parameter(d::DataRS{DataRS}, sweepid::String, sweepval::T) where T
	_ensuresweepunique(d, sweepid)
	elem = DataRS[_parameter(d.elem[i], sweepid, sweepval) for i in 1:length(d.sweep)]
	return DataRS(d.sweep, elem)
end

#Deal with leaf elements, once the sweep value is found:
function _parameter(d::DataRS, sweepid::String, sweepval::T) where T
	_ensuresweepunique(d, sweepid)
	elem = T[sweepval for i in 1:length(d.sweep)]
	return DataRS(d.sweep, elem)
end

#Main "parameter" algorithm (non-leaf elements):
function parameter(d::DataRS{DataRS}, sweepid::String)
	if sweepid == d.sweep.id #Sweep found
		elem = DataRS[_parameter(d.elem[i], sweepid, d.sweep.v[i]) for i in 1:length(d.sweep)]
	else
		elem = DataRS[parameter(d.elem[i], sweepid) for i in 1:length(d.sweep)]
	end
	return DataRS(d.sweep, elem)
end
#Main "parameter" algorithm (leaf elements):
function parameter(d::DataRS, sweepid::String)
	T = eltype(d.sweep.v)
	if sweepid == d.sweep.id #Sweep found
		return DataRS(d.sweep, d.sweep.v)
	else
		msg = "Sweep not found in DataRS: $sweepid"
		throw(ArgumentError(msg))
	end
end

#Generate DataRS from DataHR.
#-------------------------------------------------------------------------------
function _buildDataRS(d::DataHR, firstinds::Vector{Int})
	curidx = length(firstinds) + 1
	sweep = d.sweeps[curidx]
	if curidx < length(d.sweeps)
		result = DataRS{DataRS}(sweep)
		for i in 1:length(sweep.v)
			result.elem[i] = _buildDataRS(d, vcat(firstinds, i))
		end
	else #Last index.  Copy data over:
		result = DataRS{eltype(d.elem)}(sweep)
		for i in 1:length(sweep.v)
			result.elem[i] = d.elem[firstinds..., i]
		end
	end
	return result
end

function DataRS(d::DataHR)
	return _buildDataRS(d, Int[])
end

#==User-friendly show functions
===============================================================================#

#Print leaf element:
function printDataRSelem(io::IO, ds::DataRS, idx::Int, indent::String)
	if isassigned(ds.elem, idx)
		println(io, ds.elem[idx])
	else
		println(io, indent, "UNDEFINED")
	end
end
#Print next level of recursive DataRS:
function printDataRSelem(io::IO, ds::DataRS{T}, idx::Int, indent::String) where T<:DataRS
	println(io)
	if isassigned(ds.elem, idx)
		printDataRS(io, ds.elem[idx], indent)
	else
		println(io, indent, "UNDEFINED")
	end
end
#Print DataRS structure:
function printDataRS(io::IO, ds::DataRS, indent::String)
	for i in 1:length(ds.elem)
		print(io, "$indent", ds.sweep.id, "=", ds.sweep.v[i], ": ")
		printDataRSelem(io, ds, i, "$indent  ")
	end
end

function Base.show(io::IO, ds::DataRS)
	print(io, "DataRS[\n")
	printDataRS(io, ds, "  ")
	print(io, "]\n")
end


text:
include("fourierbasis.jl")

struct ConcatentateBasis{TB} <: Any
    ϕs::TB
end

struct ConcatentateBasisBuffer{TO} <: Any
    out::TO

    function ConcatentateBasisBuffer(ϕ::ConcatentateBasis)
        n = sum(length(ϕ.ϕs))
        out = zeros(n)
        return new{typeof(out)}(out)
    end
end

function (ϕ::ConcatentateBasis)(x)
    n = sum(length(ϕ.ϕs))
    out = zeros(n)
    tot = 1
    for i in 1:length(ϕ.ϕs)
        ni = length(ϕ.ϕs[i])
        out[tot:tot+ni-1] .= ϕ.ϕs[i](x)
        tot += ni
    end
    return out
end

function (ϕ::ConcatentateBasis)(buff::ConcatentateBasisBuffer, x)
    out = buff.out
    tot = 1
    for i in 1:length(ϕ.ϕs)
        ni = length(ϕ.ϕs[i])
        out[tot:tot+ni-1] .= ϕ.ϕs[i](buff.buffs[i], x)
        tot += ni
    end
    return out
end

function (ϕ::ConcatentateBasis)(buff, x)
    out = buff
    tot = 1
    for i in 1:length(ϕ.ϕs)
        y = ϕ.ϕs[i](x)
        ni = length(y)
        out[tot:tot+ni-1] .= y
        tot += ni
    end
    return out
end


text:
using Pkg;
Pkg.activate("");
using CUTEst, Dates, NLPModels, SolverBenchmark

using ARCTR, JSOSolvers, NLPModelsIpopt

nmax = 10
problems = readlines("list_problems_$nmax.dat")
cutest_problems = (CUTEstModel(p) for p in problems)

max_time = 3600.0
max_ev = typemax(Int)
max_iter = typemax(Int)
atol = 1e-5
rtol = 1e-6

solvers = Dict(
    :ipopt =>
        nlp -> ipopt(
            nlp,
            print_level = 0,
            dual_inf_tol = Inf,
            constr_viol_tol = Inf,
            compl_inf_tol = Inf,
            acceptable_iter = 0,
            max_cpu_time = max_time,
            tol = rtol,
        ),
    :lbfgs =>
        nlp -> lbfgs(
            nlp,
            verbose = 0,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_eval = max_ev,
        ),
    :trunk =>
        nlp -> trunk(
            nlp,
            verbose = 0,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_eval = max_ev,
        ),
    :tron =>
        nlp -> tron(
            nlp,
            verbose = 0,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_eval = max_ev,
        ),
)
stats = bmark_solvers(solvers, cutest_problems)

using JLD2
@save "$(today())_$(prod(String.(keys(solvers)) .* "_"))cutest_$(string(length(problems)))_$(nmax).jld2" stats


text:
##
#
# Polynomial iterators
#

"""

    PolynomialIterator(coefficients, degree)

Iterate over all polynomials with coefficients provided by array `coefficients`
and maximal degree `degree`.

Polynomials are represented by arrays of their coefficients from the lowest to
the highest order (as in the
[PolynomialRoots](https://github.com/giordano/PolynomialRoots.jl) package).
For example, the polynomial ``z^2 + 2z + 3`` corresponds to
the array `[3, 2, 1]`.

The constructor removes duplicate values from the `coefficients` array.
"""
struct PolynomialIterator{T <: Number}
  coefficients::Array{T, 1}   # array of coefficients
  degree::Int                 # maximal degree
  number_of_coefficients::Int # number of coefficients

  function PolynomialIterator(coefficients::Array{T, 1}, degree::Int) where {T <: Number}
    coeffs = unique(coefficients)

    # reject nonsensical input
    if length(coeffs) == 0 || coeffs == [T(0)]
      throw(
        ArgumentError("You have to provide at least one nonzero coefficient!")
      )
    end

    # push zero to the beginning
    if T(0) in coeffs
      filter!(x -> x != T(0), coeffs)
      prepend!(coeffs, T(0))
    end

    new{T}(coeffs, degree, length(coeffs))
  end
end

#
# Iterator interface
#
# see https://docs.julialang.org/en/v1/manual/interfaces/
#

"""

  iterate(iter)

Initialize the iterator by returning tuple consisting of a nonzero constant
polynomial and the iterator inner state.
"""
function Base.iterate(iter::PolynomialIterator{T}) where {T <: Number}
  if iter.coefficients[1] == T(0)
    return ([iter.coefficients[2]], [2])
  else
    return ([iter.coefficients[1]], [1])
  end

  return
end

"""

  iterate(iter, state)

Return the next polynomial.
"""
function Base.iterate(iter::PolynomialIterator{T}, state::Array{Int, 1}) where {T <: Number}
  for j in 1:length(state)
    if state[j] < iter.number_of_coefficients
      state[j] += 1
      return ([iter.coefficients[j] for j in state], state)
    else
      state[j] = 1
    end
  end

  # increase the degree
  if length(state) == iter.degree + 1
    return nothing
  else
    if iter.coefficients[1] == T(0) # skip zero
      append!(state, 2)
    else
      append!(state, 1)
    end
    return ([iter.coefficients[j] for j in state], state)
  end
end

"""

    length(iter)

Total number of polynomials specified by the iterator `iter`.
"""
function Base.length(iter::PolynomialIterator{T}) where {T <: Number}
  noc = iter.number_of_coefficients
  if T(0) in iter.coefficients
    return noc ^ (iter.degree + 1) - 1
  elseif noc == 1
    iter.degree + 1
  else
    return noc * div(noc ^ (iter.degree + 1) - 1, noc - 1)
  end
end


text:
"""
    discrete_negloglik::Real

Compute the negative log-likelihood of the process in argument.
"""
function discrete_negloglik(ts::TS;  μ::Real, ϕ::Real, γ::Real)::Real        # one method for point process with or without marks (model without marks)
    tst = [μ, ϕ, γ] .< 0
    any(tst) && (return Inf)

    times = first(ts.times) isa TimeType ? Dates.value.(Date.(ts.times)) : ts.times
    starttime = first(times)
    endtime = last(times)
    anytimes = starttime:oneunit(starttime-endtime):endtime
    disc_ts = TimeSeries(times)

    vol = volfunc(anytimes, disc_ts, γ)     # ν function in Li2020
    intens = μ .+ ϕ .* vol       # rate λ in Li2020
    prob = 1 .- exp.(-intens)       # probability for an event to happen
    t_idx = findall(in(times), anytimes)
    prob_1 = prob[t_idx]        # probability of the events that happened to happen
    prob_0 = 1 .- prob[findall(!in(times), anytimes)] # probability of the events that didn't happen not to happen
    term1 = sum(log.(prob_1)) + sum(log.(prob_0))

    return (-term1)
end

function discrete_negloglik(sepp::SEPPExpKern)
    ts = sepp.data
    isnothing(ts) && error("No data in model, can't compute log-likelihood")

    θ = params(sepp)
    return discrete_negloglik(ts ; θ...)
end

function discrete_negloglik(mts::MarkedTimeSeries, markdens::SupportedMarksDistributions, impact_func::Function = (d -> 1 + tanh(d)) ;  μ::Real, ϕ::Real, γ::Real, δ::Real = 0, ξ::Real, α::Real, β::Real, κ::Real = 1)
    tst = [μ, ϕ, γ, δ, β, α, κ] .< 0
    any(tst) && (return Inf)

    times = first(mts.times) isa TimeType ? Dates.value.(Date.(mts.times)) : mts.times
    starttime = first(times)
    endtime = last(times)
    anytimes= starttime:oneunit(starttime-endtime):endtime
    disc_mts = MarkedTimeSeries(times, mts.marks)

    vol = volfunc(anytimes, disc_mts, γ, δ, impact_func)     # ν function in Li2020
    intens = μ .+ ϕ .* vol       # rate λ in Li2020
    prob = 1 .- exp.(-intens)       # probability for an event to happen
    t_idx = findall(in(times), anytimes)
    prob_1 = prob[t_idx]        # probability of the events that happened to happen
    prob_0 = 1 .- prob[findall(!in(times), anytimes)] # probability of the events that didn't happen not to happen
    term1 = sum(log.(prob_1)) + sum(log.(prob_0))

    σ = β .+ α .* vol[t_idx]
    marks = mts.marks

    sig_marks = hcat(σ, marks)

    if markdens == Distributions.GeneralizedPareto
        mark_contrib = (sig_mark -> logpdf(markdens(0, sig_mark[1], ξ), sig_mark[2])).(eachrow(sig_marks))
    else        # EGPpower case
        mark_contrib = (sig_mark -> logpdf(markdens(sig_mark[1], ξ, κ), sig_mark[2])).(eachrow(sig_marks))
    end

    term2 = sum(mark_contrib)

    return -term1 - term2
end

function discrete_negloglik(sempp::SEMPPExpKern)
    mts = sempp.data
    isnothing(mts) && error("No data in model, can't compute log-likelihood")

    θ = params(sempp)
    markdens = θ[:markdens]
    impact_func = θ[:impact_function]
    delete!(θ, :markdens)
    delete!(θ, :impact_function)
    return discrete_negloglik(mts, markdens, impact_func ; θ...)
end

"""
    discrete_fit!(sepp::SEPP)

Fit a discrete self exciting point process model (whithout marks) to the time series in data, based on MLE.

Note that if the process has marks, this method ignores them, that is modelling the ground process as independent of the marks.
Returns the minimal Log-likelihood found.
"""
function discrete_fit! end

include("macro_def_discrete_fits.jl")

@def_discrete_fit_sepp [:s1]
@def_discrete_fit_sepp [:s1, :s2]
@def_discrete_fit_sepp [:s1, :s2, :s3]

function discrete_fit!(sepp::SEPPExpKern)
    return discrete_fit!(sepp, :μ, :ϕ, :γ)
end

@def_discrete_fit_sempp [:s1]
@def_discrete_fit_sempp [:s1, :s2]
@def_discrete_fit_sempp [:s1, :s2, :s3]
@def_discrete_fit_sempp [:s1, :s2, :s3, :s4]
@def_discrete_fit_sempp [:s1, :s2, :s3, :s4, :s5]
@def_discrete_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6]
@def_discrete_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6, :s7]
@def_discrete_fit_sempp [:s1, :s2, :s3, :s4, :s5, :s6, :s7, :s8]

function discrete_fit!(sempp::SEMPPExpKern)
    return discrete_fit!(sempp, :μ, :ϕ, :γ, :δ, :ξ, :α, :β)
end

#=
function discrete_fit!(sepp::SEPP, params_to_fit::Symbol...) # generic method either to fit a ts whithout marks or to fit the ground process of an mts
    ts=sepp.data
    isnothing(ts) && error("No data in model, can't fit")

    model = Model(Ipopt.Optimizer)
    θ = params(sepp)

    function to_min(μ, ϕ, γ)
        return discrete_negloglik(ts, μ = μ, ϕ = ϕ, γ = γ)
    end

    JuMP.register(model, :to_min, 3, to_min, autodiff=true)

    @variables(model, begin
        mu >= 0, (start = θ[:μ])
        phi >= 0, (start = θ[:ϕ])
        gamma >= 0, (start = θ[:γ])
    end)
    @NLobjective(model, Min, to_min(mu, phi, gamma))

    optimize!(model)

    sepp.μ = value(mu)
    sepp.ϕ = value(phi)
    sepp.γ = value(gamma)

    x = [sepp.μ, sepp.ϕ, sepp.γ]

    to_hess(θ) = to_min(θ...)

    hess = ForwardDiff.hessian(to_hess, x)

    sepp.cov_mat = inv(hess)

    return objective_value(model)
end

"""
    discrete_fit!(sempp::SEMPPExpKern, bounds::Union{Vector{<:Real}, Nothing} = nothing)

Fit a discrete self exciting marked point process model with marks distribution either GPD or EGPpower. Based on MLE.
"""
function discrete_fit!(sempp::SEMPPExpKern, bounds::Union{Vector{<:Real}, Nothing} = nothing) # default xi >= 0
    mts = sempp.data
    isnothing(mts) && error("No data in model, can't fit")

    model = Model(Ipopt.Optimizer)
    θ = params(sempp)
    markdens = θ[:markdens]

    function to_min_GPD(μ, ϕ, γ, δ, ξ, α, β)
        return discrete_negloglik(mts, Distributions.GeneralizedPareto, μ = μ, ϕ = ϕ, γ = γ, δ = δ, ξ = ξ, α = α, β = β)
    end

    function to_min_EGPD(μ, ϕ, γ, δ, ξ, α, β, κ)
        return discrete_negloglik(mts, EGPD.EGPpower, μ = μ, ϕ = ϕ, γ = γ, δ = δ, ξ = ξ, α = α, β = β, κ = κ)
    end

    if markdens == Distributions.GeneralizedPareto
        JuMP.register(model, :to_min_GPD, 7, to_min_GPD, autodiff=true)
    else
        JuMP.register(model, :to_min_EGPD, 8, to_min_EGPD, autodiff=true)
    end

    @variables(model, begin
        mu >= 0, (start = θ[:μ])
        phi >= 0, (start = θ[:ϕ])
        gamma >= 0, (start = θ[:γ])
        delta >= 0, (start = θ[:δ])
        alpha >= 0, (start = θ[:α])
        beta >= 0, (start = θ[:β])
    end)

    if markdens == EGPD.EGPpower
        @variable(model, kappa >= 0, start = θ[:κ])
    end

    if isnothing(bounds)
        @variable(model, xi >= 0, start = θ[:ξ])
    else
        @variable(model, bounds[1] <= xi <= bounds[2], start = θ[:ξ])
    end

    if markdens == Distributions.GeneralizedPareto
        @NLobjective(model, Min, to_min_GPD(mu, phi, gamma, delta, xi, alpha, beta))
    else
        @NLobjective(model, Min, to_min_EGPD(mu, phi, gamma, delta, xi, alpha, beta, kappa))
    end

    optimize!(model)

    sempp.μ = value(mu)
    sempp.ϕ = value(phi)
    sempp.γ = value(gamma)
    sempp.δ = value(delta)
    sempp.ξ = value(xi)
    sempp.α = value(alpha)
    sempp.β = value(beta)

    x = [sempp.μ, sempp.ϕ, sempp.γ, sempp.δ, sempp.ξ, sempp.α, sempp.β]

    if markdens == Distributions.GeneralizedPareto
        to_hess1(θ) = to_min_GPD(θ...)
        hess = ForwardDiff.hessian(to_hess1, x)
    else
        sempp.κ = value(kappa)
        push!(x, sempp.κ)
        to_hess2(θ) = to_min_EGPD(θ...)
        hess = ForwardDiff.hessian(to_hess2, x)
    end

    sepp.cov_mat = inv(hess)

    return objective_value(model)
end

=#


text:
using BenchmarkTools
using Lava: @scope, @autodestroy, @class

mutable struct A
    a
    function A(a)
        this = new(a)
        # println("Create ", this)
        return this
    end
end
@class A

function destroy(a::A)
    # println("Destroy ", a)
    a.a = 0
end

function testAutodestroy()
    @scope begin
        a = @autodestroy(A(1))
        a = @autodestroy(A(2))
        a = @autodestroy(A(3))
        b = @autodestroy(A(10))
        @scope begin
            local a = @autodestroy(A(4))
            local a = @autodestroy(A(5))
            local a = @autodestroy(A(6))
        end
    end
end

function testManualDestroy()
    a1 = A(1)
    a2 = A(2)
    a3 = A(3)
    b1 = A(10)
    a4 = A(4)
    a5 = A(5)
    a6 = A(6)
    destroy(a6)
    destroy(a5)
    destroy(a4)
    destroy(b1)
    destroy(a3)
    destroy(a2)
    destroy(a1)
end

function test()
    println("----Call auto destroy:")
    b = @benchmarkable testAutodestroy()
    tune!(b)
    display(run(b))

    println("\n----Call manual destroy:")
    b = @benchmarkable testManualDestroy()
    tune!(b)
    display(run(b))
    println()
end

test()


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test

source_terms = source_terms_convergence_test

# BCs must be passed as Dict
boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict(
  :all => boundary_condition
)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# Deformed rectangle that looks like a waving flag,
# lower and upper faces are sinus curves, left and right are vertical lines.
f1(s) = SVector(-1.0, s - 1.0)
f2(s) = SVector( 1.0, s + 1.0)
f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
faces = (f1, f2, f3, f4)

Trixi.validate_faces(faces)
mapping_flag = Trixi.transfinite_mapping(faces)

# Get the uncurved mesh from a file (downloads the file if not available locally)
# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                              mesh_file)

mesh = P4estMesh{2}(mesh_file, polydeg=3,
                    mapping=mapping_flag,
                    initial_refinement_level=1)

# Refine bottom left quadrant of each tree to level 2
function refine_fn(p4est, which_tree, quadrant)
  if quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 2
    # return true (refine)
    return Cint(1)
  else
    # return false (don't refine)
    return Cint(0)
  end
end

# Refine recursively until each bottom left quadrant of a tree has level 2
# The mesh will be rebalanced before the simulation starts
refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback)
###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using BinaryBuilder

name = "GnuTLS"
version = v"3.7.1"

# Collection of sources required to build GnuTLS
sources = [
    ArchiveSource("https://www.gnupg.org/ftp/gcrypt/gnutls/v$(version.major).$(version.minor)/gnutls-$(version).tar.xz",
                  "3777d7963eca5e06eb315686163b7b3f5045e2baac5e54e038ace9835e5cac6f"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gnutls-*/

# Grumble-grumble apple grumble-grumble broken linkers...
#if [[ ${target} == *-apple-* ]]; then
#    export AR=/opt/${target}/bin/ar
#fi

if [[ ${target} == *darwin* ]]; then
    # Fix undefined reference to "_c_isdigit"
    # See https://gitlab.com/gnutls/gnutls/-/issues/1033
    atomic_patch -p1 ../patches/03-undo-libtasn1-cisdigit.patch

    # We need to explicitly request a higher `-mmacosx-version-min` here, so that it doesn't
    # complain about: `Symbol not found: ___isOSVersionAtLeast`
    if [[ "${target}" == aarch64* ]]; then
        export CFLAGS="-mmacosx-version-min=11.0"
    else
        export CFLAGS="-mmacosx-version-min=10.11"
    fi
fi

GMP_CFLAGS="-I${prefix}/include" ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-included-libtasn1 \
    --with-included-unistring \
    --without-p11-kit

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# Disable windows because O_NONBLOCK isn't defined
filter!(!Sys.iswindows, platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("GMP_jll", v"6.2.0"),
    Dependency("Nettle_jll", v"3.7.2"; compat="~3.7.2"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
	xyz2grd(cmd0::String="", arg1=nothing; kwargs...)

Convert data table to a grid file.

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

Parameters
----------

- $(GMT.opt_I)
    ($(GMTdoc)xyz2grd.html#i)
- $(GMT.opt_R)
- **A** | **multiple_nodes** :: [Type => Str]      `Arg = [d|f|l|m|n|r|S|s|u|z]`

    By default we will calculate mean values if multiple entries fall on the same node.
    Use A to change this behavior.
    ($(GMTdoc)xyz2grd.html#a)
- **D** | **header** :: [Type => Str]  `Arg = [+xxname][+yyname][+zzname][+sscale][+ooffset][+ninvalid][+ttitle][+rremark]`

    Output edges
    ($(GMTdoc)xyz2grd.html#d)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdclip(....) form.
    ($(GMTdoc)xyz2grd.html#g)
- $(GMT.opt_J)
- **S** | **swap** :: [Type => Str | []]        `Arg = [zfile]`

    Swap the byte-order of the input only. No grid file is produced.
    ($(GMTdoc)xyz2grd.html#s)
- $(GMT.opt_V)
- **Z** | **flags** :: [Type => Str]

    Read a 1-column table. This assumes that all the nodes are present and sorted according to specified ordering convention contained in. ``flags``.
    ($(GMTdoc)xyz2grd.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)
"""
function xyz2grd(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("xyz2grd", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:G :RIr :J :V_params :bi :di :e :f :h :i :w :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :multiple_nodes], [:D :header], [:S :swap], [:Z :flags]])
	if (cmd0 == "" && arg1 === nothing)
		(haskey(d, :x) && haskey(d, :y) && haskey(d, :z)) && (arg1 = hcat(d[:x], d[:y], d[:z]))
        (arg1 !== nothing) && (delete!(d, :x); delete!(d, :y); delete!(d, :z))
	end
	common_grd(d, cmd0, cmd, "xyz2grd ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
xyz2grd(arg1, cmd0::String=""; kw...) = xyz2grd(cmd0, arg1; kw...)


text:
include("../ZMQ_utils.jl")

function n2b2n(ar)
    bin_data, bdds = data2bin!!(; nums = ar, kwarg_name = "nn")
    @show bdds[1]
    # bin2nums(;bin_data=bin_data, bindata_descr=bdds)
    bin2nums(; bin_data = UInt8.(collect(1:24)), bindata_descr = bdds)

end

function buildbin(sz = [24], tp = UInt8; kwn = "nn")
    lng = prod(sz)
    ar1d = tp.(1:lng)
    if tp in (ComplexF32, ComplexF64)
        ar1d = cmplxswap.(ar1d)
    end
    bin_data = reinterpret(UInt8, ar1d)
    if tp == ComplexF32
        tpstr = "ComplexF32"
    elseif tp == ComplexF64
        tpstr = "ComplexF64"
    else
        tpstr = string(tp)
    end

    bdd = BinDescr(1, length(bin_data), sz, tpstr, kwn, "numarrays")
    return (; bin_data, bdd)
end

function b2n2b(sz = [24], tp = UInt8; kwn = "nn")
    bdin, bddin = buildbin(sz, tp; kwn = kwn)
    @show bddin
    nums = bin2nums(; bin_data = bdin, bindata_descr = [bddin])[Symbol(kwn)]
    bdout, bdds = data2bin!!(; arrdata = nums, kwarg_name = kwn)
    ok = bdin == bdout
    @show ok
    return (; ok, bdin, bdout, bddin, bddout = bdds[1], nums, numsint = Int.(abs.(nums)))
end


text:
module LibTreesitter
using CEnum
using tree_sitter_c_jll

include(joinpath(@__DIR__, "..", "gen", "libtreesitter_h.jl"))
include(joinpath(@__DIR__, "..", "gen", "libtreesitter_c.jl"))

# export everything
foreach(names(@__MODULE__, all = true)) do s
    if startswith(string(s), r"ts"i)
        @eval export $s
    end
end

end # module LibTreesitter


text:
# Autogenerated wrapper script for EarCut_jll for x86_64-apple-darwin
export libearcut

JLLWrappers.@generate_wrapper_header("EarCut")
JLLWrappers.@declare_library_product(libearcut, "@rpath/libearcut.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libearcut,
        "lib/libearcut.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    solve(sso::SteadyStateObjective, [method, options, x0]) :: SteadyStateFitResult

Optimize a [`SteadyStateObjective`](@ref) `sso` using Optim.  The rest
of the positional and keyword arguments are passed to
`Optim.optimize`.
"""
function solve(
    sso::SteadyStateObjective,
    method::Optim.AbstractOptimizer = Optim.BFGS(),
    options::Optim.Options = Optim.Options(),
    x0 = get(sso.p, sso.parameterlens);
    kwargs...,
)
    result = Optim.optimize(
        Optim.only_fg!(sso),
        x0,
        method,
        setupoptions(sso, options);
        kwargs...,
    )
    return SteadyStateFitResult(result, sso)
end

"""
    SteadyStateFitResult

# Fields
- `result :: Optim.OptimizationResults`
- `objective :: SteadyStateObjective`

# Methods
- `minimum(fit::SteadyStateFitResult)`: Equivalent to `minimum(fit.result)`.
- `Optim.minimizer(fit::SteadyStateFitResult)`: Return a copy of
  `fit.objective.p` with parameters replaced by `Optim.minimizer(fit.result)`.
"""
struct SteadyStateFitResult
    result  # Optim.MultivariateOptimizationResults
    objective  # SteadyStateObjective
end

function Base.show(io::IO, fit::SteadyStateFitResult)
    maybe_default_show(io, fit) && return
    print(io, "SteadyStateFitResult: ")
    if Optim.converged(fit)
        printstyled(io, "success"; color=:green)
    else
        printstyled(io, "failure"; color=:red)
    end
    print(io, "; ")
    _summarybody(io, fit.objective)
end

function Base.show(io::IO, ::MIME"text/plain", fit::SteadyStateFitResult)
    println(io, "SteadyStateFitResult")
    _showbody(io, fit.objective)
    println(io)
    show(io, MIME("text/plain"), fit.result)
end

function Optim.minimizer(fit::SteadyStateFitResult)
    p = fit.objective.p
    parameterlens = fit.objective.parameterlens
    return set(p, parameterlens, Optim.minimizer(fit.result))
end

# Forward accessor methods to `Optim`.  See:
# https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/#complete-list-of-functions
for f in [
    :(Base.minimum),
    :(Base.summary),
    # Generic optimization:
    :(Optim.iterations),
    :(Optim.iteration_limit_reached),
    :(Optim.trace),
    :(Optim.x_trace),
    :(Optim.f_trace),
    :(Optim.f_calls),
    :(Optim.converged),
    # Multivariate optimization:
    :(Optim.g_norm_trace),
    :(Optim.g_calls),
    :(Optim.x_converged),
    :(Optim.f_converged),
    :(Optim.g_converged),
    :(Optim.initial_state),
]
    @eval $f(fit::SteadyStateFitResult) = $f(fit.result)
end

"""
    steadystates(fit::SteadyStateFitResult)

Compute steady states with the minimizer.
"""
steadystates(fit::SteadyStateFitResult) =
    steadystates(fit.objective, Optim.minimizer(fit.result))

# A hacky solution: To use the steady state of the last iteration,
# `updatesteadystates!` is called via `callback` option.

function setupoptions(sso::SteadyStateObjective, options::Optim.Options)
    @assert options.show_every == 1
    options = @set options.extended_trace = true
    options = @set options.callback =
        StateUpdater(sso, options.callback, options.store_trace)
    return options
end

struct StateUpdater{T, C}
    sso::T
    callback::C
    store_trace::Bool
end

function (su::StateUpdater)(os)
    if su.store_trace
        # os :: Vector{<:OptimizationState}
        updatesteadystates!(su.sso, os[end].metadata["x"])
    else
        # os :: OptimizationState
        updatesteadystates!(su.sso, os.metadata["x"])
    end
    su.callback === nothing ? false : su.callback(os)
end


text:
# prints solution of both parts
function x()
	f = open("4.txt")
	c1 = 0
	c2 = 0
	s1 = 0
	s2 = 0
	for i=1:954
		l = readline(f)
		if l == ""
			if s1 >= 7
				c1 += 1
			end
			if s2 >= 7
				c2 += 1
			end
			s1 = 0
			s2 = 0
			continue
		end
		for kv in split(l, " ")
			k = split(kv, ":")[1]
			v = split(kv, ":")[2]
			if k == "byr"
				s1 += 1
				y = parse(Int, v)
				if 1920 <= y <= 2002
					s2 += 1
				end
			elseif k == "iyr"
				s1 += 1
				y = parse(Int, v)
				if 2010 <= y <= 2020
					s2 += 1
				end
			elseif k == "eyr"
				s1 += 1
				y = parse(Int, v)
				if 2020 <= y <= 2030
					s2 += 1
				end
			elseif k == "hgt"
				s1 += 1
				a = findfirst(isequal('c'), v)
				b = findfirst(isequal('i'), v)
				if typeof(a) != Nothing
					cm = parse(Int, v[1:a-1])
					if 150 <= cm <= 193
						s2 += 1
					end
				elseif typeof(b) != Nothing
					inch = parse(Int, v[1:b-1])
					if 59 <= inch <= 76
						s2 += 1
					end
				end
			elseif k == "hcl"
				s1 += 1
				if v[1] == '#' && length(v) == 7
					try
						h = parse(Int, v[2:7], base=16)
						s2 += 1
					catch
						s2 += 0
					end
				end
			elseif k == "ecl"
				s1 += 1
				if v in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
					s2 += 1
				end
			elseif k == "pid"
				s1 += 1
				if length(v) == 9
					try
						h = parse(Int, v)
						s2 += 1
					catch
						s2 += 0
					end
				end
			end
		end
	end
	println(c1)
	println(c2)
end


text:
struct Shape1D <: FluxLayer end
NaiveNASflux.indim(::Shape1D) = 1
NaiveNASflux.outdim(::Shape1D) = 1
NaiveNASflux.actdim(::Shape1D) = 1
NaiveNASflux.actrank(::Shape1D) = 0

# TODO: Move to NaiveNASflux
NaiveNASflux.nin(sc::SkipConnection) = nin(sc.layers)

numpy2fluxdim(np_axis, v::AbstractVertex) = numpy2fluxdim(np_axis, 1 + NaiveNASflux.actrank(v)[1])
numpy2fluxdim(np_axis, ndims) = np_axis >= 0 ? ndims - np_axis : abs(np_axis)

flux2numpydim(dim, ndims) = ndims - dim

function shape(v::AbstractVertex)
    outshape = shape(layertype(v), nout(v))
    ismissing(outshape) && return first(unique(shape.(inputs(v))))
    return outshape
end
shape(::FluxLayer, outsize) = missing
shape(::Shape1D, outsize) = (outsize,)
shape(::Flux2D, outsize) = (outsize, missing)
shape(::FluxConvolutional{N}, outsize) where N = ((missing for _ in 1:N)..., outsize, missing)
shape(::FluxRecurrent, outsize) = (outsize, missing, missing)

aggshape(f, d::Number...) = f(d...)
aggshape(f, d...) = missing

rmdims(t::Tuple, dim::Integer) = t[1:end .!= dim]
rmdims(t::Tuple, dims) = Tuple(t[i] for i in 1:length(t) if i ∉ dims)

function guess_layertype(ndims::Integer)
    ndims <= 1 && return Shape1D()
    ndims == 2 && return GenericFlux2D()
    ndims == 3 && return GenericFluxRecurrent()
    return GenericFluxConvolutional{ndims-2}()
end

flipweights(l, w) = w
flipweights(::FluxConvolutional{N}, w) where N = w[(size(w,i):-1:1 for i in 1:N)..., :, :]
flipweights(::FluxRnn, w, hsize) = w
function flipweights(::FluxLstm, w, hsize)
    input = Flux.gate(w, hsize, 1)
    forget = Flux.gate(w, hsize, 2)
    cell = Flux.gate(w, hsize, 3)
    output = Flux.gate(w, hsize, 4)
    return vcat(input, output, forget, cell)
end

unflipweights(::FluxRnn, w, hsize) = w
function unflipweights(::FluxLstm, w, hsize)
    input = Flux.gate(w, hsize, 1)
    output = Flux.gate(w, hsize, 2)
    forget = Flux.gate(w, hsize, 3)
    cell = Flux.gate(w, hsize, 4)
    return vcat(input, forget, cell, output)
end

outshape(l, s) = outshape(layertype(l), l, s)
outshape(::FluxParLayer, l, ::Missing) = outshape(l, shape(layertype(l), nin(l)[]))
outshape(lt, l, ::Missing) = missing

function outshape(::Flux2D, l, s::Tuple)
    assertshape(s, 2, l)
    assertsize(s[1], nin(l)[], l)
    return (nout(l), s[end])
end
function outshape(::FluxRecurrent, l, s::Tuple)
    # l might be a cell and nin/nout in NaiveNASflux is only defined for Flux.Recur{cell}.
    assertshape(s, (3, 4), l)
    assertsize(s[1], nin(l)[], l)
    # ONNX wants num directions as an extra dimension to output
    # Fluxs RNNs are not bidirectional so num directions is always 1
    return (nout(l), s[2], 1, s[end])
end

function outshape(::FluxConvolutional{N}, l, s::Tuple) where N
    assertshape(s, N+2, l)
    assertsize(s[N+1], nin(l)[], l)
    p = length(l.pad) == N ? 2 .* l.pad : l.pad[1:2:end] .+ l.pad[2:2:end]
    k = size(weights(l))[1:N]
    d = l.dilation
    stride = l.stride

    o = map(zip(1:N, s)) do (i, si)
        # Conv arithmetic from https://arxiv.org/pdf/1603.07285.pdf
        aggshape(x -> (x + p[i] - k[i] - (k[i] - 1)*(d[i] - 1)) ÷ stride[i] + 1, si)
    end
    return (o..., nout(l), s[N+2])
end

outshape(l::Union{Flux.MaxPool{N}, Flux.MeanPool{N}}, ::Missing) where N = outshape(l, ntuple(i->missing, N+2))
function outshape(l::Union{Flux.MaxPool{N}, Flux.MeanPool{N}}, s::Tuple) where N
    assertshape(s, N+2, l)
    p = length(l.pad) == N ? 2 .* l.pad : l.pad[1:2:end] .+ l.pad[2:2:end]
    k = l.k
    stride = l.stride

    o = map(zip(1:N, s)) do (i, si)
        # Conv arithmetic from https://arxiv.org/pdf/1603.07285.pdf
        aggshape(x -> (x + p[i] - k[i]) ÷ stride[i] + 1, si)
    end

    return (o..., s[N+1], s[N+2])
end

function assertshape(s, expected, lmsg)
    if all(!=(length(s)), expected)
        throw(DimensionMismatch("Wrong input dimension for $(lmsg)! Expected $(join(expected,", ", " or ")) dimensions, got shape $s"))
    end
end
function assertsize(s::Integer, expected, lmsg)
    if s != expected
        throw(DimensionMismatch("Wrong input size for $(lmsg)! Expected $expected elements, got size $s"))
    end
end
function assertsize(s, expected, lmsg) end


text:
immutable Frenet
    roadind::RoadIndex
    s::Float64 # distance along lane
    t::Float64 # lane offset, positive is to left
    ϕ::Float64 # lane relative heading
end
function Frenet(roadind::RoadIndex, roadway::Roadway; t::Float64=0.0, ϕ::Float64=0.0)
    s = roadway[roadind].s
    ϕ = _mod2pi2(ϕ)
    Frenet(roadind, s, t, ϕ)
end
function Frenet(roadproj::RoadProjection, roadway::Roadway)
    roadind = RoadIndex(roadproj.curveproj.ind, roadproj.tag)
    s = roadway[roadind].s
    t = roadproj.curveproj.t
    ϕ = _mod2pi2(roadproj.curveproj.ϕ)
    Frenet(roadind, s, t, ϕ)
end
Frenet(posG::VecSE2, roadway::Roadway) = Frenet(proj(posG, roadway), roadway)
Frenet(posG::VecSE2, lane::Lane, roadway::Roadway) = Frenet(proj(posG, lane, roadway), roadway)

function get_posG(frenet::Frenet, roadway::Roadway)
    curvept = roadway[frenet.roadind]
    pos = curvept.pos + polar(frenet.t, curvept.pos.θ + π/2)
    VecSE2(pos.x, pos.y, frenet.ϕ + curvept.pos.θ)
end

const NULL_FRENET = Frenet(NULL_ROADINDEX, NaN, NaN, NaN)

Base.show(io::IO, frenet::Frenet) = print(io, "Frenet(", frenet.roadind, @sprintf(", %.3f, %.3f, %.3f)", frenet.s, frenet.t, frenet.ϕ))
function Base.isapprox(a::Frenet, b::Frenet;
    rtol::Real=cbrt(eps(Float64)),
    atol::Real=sqrt(eps(Float64))
    )

    a.roadind.tag == b.roadind.tag &&
    isapprox(a.roadind.ind.t, b.roadind.ind.t, atol=atol, rtol=rtol) &&
    isapprox(a.s, b.s, atol=atol, rtol=rtol) &&
    isapprox(a.t, b.t, atol=atol, rtol=rtol) &&
    isapprox(a.ϕ, b.ϕ, atol=atol, rtol=rtol)
end


text:
Base.@pure addtype(T) = typeof(zero(T)+log1p(exp(zero(T))))

# generic

exp(x::AnyLog) = uexp(x)

# ULogarithmic

log(x::ULogarithmic) = x.log
# *(x::ULogarithmic{T}, y::ULogarithmic{T}, rest::ULogarithmic{T}...) where {T} = uexp(T, +(log.((x, y, rest...))...))
*(x::ULogarithmic{T}, y::ULogarithmic{T}) where {T} = uexp(x.log + y.log)
/(x::ULogarithmic{T}, y::ULogarithmic{T}) where {T} = uexp(x.log - y.log)
^(x::ULogarithmic, n::Real) = uexp(x.log * n)
^(x::ULogarithmic, n::Integer) = uexp(x.log * n)
inv(x::ULogarithmic{T}) where {T} = uexp(-x.log)
+(x::ULogarithmic{T}) where {T} = uexp(addtype(T), x.log)
+(x::ULogarithmic{T}, y::ULogarithmic{T}) where {T} =
	let xlog=x.log, ylog=y.log
		uexp(addtype(T), if xlog == ylog
			xlog + log1p(one(T))
		elseif xlog ≥ ylog
			xlog + log1p(exp(ylog - xlog))
		else
			ylog + log1p(exp(xlog - ylog))
		end)
	end
-(x::ULogarithmic{T}, y::ULogarithmic{T}) where {T} =
	let xlog=x.log, ylog=y.log
		uexp(addtype(T), if xlog ≥ ylog
			if iszero(x)
				xlog
			else
				xlog + log1p(-exp(ylog - xlog))
			end
		else
			throw(DomainError((x, y), "difference is negative"))
		end)
	end

# Logarithmic

-(x::Logarithmic) = Logarithmic(x.abs, !x.signbit)
log(x::Logarithmic) = (x.signbit && !iszero(x) && throw(DomainError(x)); log(x.abs))
*(x::Logarithmic{T}, y::Logarithmic{T}) where {T} = Logarithmic(x.abs * y.abs, x.signbit ⊻ y.signbit)
/(x::Logarithmic{T}, y::Logarithmic{T}) where {T} = Logarithmic(x.abs / y.abs, x.signbit ⊻ y.signbit)
^(x::Logarithmic, n::Integer) = Logarithmic(x.abs^n, x.signbit & isodd(n))
^(x::Logarithmic, n::Real) = (x.signbit && !iszero(x) && throw(DomainError(x)); Logarithmic(x.abs^n))
inv(x::Logarithmic{T}) where {T} = Logarithmic(inv(x.abs), x.signbit)
+(x::Logarithmic) = Logarithmic(+x.abs, x.signbit)
+(x::Logarithmic{T}, y::Logarithmic{T}) where {T} = x.signbit==y.signbit ? Logarithmic(x.abs + y.abs, x.signbit) : x.abs ≥ y.abs ? Logarithmic(x.abs - y.abs, x.signbit) : Logarithmic(y.abs - x.abs, y.signbit)
-(x::Logarithmic{T}, y::Logarithmic{T}) where {T} = x.signbit==y.signbit ? x.abs ≥ y.abs ? Logarithmic(x.abs - y.abs, x.signbit) : Logarithmic(y.abs - x.abs, !y.signbit) : Logarithmic(x.abs + y.abs, x.signbit)

# CLogarithmic

log(x::CLogarithmic) = Complex(x.abs.log, x.angle)
conj(x::CLogarithmic) = CLogarithmic(x.abs, -x.angle)
*(x::CLogarithmic{T}, y::CLogarithmic{T}) where {T} =
	CLogarithmic(x.abs * y.abs, x.angle + y.angle)
/(x::CLogarithmic{T}, y::CLogarithmic{T}) where {T} =
	CLogarithmic(x.abs / y.abs, x.angle - y.angle)
^(x::CLogarithmic, n::Integer) = CLogarithmic(x.abs^n, x.angle*n)
^(x::CLogarithmic, n::Real) = CLogarithmic(x.abs^n, x.angle*n)
inv(x::CLogarithmic) = CLogarithmic(inv(x.abs), -x.angle)
+(x::CLogarithmic) = CLogarithmic(+x.abs, x.angle)
+(x::CLogarithmic{T}, y::CLogarithmic{T}) where {T} =
	if isinf(x) && isinf(y) && x.angle ≈ y.angle
		t = exp(zero(Complex{T}))
		x * t
	elseif iszero(x) && iszero(y)
		t = exp(zero(Complex{T}))
		x * t
	elseif x.abs ≥ y.abs
		t = exp(Complex(y.abs.log-x.abs.log,y.angle-x.angle))
		t += one(t)
		x * t
	else
		t = exp(Complex(x.abs.log-y.abs.log,x.angle-y.angle))
		t += one(t)
		t * y
	end
-(x::CLogarithmic{T}, y::CLogarithmic{T}) where {T} =
	if iszero(x) && iszero(y)
		t = exp(zero(Complex{T}))
		x * t
	elseif x.abs ≥ y.abs
		t = exp(Complex(y.abs.log-x.abs.log,y.angle-x.angle))
		t = one(t) - t
		x * t
	else
		t = exp(Complex(x.abs.log-y.abs.log,x.angle-y.angle))
		t = t - one(t)
		t * y
	end


text:
using FillArrays
using StaticArrays
using Distances
using OffsetArrays
using DataStructures
using Statistics

const sqrt2 = sqrt(2.0)
const distance_8 = @SMatrix[sqrt2 1 sqrt2; 1 Inf 1; sqrt2 1 sqrt2]
const distance_4 = @SMatrix[Inf 1 Inf; 1 Inf 1; Inf 1 Inf]
const Δ = CartesianIndex(1, 1)

"""
```
spread(points::Matrix{<:Real}, initial::Matrix{<:Real}, friction::Matrix{<:Real}; res=1, limit=Inf)
```
Total friction distance spread from `points` as by *Tomlin (1983)* [^tomlin1983].
This is also the method implemented by [PCRaster](https://pcraster.geo.uu.nl/pcraster/4.0.2/doc/manual/op_spread.html).

# Output
- `Array{Float64,2}` Total friction distance

# Arguments
- `points::Matrix{<:Real}` Input Array
- `initial::Matrix{<:Real}` Factor to exaggerate elevation
- `friction::Matrix{<:Real}` Resolution of cell size
- `res=1` Resolution or cell size
- `limit=Inf` Initial fill value

[^tomlin1983]: Tomlin, Charles Dana. 1983. Digital Cartographic Modeling Techniques in Environmental Planning. Yale University.
"""
function spread(points::AbstractMatrix{<:Real}, initial::AbstractMatrix{<:Real}, friction::AbstractMatrix{<:Real}; res=1, limit=Inf)

    ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction

    result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    r = @view result[1:end - 1,1:end - 1]
    locations = points .> 0
    r[locations] .= initial[locations]

    # Construct stack for locations
    mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
    mask[begin + 1:end - 1,begin + 1:end - 1] .= false

    II = CartesianIndices(size(points))
    stack = Deque{CartesianIndex}()
    for I in II[locations]
        push!(stack, I)
    end

    # Step 1: Set the distance of the starting node to 0 and the distances of all other nodes to the highest value possible.
    sdata = zeros(MMatrix{3,3})
    mcell = MMatrix{3,3}(false, false, false, false, false, false, false, false, false)

    # Step 3: For each of the active node’s adjacent neighbors, set its distance to whichever is
    # less: its current distance value or the sum of the distance of the active node plus the
    # weight of the arc from the active node to that neighbor.
    while !isempty(stack)
        spread!(stack, mask, result, ofriction, sdata, mcell, res)
    end
    r
end

function spread!(stack, mask, result, ofriction, sdata, mcell, res)
    I = popfirst!(stack)
    mask[I] = true
    patch = I - Δ:I + Δ

    rdata = view(result, patch)
    fdata = view(ofriction, patch)

    # New distance is cell_distance + average friction values
    for i ∈ eachindex(sdata)
        sdata[i] = muladd(fdata[i] + fdata[2,2], res / 2 * distance_8[i], rdata[2,2])
        mcell[i] = sdata[i] < rdata[i]  # cells where new distance is lower
    end
    rdata[mcell] .= sdata[mcell]
    result[patch] .= rdata

    # Add new cells to stack
    for I in patch[mcell]
        mask[I] || push!(stack, I)
    end
end

"""
```
spread2(points::Matrix{<:Real}, initial::Matrix{<:Real}, friction::Matrix{<:Real}; res=1, limit=Inf, iterations=3)
```

Pushbroom method for friction costs as discussed by *Eastman (1989) [^eastman1989].
This method should scale much better (linearly) than the [^tomlin1983] method, but can require more
`iterations` than set by default (3) in the case of maze-like, uncrossable obstacles.

# Output
- `Array{Float64,2}` Total friction distance

# Arguments
- `points::Matrix{<:Real}` Input Array
- `initial::Matrix{<:Real}` Factor to exaggerate elevation
- `friction::Matrix{<:Real}` Resolution of cell size
- `res=1` Resolution or cell size
- `limit=Inf` Initial fill value
- `iterations=3` Number of pushbroom iterations

[^eastman1989]: Eastman, J. Ronald. 1989. ‘Pushbroom Algorithms for Calculating Distances in Raster Grids’. In Proceedings, Autocarto, 9:288–97.
"""
function spread2(points::AbstractMatrix{<:Real}, initial::AbstractMatrix{<:Real}, friction::AbstractMatrix{<:Real}; res=1, limit=Inf, iterations=3)

    ofriction = OffsetMatrix(fill(Inf, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    ofriction[begin + 1:end - 1,begin + 1:end - 1] .= friction

    result = OffsetMatrix(fill(limit, size(friction) .+ 2), UnitRange.(0, size(points) .+ 1))
    r = @view result[1:end - 1,1:end - 1]
    locations = points .> 0
    r[locations] .= initial[locations]

    mask = OffsetMatrix(trues(size(points) .+ 2), UnitRange.(0, size(points) .+ 1))
    mask[begin + 1:end - 1,begin + 1:end - 1] .= false

    minval, minidx = [0.], [CartesianIndex(1, 1)]
    x = @MMatrix zeros(3, 3)

    indices = CartesianIndices(size(points))
    for i in 1:iterations
        II = (i % 2 == 1) ? indices : reverse(indices)
        for I ∈ II
            patch = I - Δ:I + Δ

            rdata = view(result, patch)
            fdata = view(ofriction, patch)

            x .= (fdata .+ fdata[2,2]) .* res ./ 2 .* distance_8 .+ rdata
            findmin!(minval, minidx, x)
            rdata[2,2] = min(rdata[2,2], minval[1])
        end
    end
    r
end

"""
```
spread(points::Matrix{<:Real}, initial::Real, friction::Real; distance=Euclidean(), res=1.0)
```
"""
function spread(points::AbstractMatrix{<:Real}, initial::Real, friction::Real; distance=Euclidean(), res=1.0)
    init = fill(initial, size(points))
    return spread(points, init, friction; distance, res)
end

"""
```
spread(points::Matrix{<:Real}, initial::Matrix{<:Real}, friction::Real; distance=Euclidean(), res=1.0)
```
Optimized (and more accurate) function based on the same friction everywhere.

When the friction is the same everywhere, there's no need for searching the shortest cost path,
as one can just take a direct line to the input points.

The calculated cost is more accurate, as there's no 'zigzag' from cell center to cell center.
"""
function spread(points::AbstractMatrix{<:Real}, initial::AbstractMatrix{<:Real}, friction::Real; distance=Euclidean(), res=1.0)
    locations = points .> 0
    I = CartesianIndices(size(points))

    result = fill(Inf, size(points))
    for location ∈ I[locations]
        for cell ∈ I
            result[cell] = min(evaluate(distance, location.I, cell.I) * res * friction + initial[location], result[cell])
        end
    end
    m = .~isfinite.(points)
    result[m] .= points[m]
    return result
end


text:
"""
    CoamoebaFiber2D(f::AbstractPolynomial, θ=(0., 0.))

Construct the fiber `CoamoebaFiber2D` of the bivariate polynomial `f` at `θ`.
"""
function CoamoebaFiber2D(p::MP.AbstractPolynomialLike{S}, θ=(0.0, 0.0)) where S
    vars = MP.variables(p)
    @assert (length(vars) == 2) "Expected a bivariate polynomial, but got $(p)."

    T = realified_coefficient_type(S)
    f_re, f_im = SP.Polynomial.(one(T) .* realify(p))

    sincosθ = _sincosθ(θ)
    U = zeros(2, 4)
    x = SVector(0.0, 0.0, 0.0, 0.0)
    CoamoebaFiber2D(f_re, f_im, sincosθ, U, x)
end

_sincosθ(θ) = (sin(θ[1]), sin(θ[2]), cos(θ[1]), cos(θ[2]))

function update_fiber!(F::CoamoebaFiber2D, θ)
    F.sincosθ = _sincosθ(θ)
    F
end

function phi!(fiber::CoamoebaFiber2D, w)
    s1, s2, c1, c2 = fiber.sincosθ
    v1, v2 = exp(w[1]), exp(w[2])
    x = SVector(v1 * c1, v2 * c2, v1 * s1, v2 * s2)
    fiber.x = x
    x
end

function evaluate(fiber::CoamoebaFiber2D, w)
    x = phi!(fiber, w)
    SVector(SP.evaluate(fiber.f_re, x), SP.evaluate(fiber.f_im, x))
end

function jacobian(fiber::CoamoebaFiber2D, w, precomputed=false)
    x = precomputed ? fiber.x : phi!(fiber, w)
    U = fiber.U
    @inbounds begin
        ∇f_re = SP.gradient(fiber.f_re, x)
        for i=1:4
            U[1, i] = ∇f_re[i]
        end
        U[2, 3] = U[1, 1]
        U[2, 4] = U[1, 2]
        U[2, 1] = -U[1, 3]
        U[2, 2] = -U[1, 4]
    end

    @inbounds out = begin
        ∂re∂w1 = muladd(U[1, 1], x[1], U[1, 3] * x[3])
        ∂re∂w2 = muladd(U[1, 2], x[2], U[1, 4] * x[4])
        ∂im∂w1 = muladd(U[2, 1], x[1], U[2, 3] * x[3])
        ∂im∂w2 = muladd(U[2, 2], x[2], U[2, 4] * x[4])

        SMatrix{2,2}(∂re∂w1, ∂im∂w1, ∂re∂w2, ∂im∂w2)
    end
    out
end


text:
using AbstractMCMC
using AdvancedMH
using Clustering
using Distributions
using Distributions.FillArrays
using DistributionsAD
using FiniteDifferences
using ForwardDiff
using GalacticOptim
using MCMCChains
using Memoization
using NamedArrays
using Optim
using PDMats
using ReverseDiff
using SpecialFunctions
using StatsBase
using StatsFuns
using Tracker
using Turing
using Turing.Inference
using Turing.RandomMeasures
using Zygote

using LinearAlgebra
using Pkg
using Random
using Test
using StableRNGs

using AdvancedPS: ResampleWithESSThreshold, resample_systematic, resample_multinomial
using AdvancedVI: TruncatedADAGrad, DecayedADAGrad, apply!
using Distributions: Binomial, logpdf
using DynamicPPL: getval, getlogp
using ForwardDiff: Dual
using MCMCChains: Chains
using StatsFuns: binomlogpdf, logistic, logsumexp
using TimerOutputs: TimerOutputs, @timeit
using Turing: BinomialLogit, ForwardDiffAD, Sampler, SampleFromPrior, NUTS, TrackerAD,
                Variational, ZygoteAD, getspace, gradient_logp
using Turing.Essential: TuringDenseMvNormal, TuringDiagMvNormal
using Turing.Variational: TruncatedADAGrad, DecayedADAGrad, AdvancedVI

setprogress!(false)

include(pkgdir(Turing)*"/test/test_utils/AllUtils.jl")

# Collect timing and allocations information to show in a clear way.
const TIMEROUTPUT = TimerOutputs.TimerOutput()
macro timeit_include(path::AbstractString) :(@timeit TIMEROUTPUT $path include($path)) end

@testset "Turing" begin
    @testset "essential" begin
        @timeit_include("essential/ad.jl")
    end

    @testset "samplers (without AD)" begin
        @timeit_include("inference/AdvancedSMC.jl")
        @timeit_include("inference/emcee.jl")
        @timeit_include("inference/ess.jl")
        @timeit_include("inference/is.jl")
    end

    Turing.setrdcache(false)
    for adbackend in (:forwarddiff, :tracker, :reversediff)
        @timeit TIMEROUTPUT "inference: $adbackend" begin
            Turing.setadbackend(adbackend)
            @info "Testing $(adbackend)"
            @testset "inference: $adbackend" begin
                @testset "samplers" begin
                    @timeit_include("inference/gibbs.jl")
                    @timeit_include("inference/gibbs_conditional.jl")
                    @timeit_include("inference/hmc.jl")
                    @timeit_include("inference/Inference.jl")
                    @timeit_include("contrib/inference/dynamichmc.jl")
                    @timeit_include("contrib/inference/sghmc.jl")
                    @timeit_include("inference/mh.jl")
                end
            end

            @testset "variational algorithms : $adbackend" begin
                @timeit_include("variational/advi.jl")
            end

            @testset "modes : $adbackend" begin
                @timeit_include("modes/ModeEstimation.jl")
                @timeit_include("modes/OptimInterface.jl")
            end

        end
    end

    @testset "variational optimisers" begin
        @timeit_include("variational/optimisers.jl")
    end

    Turing.setadbackend(:forwarddiff)
    @testset "stdlib" begin
        @timeit_include("stdlib/distributions.jl")
        @timeit_include("stdlib/RandomMeasures.jl")
    end

    @testset "utilities" begin
        # include("utilities/stan-interface.jl")
        @timeit_include("inference/utilities.jl")
    end
end

show(TIMEROUTPUT; compact=true, sortby=:firstexec)


text:
# I/O functions
export readmusicxml, parsemusicxml, writemusicxml, printmusicxml, pprint, fwritemusicxml, readmusicxml_partial, parsemusicxml_partial

"""
    extractdata(doc)

Extracts musicxml data, builds all the types and stores them in proper format.

This function is not exported. Use readmusicxml and parsemusicxml instead.

# Examples
```julia
data = extractdata(doc)
```
"""
extractdata(doc::Document) = ScorePartwise(doc)
################################################################
"""
    readmusicxml(filepath)

Reads musicxml file and then extracts the data, builds all the types and stores them in proper format.

# Examples
```julia
data = readmusicxml(joinpath("examples", "musescore.musicxml"))
```
"""
function readmusicxml(filepath::String)
    doc = readxml(filepath) # read an XML document from a file
    data = extractdata(doc)
    return data
end
################################################################

"""
Similar to [`readmusicxml`](@ref) but without converting it to a Julia type. Use the appropriate type to convert it

# Examples
```julia
xml_note = readmusicxml_partial(path_to_file)
Note(xml_note)
```
"""
function readmusicxml_partial(filepath::String)
    doc = readxml(filepath) # read an XML document from a file
    return doc
end
################################################################
"""
    parsemusicxml(s)

Parses musicxml from a string and then extracts the data, builds all the types and stores them in proper format.

# Examples
```julia
data = parsemusicxml(s)
```
"""
function parsemusicxml(s::String)
    doc = parsexml(s) # Parse an XML string
    data = extractdata(doc)
    return data
end

################################################################
"""
Similar to [`parsemusicxml`](@ref) but without converting it to a Julia type. Use the appropriate type to convert it

# Examples
```julia
xml_note = parsemusicxml_partial(\"\"\"
<note>
  <grace/>
  <pitch>
    <step>G</step>
    <octave>4</octave>
  </pitch>
  <voice>1</voice>
  <type>16th</type>
  <stem>up</stem>
  <beam number="1">end</beam>
  <beam number="2">end</beam>
</note>
\"\"\")

Note(xml_note)
```
"""
function parsemusicxml_partial(s::String)
    doc = parsexml(s) # Parse an XML string
    return doc
end

################################################################
"""
    writemusicxml(filename::AbstractString, x)

Writes musicxml score into a file.

# Examples
```julia
writemusicxml("myscore.musicxml", score)
```
"""
function writemusicxml end

################################################################
"""
    printmusicxml(x)
    printmusicxml(io, x)

Print musicxml score

# Examples
```julia
printmusicxml(score)
```
"""
function printmusicxml end

writemusicxml(args...) = pprint(args...)
printmusicxml(args...) = pprint(args...)

################################################################
"""
 fwritemusicxml(filename::AbstractString, x)

Fast write
"""
fwritemusicxml(io_filename, score) = Base.write(io_filename, score.aml)


text:
using Celerite2

function logp(term::Celerite2.Term, x, y, sigma)

  N = size(x, 1)
  y0 = copy(y)
  a, U, V, P = Celerite2.get_matrices(term, x, sigma.^2 .+ zeros(N))
  Celerite2.factor!(U, P, a, V)
  Celerite2.solve!(U, P, a, V, y0)
  return -0.5 * (transpose(y) * y0 + sum(log.(a)) + N * log(2*pi))

end