text:
using Test
using LIKWID
using LinearAlgebra

# julia must have been started with `likwid-perfctr -C ... -g ... -m`
@assert Marker.isactive()

# init
@test isnothing(Marker.init())

N = 100_000_000
a = 3.141f0
z = zeros(Float32, N)
x = rand(Float32, N)
y = rand(Float32, N)

saxpy!(z,a,x,y) = z .= a .* x .+ y

@test @region "saxpy!" saxpy!(z,a,x,y)
@test @region "saxpy!" saxpy!(z,a,x,y)

A = rand(100,100)
B = rand(100,100)

@test @region "mul" for _ in 1:10
    A * B
end

# close
@test isnothing(Marker.close())


text:
@testset "Genetic Programming" begin
    rng = StableRNG(42)

    pop = 10
    terms = Terminal[:x, :y, rand]
    funcs = Function[+,-,*,/]

    t = TreeGP(pop, terms, funcs, maxdepth=2)
    @test Evolutionary.population_size(t) == pop
    @test sort(Evolutionary.terminals(t)) == [:x, :y]
    @testset for (term, dim) in t.terminals
        @test dim == 1
    end
    @testset for (func, arity) in t.functions
        @test arity == 2
    end
    show(IOBuffer(), t)

    # population initialization
    popexp = Evolutionary.initial_population(t, rng=rng);
    @test length(popexp) == pop
    popexp = Evolutionary.initial_population(t, :(x + 1));
    @test popexp[1] == :(x + 1)

    # recursive helper functions
    h = 4
    gtr = TreeGP(pop, terms, funcs, maxdepth=2, initialization=:grow)
    Random.seed!(rng, 8237463746)
    tmp = rand(rng, gtr, h)
    Random.seed!(rng, 8237463746)
    gt = rand(rng, gtr, h)
    @test tmp == gt
    @test Evolutionary.nodes(gt) < 2^(h+1)-1
    @test Evolutionary.height(gt) <= h
    @test length(gt) < 2^(h+1)-1
    ft = rand(rng, TreeGP(pop, terms, funcs, maxdepth=2, initialization=:full), h)
    @test Evolutionary.nodes(ft) == 2^(h+1)-1
    @test Evolutionary.height(ft) == h
    @test length(ft) == 2^(h+1)-1
    @test Evolutionary.depth(ft, :x) == 4
    ft[3] = :z
    @test Evolutionary.depth(ft, :z) == 4
    @test Evolutionary.depth(ft, ft) == 0
    @test Evolutionary.depth(ft, ft[3]) > 0
    @test Evolutionary.depth(ft, :w) == -1
    @test Evolutionary.evaluate(:y, Dict(:y=>1, :z=>2), 1.0, 2.0) == 1.0
    copyto!(ft, gt)
    @test ft == gt
    @test Evolutionary.symbols(ft) |> sort == [:x, :y]

    # simplification
    using Evolutionary: simplify!
    @test Expr(:call, -, :x, :x) |> simplify! == 0
    @test Expr(:call, /, :x, :x) |> simplify! == 1
    @test Expr(:call, *, 0, :x) |> simplify! == 0
    @test Expr(:call, *, :x, 0) |> simplify! == 0
    @test Expr(:call, /, 0, :x) |> simplify! == 0
    @test Expr(:call, /, 0,  1) |> simplify! == 0
    @test Expr(:call, +, 0, :x) |> simplify! == :x
    @test Expr(:call, +, :x, 0) |> simplify! == :x
    @test Expr(:call, -, :x, 0) |> simplify! == :x
    @test Expr(:call, +, :x, :x) |> simplify! == Expr(:call, *, 2, :x)
    @test Expr(:call, -, 2, 1) |> simplify! == 1
    @test Expr(:call, exp, Expr(:call, log, 1)) |> simplify! == 1
    @test Expr(:call, log, Expr(:call, exp, 1)) |> simplify! == 1
    @test Expr(:call, -, Expr(:call, +, :x, 1), 2) |> simplify! == Expr(:call, +, :x, -1)
    @test Expr(:call, -, Expr(:call, +, 1, :x), 2) |> simplify! == Expr(:call, +, :x, -1)
    @test Expr(:call, +, Expr(:call, +, :x, 1), 2) |> simplify! == Expr(:call, +, :x, 3)
    @test Expr(:call, +, Expr(:call, +, 1, :x), 2) |> simplify! == Expr(:call, +, :x, 3)
    @test Expr(:call, +, Expr(:call, -, 1, :x), 2) |> simplify! == Expr(:call, -,  3, :x)
    @test Expr(:call, -, Expr(:call, -, 1, :x), 2) |> simplify! == Expr(:call, -, -1, :x)
    @test Expr(:call, +, Expr(:call, -, :x, 1), 2) |> simplify! == Expr(:call, +, :x, 1)
    @test Expr(:call, -, Expr(:call, -, :x, 1), 2) |> simplify! == Expr(:call, -, :x, 3)
    @test Expr(:call, +, :x, Expr(:call, -, 1, :x)) |> simplify! == 1
    @test Expr(:call, +, Expr(:call, -, 2, :x), :x) |> simplify! == 2
    @test Expr(:call, -, :x, Expr(:call, +, :x, :y)) |> simplify! == Expr(:call, -, :y)
    @test Expr(:call, -, :x, Expr(:call, -, :x, :y)) |> simplify! == :y
    @test Expr(:call, -, :x, Expr(:call, +, :y, :x)) |> simplify! == Expr(:call, -, :y)
    @test Expr(:call, -, Expr(:call, -, :x, :y), :x) |> simplify! == Expr(:call, -, :y)
    @test Expr(:call, -, Expr(:call, +, :x, :y), :x) |> simplify! == :y
    @test Expr(:call, +, 2, Expr(:call, +, 1, :x)) |> simplify! == Expr(:call, +, 3, :x)
    @test Expr(:call, +, 2, Expr(:call, -, 1, :x)) |> simplify! == Expr(:call, -, 3, :x)
    @test Expr(:call, +, 2, Expr(:call, +, :x, 1)) |> simplify! == Expr(:call, +, 3, :x)
    @test Expr(:call, +, 2, Expr(:call, -, :x, 1)) |> simplify! == Expr(:call, +, 1, :x)
    @test Expr(:call, -, 2, Expr(:call, +, 1, :x)) |> simplify! == Expr(:call, -, 1, :x)
    @test Expr(:call, -, 2, Expr(:call, +, :x, 1)) |> simplify! == Expr(:call, -, 1, :x)
    @test Expr(:call, -, 1, Expr(:call, -, 2, :x)) |> simplify! == Expr(:call, +, -1, :x)
    @test Expr(:call, -, 2, Expr(:call, -, :x, 1)) |> simplify! == Expr(:call, -, 3, :x)

    # evaluation
    ex = Expr(:call, +, 1, :x) |> Evolutionary.Expression
    xs = rand(rng, 10)
    @test ex(xs[1]) == xs[1]+1
    @test ex.(xs) == xs.+1
    io = IOBuffer()
    show(io, "text/latex", ex)
    @test String(take!(io)) == "\\left(1.0+x\\right)"
    Evolutionary.infix(io, ex)
    @test String(take!(io)) == "(1.0+x)"

    # protected operations
    @test Evolutionary.pdiv(1,0) == 1e7+1
    @test Evolutionary.aq(1,0) == 1.0
    @test Evolutionary.pexp(40) == 1e16+40
    @test Evolutionary.plog(0) == -1e7
    @test Evolutionary.psin(Inf) == 1e7
    @test Evolutionary.pcos(Inf) == 1e7
    @test Evolutionary.ppow(10,20) == 1e7+30
    @test Evolutionary.ppow(-4,0.5) == 2.0
    @test Evolutionary.ppow(2,2) == 4.0
    @test Evolutionary.cond(3,1,2) == 1
    @test Evolutionary.cond(-3,1,2) == 2

    # symreg
    fitfun(x) = x*x + x + 1.0
    rg = -5.0:0.1:5.0
    ys = fitfun.(rg)
    function fitobj(expr)
        ex = Evolutionary.Expression(expr)
        yy = ex.(rg)
        sum(v->isnan(v) ? 1.0 : v, abs2.(ys .- yy) )/length(rg)
    end

    Random.seed!(rng, 42)
    res = Evolutionary.optimize(fitobj,
        TreeGP(50, Terminal[:x, randn], Function[+,-,*,Evolutionary.aq],
            mindepth=1,
            maxdepth=3,
            simplify = Evolutionary.simplify!,
            selection = tournament(3),
        ),
        Evolutionary.Options(show_trace=false, rng=rng, iterations=50)
    )
    @test minimum(res) < 1.5

end

ext: jl


text:
ENV["JUPYTER"] = "F:/Anaconda/envs/ML/Scripts/jupyter.exe"
try
    @eval using Revise
catch e
    @warn "Error initializing Revise" exception = (e, catch_backtrace())
end


text:
module MVTs

import Base: @pure, @propagate_inbounds, @_pure_meta, @_inline_meta,
    @_propagate_inbounds_meta, get, getindex, setindex!, view, indices, size, length,
    summary, show, length, start, next, done, transpose, permutedims, ndims, show, print,
    endof, first, last, names

import Base: ==

export Label, @l_str, labels
export Row, Table

include("Label.jl")
include("symbolic.jl")

include("Relation.jl")
include("Row.jl")
include("Index.jl")
include("Table.jl")

include("projection.jl")
include("join.jl")
include("broadcast.jl")

include("isequal.jl")
#include("symbolic.jl")

end # module


text:
module Module1

using HTTP
using JSON

export getcoinprices

function getcoinprices(dateFrom::String,dateTo::String)
    # dateFrom > dateTo && return Dict()
    url = string("https://api.coindesk.com/v1/bpi/historical/close.json?currency=USD&start=",
                 dateFrom, "&end=", dateTo)
    res = HTTP.request("GET", url ,verbose=0)
    dat = JSON.parse(join(readlines(IOBuffer(res.body)), " "))
    haskey(dat, "bpi") ? dat["bpi"] : Dict()
end

end # module

max_stars_repo_path: Chapter08/01. Julia development workflow with Revise.jl/Module1.jl


text:
@testset "Multi-objective EA" begin

    rng = StableRNG(42)
    opts = Evolutionary.Options(rng=rng, iterations=500)

    # domination
    P =[ [0,0,0], [0,1,0], [-1,0,0], [0,-1,0] ]
    @test Evolutionary.dominate(P[1], P[2]) == 1
    @test Evolutionary.dominate(P[1], P[1]) == 0
    @test Evolutionary.dominate(P[1], P[3]) == -1
    @test Evolutionary.dominate(P[3], P[4]) == 0
    @test Evolutionary.dominations(P)[:, 1] == [0, -1, 1, 1]

    # convergence
    R = reshape([10,0,6,1,2,2,1,6,0,10],2,5)
    A = reshape([4,2,3,3,2,4], 2, 3)
    B = reshape([8,2,4,4,2,8], 2, 3)
    @test Evolutionary.igd(A,R) ≈ 3.707092031609239
    @test Evolutionary.igd(B,R) ≈ 2.59148346584763
    @test Evolutionary.spread(R) ≈ 0.0
    @test Evolutionary.spread(R[1:1,1:4]) ≈ 3.75

    # Schaffer F2

    schafferf2(x::AbstractVector) = [ x[1]^2,  (x[1]-2)^2 ]
    Random.seed!(rng, 1)
    result = Evolutionary.optimize(schafferf2, ()->100randn(rng,1), NSGA2(), opts)
    println("NSGA2:2RLT:SBX:PLM => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test isnan(Evolutionary.minimum(result))
    mvs = vcat(Evolutionary.minimizer(result)...)
    @test sum(0 .<= mvs .<= 2)/length(mvs) >= 0.8 # 80% in PO ∈ [0,2]
    #println(result)
    #println(extrema(mvs))

    function schafferf2!(F, x::AbstractVector) # in-place update
        F[1] = x[1]^2
        F[2] = (x[1]-2)^2
        F
    end
    Random.seed!(rng, 42)
    result = Evolutionary.optimize(schafferf2!, zeros(2), ()->10randn(rng,1), NSGA2(), opts)
    println("NSGA2:2RLT:SBX:PLM => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test isnan(Evolutionary.minimum(result))
    mvs = vcat(Evolutionary.minimizer(result)...)
    @test sum(0 .<= mvs .<= 2)/length(mvs) >= 0.8 # 80% in PO ∈ [0,2]

end


text:
# Pkg.clone("https://github.com/wookay/GlyphSheets.jl.git")
using GlyphSheets
import GlyphSheets.DSL: Move, Line, Cubic, Conic

library = create_library()
typeface = create_typeface(library, "Butler_Regular.otf")

glyph = create_glyph(typeface, 'B', 5pt)
outlines = outline_decompose(glyph)

close_typeface(typeface)
close_library(library)

using Cairo
using Colors
using ColorTypes

c = CairoRGBSurface(300, 300)
cr = CairoContext(c)

save(cr)
set_source_rgb(cr, 0.8, 0.8, 0.8)
rectangle(cr, 0, 0, 300, 300)
fill(cr)
restore(cr)

translate(cr, 0, glyph.slot.metrics.height+20)
scale(cr, 1, -1)

function operate(cr, op)
    if op isa Move
        move_to(cr, op.to...)
    elseif op isa Line
        line_to(cr, op.to...)
    elseif op isa Conic
        curve_to(cr, op.control..., op.to..., op.to...)
    elseif op isa Cubic
        curve_to(cr, op.control1..., op.control2..., op.to...)
    end
end

function fill_paths(cr, outlines)
    save(cr)
    operate.(cr, outlines)
    set_source_rgba(cr, 0.9, 1, 0.9, 1)
    fill(cr)
    restore(cr)
end

function stroke_paths(cr, outlines)
    save(cr)
    set_line_width(cr, 3)
    colors = distinguishable_colors(length(outlines))
    last_pos = (0, 0)
    for (color, op) in zip(colors, outlines)
        move_to(cr, last_pos...)
        operate(cr, op)
        set_source_rgba(cr, red(color), green(color), blue(color), alpha(color))
        stroke(cr)
        last_pos = op.to
    end
    restore(cr)
end

fill_paths(cr, outlines)
stroke_paths(cr, outlines)

write_to_png(c, "cairo.png")


text:
using HTTP

"""
    get_bow_indices(bow_key_or_path_list::Vector{String}, tokenizer)

Returns a list of `list of indices` of words from each Bag of word in the `bow_key_or_path_list`, after tokenization. The functions looks for provided BoW key in the registered artifacts `Artifacts.toml` file. In case not present there, function expects that bow_key is provided as the complete path to the file the URL to download .txt file.

# Example
``` julia
get_bow_indices(["legal", "military"])
```
"""
function get_bow_indices(bow_key_or_path_list::Vector{String}, tokenizer)
    bow_indices = []
    bow_words = []
    for bow_key_or_path in bow_key_or_path_list
        if !isnothing(artifact_hash(bow_key_or_path, ARTIFACTS_TOML))
            file_path = joinpath(get_registered_file(bow_key_or_path), string(bow_key_or_path, ".txt"))
        else
            if isfile(bow_key_or_path)
                file_path = bow_key_or_path
            elseif isurl(bow_key_or_path)
                isdir("BoW")||mkdir("BoW")
                HTTP.download(bow_key_or_path, joinpath("BoW", split(bow_key_or_path, "/")[end]))
            else
                error("$bow_key_or_path not a valid entry for BoW list source")
            end
        end
        f = open(file_path, "r")
        words = split(strip(read(f, String)), "\n")
        encoded = map(x -> tokenizer(x; add_prefix_space=true), words)
        push!(bow_indices, encoded)
        push!(bow_words, words)
        close(f)
    end
    return bow_indices, bow_words
end

"""
    build_bow_ohe(bow_indices, tokenizer)

Build and return a list of `one_hot_matrix` for each Bag Of Words list from indices. Each item of the list is of dimension `(num_of_BoW_words, tokenizer.vocab_size)`.

Note: While building the OHE of word indices, it only keeps those words, which have length `1` after tokenization and discard the rest.

"""
function build_bow_ohe(bow_indices, tokenizer)
    global device
    if bow_indices==nothing
        return
    end
    onehot =[]
    onehot_single = zeros(length(tokenizer.vocab))
    for bow in bow_indices
        bow = filter(x->(length(x)<=1),bow)
        x = zeros(length(bow), length(tokenizer.vocab))
        for (i, word_id) in enumerate(bow)
            x[i, word_id[1]] = 1
            onehot_single[word_id[1]]=1
        end
        x = x|> device
        push!(onehot, x)
    end
    onehot_single = onehot_single |> device
    return onehot_single, onehot
end


text:
"""
This is much like the exact burden, except for replacing the burden adjustment with
a similar adjustment. The original was `1 - prod(1 - b_i)`. This one is
`1 - exp(-sum(b_i))`. We try this approximation because this series, as a whole series,
can be integrated, term-by-term, in a way that will allow us to calculate it in a
closed form. You see the long calculation here, but, if this does well, then we
implement the much faster version.
"""
function exponential_burden_term(weights, prevalences, which)
    w = weights[which]
    yes = prevalences[which]
    no = prevalences[.!which]
    results = zeros(Float64, length(which))
    results[which] += w * ((1 - exp(-sum(w))) * prod(yes) * prod(1 .- no) / sum(w))
    results
end

function exponential_burden_term!(weights, prevalences, which, running_sum)
    w = weights[which]
    yes = prevalences[which]
    no = prevalences[.!which]
    running_sum[which] += w * ((1 - exp(-sum(w))) * prod(yes) * prod(1 .- no) / sum(w))
end

"""
    exponential_burden(weights, prevalences)

Calculate the burden using an alternative formula that has exponentials in it.
"""
function exponential_burden(weights, prevalences)
    n = length(prevalences)
    which = zeros(Bool, n)
    b = zeros(Float64, n)
    for i in 1:(1<<n - 1)
        to_flip = grey_next_flip(i)
        which[to_flip] = !which[to_flip]
        exponential_burden_term!(weights, prevalences, which, b)
    end
    b
end

"""
As the number of causes gets larger, the error of this method decreases.
It's a lot for 5 causes but around 1% for 20 causes. That's the maximum
relative error for any one burden.
"""
function exponential_burden_random(cnt)
    w = rand(cnt)
    p = rand(cnt)
    exact = exact_burden(w, p)
    exp_burden = exponential_burden(w, p)
    total = total_burden(w, p)
    trial = exp_burden * total / sum(exp_burden)
    [sum(exp_burden), total, sum(exact), maximum(abs.((exp_burden - exact) ./ exact))]
end


text:
# PNML Intermediate Representation
# Defines types instantiated on the upwards parse step using PnmlDicts
# created on the downward (leaf-ward) parsing step.
# Enables use of dispatch to create higher-level constructs.

#module IR

include("graphics.jl")
include("structure.jl")
include("toolinfos.jl")
include("common.jl")
include("labels.jl")
include("declarations.jl")
include("terms.jl")
include("sorts.jl")
include("markings.jl")
include("conditions.jl")
include("inscriptions.jl")
include("nodes.jl")
include("page.jl")
include("net.jl")
include("model.jl")

include("flatten.jl")
include("show.jl")

#end # module IR


text:
@doc raw"""
    gaussjacobi(n::Integer) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Jacobi quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature).

```math
\int_{-1}^{1} f(x) (1-x)^\alpha(1+x)^\beta dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gaussjacobi(3, 1/3, -1/3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 268π/729(√3)
true
```
"""
function gaussjacobi(n::Integer, α::Real, β::Real)
    #GAUSS-JACOBI QUADRATURE NODES AND WEIGHTS
    if n < 0
        throw(DomainError(n, "gaussjacobi($n,$α,$β) not defined: n must be non-negative."))
    elseif α == 0. && β == 0.
        gausslegendre(n)
    elseif α == -0.5 && β == -0.5
        gausschebyshev(n, 1)
    elseif α == 0.5 && β == 0.5
        gausschebyshev(n, 2)
    elseif α == -0.5 && β == 0.5
        gausschebyshev(n, 3)
    elseif α == 0.5 && β == -0.5
        gausschebyshev(n, 4)
    elseif n == 0
        Float64[], Float64[]
    elseif n == 1
        [(β - α) / (α + β + 2)], [2^(α + β + 1) * beta(α + 1, β + 1)]
    elseif min(α,β) ≤ -1.
        throw(DomainError((α,β), "The Jacobi parameters correspond to a nonintegrable weight function"))
    elseif n ≤ 100 && max(α,β) < 5.
        jacobi_rec(n, α, β)
    elseif n > 100 && max(α,β) < 5.
        jacobi_asy(n, α, β)
    elseif n ≤ 4000 && max(α,β) ≥ 5.
        jacobi_gw(n, α, β)
    else
        error("gaussjacobi($n,$α,$β) is not implemented: n must be ≤ 4000 for max(α,β) ≥ 5.")
    end
end

# Convenience function: convert any kind of numbers a and b to a joint floating point type
jacobi_rec(n::Integer, α::Real, β::Real) = jacobi_rec(n, promote(float(α), float(β))...)

function jacobi_rec(n::Integer, α::T, β::T) where {T <: AbstractFloat}
    #Compute nodes and weights using recurrrence relation.
    x11, x12 = HalfRec(n, α, β, 1)
    x21, x22 = HalfRec(n, β, α, 0)

    x = Array{T}(undef,n)
    w = Array{T}(undef,n)
    m1 = length(x11)
    m2 = length(x21)
    sum_w = zero(T)
    @inbounds for i in 1:m2
        idx = m2 + 1 - i
        xi = -x21[i]
        der = x22[i]
        wi = 1 / ((1 - xi^2) * der^2)
        w[idx] = wi
        x[idx] = xi
        sum_w += wi
    end
    @inbounds for i in 1:m1
        idx = m2 + i
        xi = x11[i]
        der = x12[i]
        wi = 1 / ((1 - xi^2) * der^2)
        w[idx] = wi
        x[idx] = xi
        sum_w += wi
    end
    c = (2^(α+β+1)*gamma(2+α)*gamma(2+β)/(gamma(2+α+β)*(α+1)*(β+1)))
    rmul!(w, c / sum_w)
    return x, w
end

function HalfRec(n::Integer, α::T, β::T, flag) where {T <: AbstractFloat}
    # HALFREC  Jacobi polynomial recurrence relation.
    # Asymptotic formula - only valid for positive x.
    r = (flag == 1) ? (ceil(n / 2):-1:1) : (floor(n / 2):-1:1)
    m = length(r)
    c1 = 1 / (2 * n + α + β + 1)
    a1 = one(T)/4 - α^2
    b1 = one(T)/4 - β^2
    c1² = c1^2
    x = Array{T}(undef,m)
    @inbounds for i in 1:m
        C = muladd(2*one(T), r[i], α - one(T)/2) * (T(π) * c1)
        C_2 = C/2
        x[i] = cos(muladd(c1², muladd(-b1, tan(C_2), a1 * cot(C_2)), C))
    end

    P1 = Array{T}(undef,m)
    P2 = Array{T}(undef,m)
    # Loop until convergence:
    for _ in 1:10
        innerjacobi_rec!(n, x, α, β, P1, P2)
        dx2 = 0.0
        @inbounds for i in 1:m
            dx = P1[i] / P2[i]
            _dx2 = abs2(dx)
            dx2 = ifelse(_dx2 > dx2, _dx2, dx2)
            x[i] = x[i] - dx
        end
        dx2 > eps(T) / 1e6 || break
    end
    # Once more for derivatives:
    innerjacobi_rec!(n, x, α, β, P1, P2)
    return x, P2
end

function innerjacobi_rec!(n, x, α::T, β::T, P, PP) where {T <: AbstractFloat}
    # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE.
    N = length(x)
    @inbounds for j = 1:N
        xj = x[j]

        Pj = (α - β + (α + β + 2) * xj)/2
        Pm1 = one(T)
        PPj = (α + β + 2)/2
        PPm1 = zero(T)
        for k in 1:n-1
            k0 = muladd(2*one(T), k, α + β)
            k1 = k0 + 1
            k2 = k0 + 2
            A = 2 * (k + 1) * (k + (α + β + 1)) * k0
            B = k1 * (α^2 - β^2)
            C = k0 * k1 * k2
            D = 2 * (k + α) * (k + β) * k2
            c1 = muladd(C, xj, B)
            Pm1, Pj = Pj, muladd(-D, Pm1, c1 * Pj) / A
            PPm1, PPj = PPj, muladd(c1, PPj, muladd(-D, PPm1, C * Pm1)) / A
        end
        P[j] = Pj
        PP[j] = PPj
    end
    nothing
end

function innerjacobi_rec(n, x, α::T, β::T) where {T <: AbstractFloat}
    # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE.
    N = length(x)
    P = Array{T}(undef,N)
    PP = Array{T}(undef,N)
    innerjacobi_rec!(n, x, α, β, P, PP)
    return P, PP
end

function weightsConstant(n, α, β)
    # Compute the constant for weights:
    M = min(20, n - 1)
    C = 1.0
    p = -α * β / n
    for m = 1:M
        C += p
        p *= -(m + α) * (m + β) / (m + 1) / (n - m)
        abs(p / C) < eps(Float64) / 100 && break
    end
    return 2^(α + β + 1) * C
end

function jacobi_asy(n, α, β)
    # ASY  Compute nodes and weights using asymptotic formulae.

    # Determine switch between interior and boundary regions:
    nbdy = 10
    bdyidx1 = n - (nbdy - 1):n
    bdyidx2 = nbdy:-1:1

    # Interior formula:
    x, w = asy1(n, α, β, nbdy)

    # Boundary formula (right):
    xbdy = boundary(n, α, β, nbdy)
    x[bdyidx1], w[bdyidx1] = xbdy

    # Boundary formula (left):
    if α ≠ β
        xbdy = boundary(n, β, α, nbdy)
    end
    x[bdyidx2] = -xbdy[1]
    w[bdyidx2] = xbdy[2]

    rmul!(w, weightsConstant(n, α, β))
    return x, w
end

function asy1(n::Integer, α::Float64, β::Float64, nbdy::Integer)
    # Algorithm for computing nodes and weights in the interior.

    # Approximate roots via asymptotic formula: (Gatteschi and Pittaluga, 1985)
    K = π*(2(n:-1:1).+α.-0.5)/(2n+α+β+1)
    tt = K .+ (1/(2n+α+β+1)^2).*((0.25-α^2).*cot.(K/2).-(0.25-β^2).*tan.(K/2))

    # First half (x > 0):
    t = tt[tt .≤ π/2]
    mint = t[end-nbdy+1]
    idx = 1:max(findfirst(t .< mint)-1, 1)

    # Newton iteration
    for _ in 1:10
        vals, ders = feval_asy1(n, α, β, t, idx)  # Evaluate
        dt = vals./ders
        t += dt  # Next iterate
        if norm(dt[idx],Inf) < sqrt(eps(Float64))/100
            break
        end
    end
    vals, ders = feval_asy1(n, α, β, t, idx)  # Once more for luck
    t += vals./ders

    # Store
    x_right = cos.(t)
    w_right = 1 ./ ders.^2

    # Second half (x < 0):
    α, β = β, α
    t = π .- tt[1:(n-length(x_right))]
    mint = t[nbdy]
    idx = max(findfirst(t .> mint), 1):length(t)

    # Newton iteration
    for _ in 1:10
        vals, ders = feval_asy1(n, α, β, t, idx)  # Evaluate.
        dt = vals./ders  # Newton update.
        t += dt
        if norm(dt[idx],Inf) < sqrt(eps(Float64))/100
            break
        end
    end
    vals, ders = feval_asy1(n, α, β, t, idx)  # Once more for luck.
    t += vals./ders  # Newton update.

    # Store
    x_left = cos.(t)
    w_left = 1 ./ ders.^2

    return vcat(-x_left, x_right), vcat(w_left, w_right)
end

"""
Evaluate the interior asymptotic formula at x = cos(t).
Assumption:
* `length(t) == n ÷ 2`
"""
function feval_asy1(n::Integer, α::Float64, β::Float64, t::AbstractVector, idx)
    # Number of terms in the expansion:
    M = 20

    # Number of elements in t:
    N = length(t)

    # Some often used vectors/matrices:
    onesM = ones(M)

    # The sine and cosine terms:
    A = repeat((2n+α+β).+(1:M),1,N).*repeat(t',M)/2 .- (α+1/2)*π/2  # M × N matrix
    cosA = cos.(A)
    sinA = sin.(A)

    sinT = repeat(sin.(t)',M)
    cosT = repeat(cos.(t)',M)
    cosA2 = cosA.*cosT .+ sinA.*sinT
    sinA2 = sinA.*cosT .- cosA.*sinT

    sinT = hcat(ones(N), cumprod(repeat((csc.(t/2)/2),1,M-1), dims=2))  # M × N matrix
    secT = sec.(t/2)/2

    _vec = [(α+j-1/2)*(-α+j-1/2)/(2n+α+β+j+1)/j for j in 1:M-1]
    P1 = [1;cumprod(_vec)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = Matrix(1.0I, M, M)
    for l in 1:M
        _vec = [(β+j-1/2)*(-β+j-1/2)/(2n+α+β+j+l)/j for j in 1:M-l-2]
        P2[l,l+1:M-2] = cumprod(_vec)
    end
    PHI = repeat(P1,1,M).*P2

    _vec = [(α+j-1/2)*(-α+j-1/2)/(2n+α+β+j-1)/j for j in 1:M-1]
    P1 = [1;cumprod(_vec)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = Matrix(1.0I, M, M)
    for l in 1:M
        _vec = [(β+j-1/2)*(-β+j-1/2)/(2n+α+β+j+l-2)/j for j in 1:M-l-2]
        P2[l,l+1:M-2] = cumprod(_vec)
    end
    PHI2 = repeat(P1,1,M).*P2

    S = zeros(N)
    S2 = zeros(N)
    for m in 1:M
        l = 1:2:m
        phi = PHI[l, m]
        dS1 = (sinT[:, l]*phi) .* cosA[m, :]
        phi2 = PHI2[l, m]
        dS12 = (sinT[:, l]*phi2) .* cosA2[m, :]
        l = 2:2:m
        phi = PHI[l, m]
        dS2 = (sinT[:, l]*phi) .* sinA[m, :]
        phi2 = PHI2[l, m]
        dS22 = (sinT[:, l]*phi2) .* sinA2[m, :]
        if m - 1 > 10 && norm(dS1[idx] + dS2[idx], Inf) < eps(Float64) / 100
            break
        end
        S .+= dS1
        S .+= dS2
        S2 .+= dS12
        S2 .+= dS22
        sinT[:,1:m] .*= secT
    end

    # Constant out the front:
    dsa = α^2/2n
    dsb = β^2/2n
    dsab = (α+β)^2/4n
    ds = dsa + dsb - dsab
    s = ds
    i = 1
    dsold = ds # to fix α = -β bug.
    while abs(ds/s)+dsold > eps(Float64)/10
        dsold = abs(ds/s)
        i += 1
        tmp = -(i-1)/(i+1)/n
        dsa = tmp*dsa*α
        dsb = tmp*dsb*β
        dsab = tmp*dsab*(α+β)/2
        ds = dsa + dsb - dsab
        s = s + ds
    end
    p2 = exp(s)*sqrt(2π*(n+α)*(n+β)/(2n+α+β))/(2n+α+β+1)
    # g is a vector of coefficients in
    # ``\Gamma(z) = \frac{z^{z-1/2}}{\exp(z)}\sqrt{2\pi} \left(\sum_{i} B_i z^{-i}\right)``, where B_{i-1} = g[i].
    # https://math.stackexchange.com/questions/1714423/what-is-the-pattern-of-the-stirling-series
    g = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880,
         5246819/75246796800, -534703531/902961561600,
         -4483131259/86684309913600, 432261921612371/514904800886784000]
    f(g,z) = dot(g, [1;cumprod(ones(9)./z)])

    # Float constant C, C2
    C = 2*p2*(f(g,n+α)*f(g,n+β)/f(g,2n+α+β))/π
    C2 = C*(α+β+2n)*(α+β+1+2n)/(4*(α+n)*(β+n))

    vals = C*S

    # Use relation for derivative:
    ders = (n*((α-β).-(2n+α+β)*cos.(t)).*vals .+ (2*(n+α)*(n+β)*C2).*S2)./(2n+α+β)./sin.(t)
    denom = 1 ./ (sin.(abs.(t)/2).^(α+0.5).*cos.(t/2).^(β+0.5))
    vals .*= denom
    ders .*= denom

    return vals, ders
end

function boundary(n::Integer, α::Float64, β::Float64, npts::Integer)
# Algorithm for computing nodes and weights near the boundary.

    # Use Newton iterations to find the first few Bessel roots:
    smallK = min(30, npts)
    jk = approx_besselroots(α, smallK)

    # Approximate roots via asymptotic formula: (see Olver 1974)
    phik = jk/(n + .5*(α + β + 1))
    x = cos.( phik .+ ((α^2-0.25).*(1 .-phik.*cot.(phik))./(8*phik) .- 0.25.*(α^2-β^2).*tan.(0.5.*phik))./(n + 0.5*(α + β + 1))^2 )

    # Newton iteration:
    for _ in 1:10
        vals, ders = innerjacobi_rec(n, x, α, β)  # Evaluate via asymptotic formula.
        dx = -vals./ders  # Newton update.
        x += dx  # Next iterate.
        if norm(dx,Inf) < sqrt(eps(Float64))/200
            break
        end
    end
    vals, ders = innerjacobi_rec(n, x, α, β)  # Evaluate via asymptotic formula.
    dx = -vals./ders
    x += dx

    # flip:
    x = reverse(x)
    ders = reverse(ders)

    # Revert to x-space:
    w = 1 ./ ((1 .- x.^2) .* ders.^2)
    return x, w
end

function jacobi_jacobimatrix(n, α, β)
    s = α + β
    ii = 2:n-1
    si = 2*ii .+ s
    aa = [(β - α)/(2 + s);
          (β^2 - α^2) ./ ((si .- 2).*si);
          (β^2 - α^2) ./ ((2n - 2+s).*(2n+s))]
    bb = [2*sqrt( (1 + α)*(1 + β)/(s + 3))/(s + 2) ;
          2 .*sqrt.(ii.*(ii .+ α).*(ii .+ β).*(ii .+ s)./(si.^2 .- 1))./si]
    return SymTridiagonal(aa, bb)
end

function jacobimoment(α, β)
    s = α + β
    T = float(typeof(s))
    # Same as 2^(α+β+1) * beta(α+1,β+1)
    return exp((s+1)*log(convert(T,2)) + loggamma(α+1)+loggamma(β+1)-loggamma(2+s))
end

function jacobi_gw(n::Integer, α, β)
    # Golub-Welsh for Gauss--Jacobi quadrature. This is used when max(α,β)>5.
    x, V = eigen(jacobi_jacobimatrix(n, α, β))  # Eigenvalue decomposition.
    # Quadrature weights:
    w = V[1,:].^2 .* jacobimoment(α, β)
    return x, w
end


text:
using JLD2
using DataFrames
using StatsPlots

@load "paper/box.jld2" box_df

p = @df box_df[.!ismissing.(box_df[:, :mean]), :] plot(:n, :mean, group=(:Type, :N), marker=4,
  xlabel="T", ylabel="mean coupling time", xticks=[500;1000;2000;4000], size=(400,250))
Plots.savefig(p, "paper/boxDemoPlot.pdf")

@load "paper/boxInd.jld2" boxInd_df

p = @df boxInd_df[.!ismissing.(boxInd_df[:, :mean]), :] plot(:n, :mean, group=(:Type, :N), marker=4,
  xlabel="T", ylabel="mean coupling time", xticks=[500;1000;2000;4000], size=(450,250))
Plots.savefig(p, "paper/boxDemoIndPlot.pdf")


text:
function gauss_quadrature(y::Array{<:Number},GLweights::AbstractArray{<:Real,1})
	# Assume that y is evaluated at the zeros
	# I = Σ wi*y(zi) where zi are the roots of Pl of the appropriate order
	return squeeze(sum(GLweights.*y,1),1) :: Array{<:Number}
end

function clenshaw_curtis_quadrature(y::Array{<:Real})::Array{<:Real}
	N = size(y,1)::Int64-1 # number of divisions

	g = (FFTW.r2r(y, FFTW.REDFT00, 1)/N) :: Array{<:Real}

	inty = zeros(size(y)[2:end]...) :: Array{Float64}

	for ind in eachindex(inty)
		remaining_indices = ind2sub(size(inty),ind)
		inty[ind] = g[1,remaining_indices...]
		for k=1:div(N+1,2)-1
			inty[ind] += 2g[2k+1,remaining_indices...]/(1-(2k)^2)
		end
	end
	return inty
end

function clenshaw_curtis_quadrature(y::Array{<:Complex})::Array{<:Complex}
	Ny = size(y,1)::Int64-1

	# compute the fft of the even, complex sequence
	slice = [1:size(y,dim) for dim in 2:ndims(y)]
	g = (fft(vcat(y,y[Ny:-1:2,slice...]),1)/2Ny ) :: Array{<:Complex}

	inty = zeros(Complex128,size(y)[2:end]...) :: Array{Complex128}

	for ind in eachindex(inty)
		remaining_indices = ind2sub(size(inty),ind)
		inty[ind] = 2g[1,remaining_indices...]
		for k=1:div(Ny+1,2)-1
			inty[ind] += 2(g[2k+1,remaining_indices...]+
				g[end-2k+1,remaining_indices...])/(1-(2k)^2)
		end
	end
	return inty
end

#####################################################################################
#####################################################################################
# Integral on and in sphere
#####################################################################################

function sphericalIntegral3D(y::Array{<:Number,3},
	r::Union{Real,AbstractArray{<:Real,1}}=1,
	θ::AbstractArray{<:Real,1}=linspace(0,π,size(y,2)),
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,3));
	θint::String="clenshaw_quadrature",even::String="avg",
	GLweights::Union{Void,AbstractArray{<:Real,1}} = nothing,
	ϕ2π::Bool=true):: Array{<:Number,0}

	if typeof(r) <: AbstractArray
		z = simps(r.^2 .* y,x=r,even=even) :: Array{<:Number,2}
	else
		dr = r
		z = dr^2 .* simps((1:size(y,1)).^2 .* y,dx=dr,even=even) :: Array{<:Number,2}
	end

	trapzSpherical2D(z,θ,ϕ,θint=θint,GLweights=GLweights,ϕ2π=ϕ2π) :: Array{<:Number,0}
end

function sphericalIntegral2D(y::Array{<:Number,2},
	θ::AbstractArray{<:Real,1}=linspace(0,π,size(y,1)),
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,2));
	θint::String="clenshaw_quadrature",
	GLweights::Union{Void,AbstractArray{<:Real,1}} = nothing,
	ϕ2π::Bool=true):: Array{<:Number,0}

	trapzSpherical2D(y,θ,ϕ,θint=θint,GLweights=GLweights,ϕ2π=ϕ2π)
end

function sphericalIntegral1D(y::Vector{<:Number},
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,1));
	ϕ2π::Bool=true):: Array{<:Number,0}

	trapzCircle(y,ϕ,ϕ2π=ϕ2π)
end

function sphericalIntegral(y::Array{<:Number,N},args...;kwargs...) where N
	@assert(N<=3,"Integrals above 3D are not defined")
	f = eval(parse("sphericalIntegral"*string(N)*"D"))
	return f(y,args...,kwargs...) :: Array{<:Number,0}
end

function cylindricalIntegral1D(y::Vector{<:Number},
	ϕ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,1));
	ϕ2π::Bool=true):: Array{<:Number,0}

	trapzCircle(y,ϕ,ϕ2π=ϕ2π)
end

function cylindricalIntegral2D(y::Array{<:Number,2},
	r::Union{Real,AbstractArray{<:Real,1}}=1,
	θ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,2));
	even::String="avg",ϕ2π::Bool=true):: Array{<:Number,0}

	if typeof(r) <: AbstractArray
		z = simps(r.* y,x=r,even=even) :: Vector{<:Number}
	else
		dr = r
		z = dr*simps((1:size(y,1)).* y,dx=dr,even=even) :: Vector{<:Number}
	end

	trapzCircle(z,θ,ϕ2π=ϕ2π) :: Array{<:Number,0}
end

function cylindricalIntegral3D(y::Array{<:Number,3},
	r::Union{Real,AbstractArray{<:Real,1}}=1,
	θ::AbstractArray{<:Real,1}=linspace(0,2π,size(y,2)),
	z::Union{Real,AbstractArray{<:Real,1}}=1;
	even::String="avg",ϕ2π::Bool=true):: Array{<:Number,0}

	if typeof(r) <: AbstractArray
		yθz = simps(r.* y,x=r,even=even) :: Array{<:Number,2}
	else
		dr = r
		yθz = dr*simps((1:size(y,1)).* y,dx=dr,even=even) :: Array{<:Number,2}
	end

	T = promote_type(eltype(y), eltype(θ))
	yz = zeros(T,size(y,3))  :: Array{<:Number,1}

	for ind in eachindex(yz)
		yz[ind] = trapzCircle(yθz[1:end,ind],θ,ϕ2π=ϕ2π)[1]
	end

	if typeof(z) <: AbstractArray
		y_int = simps(z.* yz,x=z,even=even) :: Array{<:Number,0}
	else
		dz = z
		y_int = dz*simps((1:size(yz,1)).* yz,dx=dz,even=even) :: Array{<:Number,0}
	end

	return y_int :: Array{<:Number,0}
end

function cylindricalIntegral(y::Array{<:Number,N},args...;kwargs...) where N
	@assert(N<=3,"Integrals above 3D are not defined")
	f = eval(parse("cylindricalIntegral"*string(N)*"D"))
	return f(y,args...,kwargs...) :: Array{<:Number,0}
end


text:
function render!(rendermodel::RenderModel, roadway::StraightRoadway;
    color_asphalt::Colorant=COLOR_ASPHALT,
    lane_width::Float64 = DEFAULT_LANE_WIDTH,
    extra_length::Float64 = 50.0, # [m]
    lane_marking_width::Float64 = 0.15, # [m]
    )

    pts = Array(VecE2, 2)
    pts[1] = VecE2(-extra_length, 0)
    pts[2] = VecE2( extra_length + roadway.length, 0)

    add_instruction!(rendermodel, render_line, (pts, color_asphalt, lane_width))
    add_instruction!(rendermodel, render_line, ([p + VecE2(0, -lane_width/2) for p in pts], COLOR_LANE_MARKINGS_WHITE, lane_marking_width))
    add_instruction!(rendermodel, render_line, ([p + VecE2(0,  lane_width/2) for p in pts], COLOR_LANE_MARKINGS_WHITE, lane_marking_width))
    return rendermodel
end
function render(roadway::StraightRoadway;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel = RenderModel(),
    cam::Camera = FitToContentCamera(),
    )

    s = CairoRGBSurface(canvas_width, canvas_height)
    ctx = creategc(s)
    clear_setup!(rendermodel)

    render!(rendermodel, roadway)

    camera_set!(rendermodel, cam, canvas_width, canvas_height)
    render(rendermodel, ctx, canvas_width, canvas_height)
    return s
end
Base.show(io::IO, ::MIME"image/png", roadway::StraightRoadway) = show(io, MIME"image/png"(), render(roadway))


text:
include("EventsWiring.jl")

using Cairo, Gtk, Gtk.ShortNames, Graphics

export AttatchEvents

function AttatchEvents(document, canvas)
  #= Example of keypress output:
Gtk.GdkEventKey(
      event_type::GEnum          8,
      gdk_window::Ptr{Void}      Ptr{Void} @0x0000000006f761d0,
      send_event::Int8           0,
      time::UInt32               0x2bd9c220,
      state::UInt32              0x00000000,
      keyval::UInt32              0x00000061,
      length::Int32              1,
      string::Ptr{UInt8}         Ptr{UInt8} @0x000000001d32c310,
      hardware_keycode::UInt16   0x0041,
      group::UInt8               0x00,
      flags::UInt32              0x00000000
)
=#
# for signals: https://developer.gnome.org/gtk3/stable/GtkWidget.html#GtkWidget-accel-closures-changed
    signal_connect(document.win, "key-press-event") do widget, event
      println("You pressed key ", event) #.keyval
    end
    signal_connect(document.win, "key-release-event") do widget, event
      println("You pressed key ", event) #.keyval
    end
    #.............................................................
    canvas.mouse.button1motion = @guarded (widget, event) -> begin
        #DrawEvent(document, widget, event)
        push!(document.event.pressed, Point(event.x, event.y) )
        #println("move x$(event.x) y$(event.y)")
        #println.(document.event.pressed)
        #DragEvent(document, widget, event)
        pressed, released = document.event.pressed[1], document.event.pressed[end]
        ctx = getgc(widget)
        selectText(ctx, document, pressed, released)
        reveal(widget)
    end
    #.............................................................
    canvas.mouse.motion = @guarded (widget, event) -> begin
        # mouseover #####################################
        hovers = document.eventsList.mouseover
        for node in hovers
            if onArea( event.x, event.y, node.shape) && insideParentClip(event.x, event.y, node) && node !== document.focusNode
                ctx = getgc(widget)
                page = document.children[1].children[3]
                node.temp = deepcopy(node.shape)

                #"hover":{"color":"yellow"},
                if haskey(node.DOM, "hover")
                        hover = node.DOM["hover"]
                        if haskey(hover, "color")
                            node.shape.color = hover["color"]
                                #node.shape.color = [0.996094, 0.5, 0.5]
                        end
                end

                drawNode(ctx, document, node)
                reveal(widget)
                # (1) register hovered object
                document.focusNode = node
            end
        end
        # mousemove #####################################
        # hovers = document.eventsList.mousemove
        # for node in hovers
        #     if onArea(node.shape, event.x, event.y)
        #         println("X: $(event.x), Y: $(event.y), ")
        #     end
        # end
        # mouseout #####################################
        # (1) get registered hover object
         #drawNode(ctx, document, document.focusNode)
         #reveal(widget)
        # # () unregister un-hovered object
        node = document.focusNode
        if node != nothing
                if !onArea( event.x, event.y, node.shape) # && insideParentClip( node , event.x, event.y)
                        if haskey(node.DOM, "hover")
                                hover = node.DOM["hover"]
                                if haskey(hover, "color")
                                    if haskey(node.DOM, "color")
                                           node.shape.color = GetTheColor(node.shape, node.DOM["color"])
                                    else
                                           node.shape.color = inheritColor(node.parent)
                                    end
                                end
                        end
                        ctx = getgc(widget)
                        drawNode(ctx, document, node)
                        reveal(widget)
                        document.focusNode = document.parent
                end
        end

        #  unhovers = document.eventsList.mouseout
        #  found = false
        #  for n in unhovers
        #          if node !== n
        #          end
        #
        # #     if onArea( event.x, event.y,node.shape)
        # #         println("X: $(event.x), Y: $(event.y), ")
        # #     end
    end
    # Record to later test for click event
    #.............................................................
    canvas.mouse.button1press = @guarded (widget, event) -> begin
        document.event.pressed = [] # Test code to be removed!
        document.event.pressed = [Point(event.x, event.y)]
        #push!(document.event.pressed, Point(event.x, event.y) )
        MouseDownEvent(document, widget, event)
    end

    #.............................................................
    canvas.mouse.button1release = @guarded (widget, event) -> begin
        push!(document.event.pressed, Point(event.x, event.y) )
        pressed, released = document.event.pressed[1], document.event.pressed[end]
        # Click Event fired!
           if -5 < (pressed.x - released.x) < 5 && -5 < (pressed.y - released.y) < 5
             ClickEvent(document, widget, event)
           end
           document.event.pressed = []
    end
    #.............................................................
    canvas.mouse.scroll = @guarded (widget, event) -> begin
        ScrollEvent(document, widget, event)
    end

end

#======================================================================================#
#    Similar to MoveAll() in LayoutBuild.jl
#======================================================================================#
function VmoveAll(node,y)
  shape = node.shape
  if shape.flags[Fixed]
    return
  end
  shape.top  += y
    if isdefined(node, :rows) # ..it has rows of children so let's move them!
      for i in 1:length(node.rows)
        row = node.rows[i]
        row.y += y
        for j in 1:length(row.nodes)
            VmoveAll(row.nodes[j],y) # do the same for each child
        end
      end
    end
end


text:
using Base.Test
using SapphireORM

# Testing constructor
global conn = Sapphire(dsn=dsnTest,username=root,password="")

max_stars_repo_name: PhoenixMachina/SapphireORM

max_issues_repo_name: PhoenixMachina/SapphireORM

max_forks_repo_name: PhoenixMachina/SapphireORM


text:
struct TPU_SGD
    # Learning rate; the only data this optimizer needs to bundle with itself
    η::XRTArray{Float32,(),0}
end

# Simplest update step in existence.
sgd_update!(model::XRTArray, Δ::XRTArray, η) = model - (Δ .* η)

# If this leaf node had no updates calculated for it, then skip out early.
sgd_update!(model, Δ::Nothing, η) = model

function sgd_update!(model, Δ, η)
    # Base condition; if we have reached a leaf node return the inputs unchanged.
    # Note that if `model` is an XRTArray, we will hit the override above that actually
    # updates the model rather than this generic update!(), same for if Δ is `nothing`.
    if nfields(model) == 0
        return model
    end

    # Recursively pass the fields of this model through the update machinery.  We use
    # this strange ntuple() do-block because we cannot perform any kind of mutation
    # (such as push!()'ing onto a list) and so we adopt this more functional-style of
    # programming.
    new_fields = ntuple(Val(nfields(model))) do i
        return sgd_update!(getfield(model, i), getfield(Δ, i), η)
    end

    # Return something of the same type as `model`, but with the new fields
    if isa(model, Tuple)
        return tuple(new_fields...)
    else
        return typeof(model)(new_fields...)
    end
end

# Main entry point for this optimizer's update steps
update!(opt::TPU_SGD, model, Δ) = sgd_update!(model, Δ, opt.η)

struct TPU_ADAM{T<:Tuple}
    # This `state` is a recursive named tuple that mimicks the structure
    # of the actual model weights, only for each parameter, this `state`
    # holds a tuple of state for that weight (mt, vt, βp)
    state::T
    η::XRTArray{Float32,(),0}
    β::Tuple{XRTArray{Float32,(),0},XRTArray{Float32,(),0}}
end

function TPU_ADAM(model::ImmutableChain, η, β)
    # Take in the model, flatten it into a tuple and map each parameter to (w1, w2, β)
    flat_model = XLA.flatten_tuple(model)
    opt_state = Zygote.map(x -> (zero(x), zero(x), β), flat_model)

    # We unflatten so that the optimizer state retains the same structure as the model,
    # so that our update!() function can simply walk the
    opt_state = unflatten_opt_state(model, opt_state)
    return TPU_ADAM(opt_state, η, β)
end

function adam_update!(state::Tuple, model::XRTArray, Δ::XRTArray, η, β)
    # Destructure the state stored within `state`
    mt, vt, βp = state

    # Update `mt` and `vt` with statistics of `Δ`
    mt = β[1] .* mt .+ (XRTArray(1f0) .- β[1]) .* Δ
    vt = β[2] .* vt .+ (XRTArray(1f0) .- β[2]) .* Δ.^XRTArray(2f0)

    # Package these statistics up as the new state
    new_state = (mt, vt, βp .* β)

    # Calculate new weights by taking the statistics and subtracting them from the model
    vt_denom = sqrt.(vt ./ (XRTArray(1f0) .- βp[2])) .+ XRTArray(eps(Float32))
    Δ_model = mt ./ ((XRTArray(1f0) .- βp[1]) .* vt_denom)
    new_model = model .- η .* Δ_model

    # Return a tuple of updated (state, model)
    return (new_state, new_model)
end

# If this leaf node had no updates calculated for it, then skip out early.
adam_update!(state::Tuple, model, Δ::Nothing, η, β) = (state, model)

function adam_update!(state::Tuple, model, Δ, η, β)
    # Base condition; if we have reached a leaf node return the inputs unchanged.
    # Note that if `model` is an XRTArray, we will hit the override above that actually
    # updates the model rather than this generic update!(), same for if Δ is `nothing`.
    if nfields(model) == 0
        return (state, model)
    end

    # Recursively call update!(), delving into the state, model and updates,
    # which should have very similar structures.
    new_fields_and_state = ntuple(Val(nfields(model))) do i
        return adam_update!(getfield(state, i), getfield(model, i), getfield(Δ, i), η, β)
    end

    # Destructure combined (new_state, new_weights) return value from above.  We use
    # this strange ntuple() do-block syntax because we cannot push!() elements onto
    # lists and create tuples from them, we must program in a more functional style,
    # honoring the immutable structure of the program.
    new_state = ntuple(Val(nfields(model))) do i
        return new_fields_and_state[i][1]
    end
    new_fields = ntuple(Val(nfields(model))) do i
        return new_fields_and_state[i][2]
    end

    # Construct tuples that contain the model and state that we have updated.
    if isa(model, Tuple)
        return (new_state, new_fields)
    else
        return (tuple(new_state...), typeof(model)(new_fields...))
    end
end

# Entry point for optimizer update step; immediately start recursing over
# the values held within our optimizer state, and return (new_state, new_model)
function update!(opt::TPU_ADAM, model::ImmutableChain, Δ)
    new_opt_state, model = adam_update!(opt.state, model, Δ, opt.η, opt.β)
    return TPU_ADAM(new_opt_state, opt.η, opt.β), model
end

# Helper function to unflatten the optimizer state.  Very similar to unflatten_tuple(),
# but emits tuples and recurses slightly differently.
@generated function unflatten_opt_state(pattern, t)
    function expand_arg(pattern, idx)
        fields = Any[]
        for i = 1:fieldcount(pattern)
            fT = fieldtype(pattern, i)
            if fT === Nothing
                push!(fields, nothing)
                continue
            end
            if !(fT <: AbstractArray)
                arg, idx = expand_arg(fT, idx)
                push!(fields, arg)
            else
                push!(fields, Expr(:call, :getfield, :t, idx))
                idx += 1
            end
        end
        Expr(:call, tuple, fields...), idx
    end
    res = expand_arg(pattern, 1)[1]
    return res
end


text:
abstract type Point3D{T} <: FieldVector{3,T} end

struct SPoint3D{T} <: Point3D{T}
    x::T
    y::T
    z::T
end

@inline function SPoint3D{T}(A::AbsMat, colidx::Integer) where {T}
    @boundscheck _checkbounds_point3d(A, colidx)
    @inbounds SPoint3D{T}(A[1, colidx], A[2, colidx], A[3, colidx])
end
@inline SPoint3D(A::AbsMat, colidx::Integer) = SPoint3D{eltype(A)}(A, colidx)

mutable struct MPoint3D{T} <: Point3D{T}
    x::T
    y::T
    z::T
end

@inline function MPoint3D{T}(A::AbsMat, colidx::Integer) where {T}
    @boundscheck _checkbounds_point3d(A, colidx)
    @inbounds MPoint3D{T}(A[1, colidx], A[2, colidx], A[3, colidx])
end
@inline MPoint3D(A::AbsMat, colidx::Integer) = MPoint3D{eltype(A)}(A, colidx)

@inline function _checkbounds_point3d(A::AbsMat, colidx::Integer)
    axes(A, 1) === Base.OneTo(3) || throw(ArgumentError("`A` must be a (3, N) matrix"))
    if !checkbounds(Bool, axes(A, 2), colidx)
        Base.throw_boundserror(A, (Base.OneTo(3), colidx))
    end
    nothing
end


text:
using Random
include(joinpath(@__DIR__, "pedagogy_spatial.jl"))

β₀ = 0.1
r = 0.05  # interaction range = 2*radius
β, γ = β₀ / (π * (r/2)^2), 0.1  # γ = 0 ≡ Hard core

β₀ = 0.3
r = 0.07  # interaction range = 2*radius
β, γ = β₀ / (π * (r/2)^2), 0.1  # γ = 0 ≡ Hard core

c, w = [0.0, 0.0], 1.0
win = PRS.SquareWindow(c, w)

strauss = PRS.StraussPointProcess(β, γ, r, win)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_grid_prs(rng, strauss)

p = pedagogy_plot(sample,
    true,
    strauss.r/2,
    "white",
    win
)

path = joinpath(@__DIR__, "output", "strauss", "strauss_grid_prs.pdf")
Plots.savefig(p, path)


text:
@testset "experiment.jl" begin

    include("launch_experiment.jl")
    include("training.jl")
    include("benchmark_solving.jl")
    include("evaluation.jl")
    include("basicmetrics.jl")
end


text:
@everywhere using CSV, GLMNet, Combinatorics, DataFrames, Plots
#@everywhere using ScikitLearn: fit!, predict, @sk_import, fit_transform!
#@everywhere @sk_import ensemble: RandomForestRegressor

include("../forc-helpers/OLSestimator.jl")
include("../forc-helpers/create_matrices.jl")
include("../forc-models/all_models.jl")

function prepare(level, h, J)
    h_lag = 12 ## Number of months over which to calculate inflation
    J_cv = (J - (1 + h_lag + h)) ÷ 2 + (1 + h_lag + h) ## Half of the within sample size used as for cross validation
    head_line = true ## Choose the head line price to construct the inflation rates
    println("h = $h, J = $J, J_cv = $J_cv")
    Y, Y_lag, X_lag, X_lag2 = create_matrices(level, h, h_lag, head_line)
    return h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
function bicgstab(a::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64};
  tol=1.0e-10, limit=1000, trace=false)

  n = size(a, 1)
  x = copy(x0)
  local converged

  r = b - a * x
  r0_hat = copy(r)
  rho0 = w = α = 1.0
  v = zeros(n)
  p = zeros(n)
  rho1 = dot(r0_hat, r)
  trace && println("\nFirst few iterations:")
  iters = 0
  while true
    iters += 1
    converged = norm(r) < tol * norm(b)
    (converged || (iters == limit)) && break
    β = (rho1 / rho0) * (α / w)
    p = r + β * (p - w * v)
    v = a * p
    α = rho1 / dot(r0_hat, v)
    s = r - α * v
    t = a * s
    w = dot(t, s) / dot(t, t)
    rho0 = rho1
    rho1 = -w * dot(r0_hat, t)
    x += α * p + w * s
    r = s - w * t
    if trace && iters <= 5
      println("After $iters iterations: x = $x")
    end
  end
  iters, converged, x
end


text:
# optimization with potential
using ColorSchemes, Revise, FourierGPE

import FourierGPE:Lgp!, nlin!, V
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

# prettier plot
import FourierGPE.showpsi

bone = cgrad(ColorSchemes.bone_1.colors)
turbo = cgrad(ColorSchemes.turbo.colors)
function showpsi(x,y,ψ)
    p1 = heatmap(x,y,abs2.(ψ),aspectratio=1,c=bone)
    xlims!(x[1],x[end]);ylims!(y[1],y[end])
    xlabel!(L"x");ylabel!(L"y")
    title!(L"|\psi|^2")
    p2 = heatmap(x,y,angle.(ψ),aspectratio=1,c=turbo)
    xlims!(x[1],x[end]);ylims!(y[1],y[end])
    xlabel!(L"x");ylabel!(L"y")
    title!(L"\textrm{phase} (\psi)")
    p = plot(p1,p2,size=(600,300))
    return p
end

function Lgp!(dϕ,ϕ,sim,t)
    @unpack γ,μ,espec = sim
    nlin!(dϕ,ϕ,sim,t)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
    return nothing
end

function nlin!(dϕ,ϕ,sim::Sim{2},t)
    @unpack g,X,V0 = sim; x,y = X
    y = y'
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,y,t) + g*abs2(dϕ)
    # @. dϕ *= V(x,y,t) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

# ==== Initialize simulation
L = (25.0,25.0)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

# ==== set simulation parameters
μ = 15.0

# Time-dependent potential function (here trivial t dep)
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# Or, as a static Potential
# x,y = X
# Vs(x,y) = 0.5*(x^2 + y^2)
# ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-Vs(x,y)/μ,0.0)+im*0.0)
# sim = Sim(sim,V0=Vs.(x,y'))

# ==== make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

# ==== evolve
@time sol = runsim(sim)

# ==== ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

# ==== set simulation parameters
γ = 0.0
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψg,sim)
# reltol = 1e-7
# alg = DP5()

using VortexDistributions
R(w) = sqrt(2*μ/w^2)
R(1)
rv = 3.
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
ξ0 = healinglength.(0.,0.,μ,g)
ξ = healinglength(rv,0.,μ,g)

# vcore = Exact(VortexDistributions.ψi,ξ)
vcore = Exact(ξ)
pv = PointVortex(rv,0.,1)
vi = ScalarVortex(vcore,pv)
psi = Torus(copy(ψg),x,y)
vortex!(psi,vi)
showpsi(x,y,psi.ψ)
ψi .= psi.ψ
ϕi = kspace(ψi,sim)

# ===== compare with Fetter JLTP 2010
ξ = 1/sqrt(μ)
Rtf = R(1)
Ωm = 3*log(Rtf/ξ/sqrt(2))/2/Rtf^2
Ωv = Ωm/(1-rv^2/Rtf^2)

# or a precession period
Tv = 2*pi/Ωv

ti = 0.; tf = Tv
t = LinRange(ti,tf,Nt)

@pack_Sim! sim

# ==== evolve
@time solv = runsim(sim)

# ==== analyse
ϕf = solv[end]
ψf = xspace(ϕf,sim)
showpsi(x,y,ψf)

anim = @animate for i in eachindex(t)
    ψ = xspace(solv[i],sim)
    showpsi(x,y,ψ)
end

saveto = joinpath(@__DIR__,"vortexprecession.gif")
gif(anim,saveto,fps=20)


text:
function load(fs::Stream{format"IFS"}, MeshType = GLNormalMesh)
    io = stream(fs)
    function str()
        n = read(io, UInt32)
        String(read(io, UInt8, n))
    end
    name = str() # just skip name for now
    vertices = str()
    if vertices != "VERTICES\0"
        error("$(filename(fs)) does not seem to be of format IFS")
    end
    nverts = read(io, UInt32)
    verts_float = read(io, Float32, nverts * 3)
    verts = reinterpret(Point3f0, verts_float)
    tris = str()
    if tris != "TRIANGLES\0"
        error("$(filename(fs)) does not seem to be of format IFS")
    end
    nfaces = read(io, UInt32)
    faces_int = read(io, UInt32, nfaces * 3)
    faces = reinterpret(GLTriangle, faces_int)
    MeshType(vertices = verts, faces = faces)
end

function save(fs::Stream{format"IFS"}, msh::AbstractMesh; meshname = "mesh")
    io = stream(fs)
    function write0str(s)
        s0 = s * "\0"
        write(io, UInt32(length(s0)))
        write(io, s0)
    end
    vts = decompose(Point3f0, msh)
    fcs = decompose(GLTriangle, msh)

    # write the header
    write0str("IFS")
    write(io, 1f0)
    write0str(meshname)

    write0str("VERTICES")
    write(io, UInt32(length(vts)))
    write(io, reinterpret(Float32, vts))

    write0str("TRIANGLES")
    write(io, UInt32(length(fcs)))
    write(io, reinterpret(UInt32, fcs))
end

max_stars_repo_name: JuliaIO/MeshIO

max_issues_repo_name: JuliaIO/MeshIO

max_forks_repo_name: JuliaIO/MeshIO


text:
"""
nested_eltype(x)

returns the element type of a variable
on numbers, returns `typeof(x)`
on arrays, works recursively (first element) until reaching an array of numbers, and then returns `eltype(x)`
"""
function nested_eltype end

nested_eltype(x::Number) = typeof(x)
nested_eltype(x::AbstractArray{<:Number}) = eltype(x)
nested_eltype(x::AbstractArray) = nested_eltype(first(x))

"""
nested_copyto!(dest,src)

assigns the contents of `src` on `dest`, works on deeply nested arrays

"""
function nested_copyto! end

function nested_copyto!(dest::Number,src::Number)
return dest
end

function nested_copyto!(dest::T,src::T) where T<:AbstractArray{<:Number}
copyto!(dest,src)
return dest
end

function nested_copyto!(dest::T,src::T) where T<:AbstractArray
nested_copyto!.(dest,src)
return dest
end


text:
#!/usr/bin/env julia

module ResEval

export readtodict, sumdict, mindict, argmindict, meanstats

function readfix(filename)
    a = open(filename) do f
        readlines(f)[2:end]
    end
    return a
end

function readtodict(filename, genfirst=false)
    mixed = readfix(filename)
    mixed = map(x -> split(x, ','), mixed)
    mixed = collect(Iterators.flatten(mixed))
    mixed = reshape(mixed, 2, :)

    combined = Array{Union{String, Float64}, 2}(undef, 2, size(mixed, 2))
    combined[1, :] = convert(Array{String, 1}, @view mixed[1, :])
    combined[2, :] = map(x -> parse(Float64, x), @view mixed[2, :])

    dict = Dict{Int16, Dict{Char, Dict{Char, Float64}}}()
    name_re = r"(\d+)-([a-z])-([a-z])-.*"
    for (n, v) in eachcol(combined)
        m = match(name_re, n)
        @assert m !== nothing "$n is a $(typeof(n))"
        imsize = parse(Int16, m.captures[1])
        discnorm, gennorm = map(x -> Char(x[1]), m.captures[2:3])
        if genfirst
            discnorm, gennorm = gennorm, discnorm
        end

        if !haskey(dict, imsize)
            dict[imsize] = Dict{Char, Dict{Char, Float64}}()
        end
        if !haskey(dict[imsize], discnorm)
            dict[imsize][discnorm] = Dict{Char, Float64}()
        end

        dict[imsize][discnorm][gennorm] = v
    end
    return dict
end

function sumdict(dict::Dict{Char, Float64})
    reduce(+, values(dict))
end

function sumdict(dict::Dict)
    mapreduce(sumdict, +, values(dict))
end

function mindict(dict::Dict{Char, Float64})
    minimum(values(dict))
end

function mindict(dict::Dict)
    minimum(map(mindict, values(dict)))
end

function argmindict(dict::Dict{Char, Float64})
    argmin(dict)
end

function argmindict(dict::Dict)
    minargs = Dict(k * ' ' * argmindict(v) => minimum(values(v)) for (k, v) in dict)
    res = argmin(minargs)
    println("min: ", minargs[res], " at ", res)
    return res
end

function getcombmeans(a, b)
    means = Vector{Float64}()
    i = 0
    for (c, d) in zip(a, b)
        cs, cv = split(c, ',')[1:2]
        ds, dv = split(d, ',')[1:2]
        cv = parse(Float64, cv)
        dv = parse(Float64, dv)
        i += 1
        mean = (cv + dv) / 2
        println(i, ' ', cs, "\t\t", mean)
        push!(means, mean)
    end
    println("min: ", minimum(means), '\t', argmin(means))
    println("max: ", maximum(means), '\t', argmax(means))
end

function meanstats()
    a = readfix("ganerator-64-results.csv")
    b = readfix("ganerator-128-results.csv")
    getcombmeans(a, b)
end

end

max_stars_repo_name: janEbert/GANerator

max_issues_repo_name: janEbert/GANerator

max_forks_repo_name: janEbert/GANerator


text:
using LabelledArrays
using Test
using StaticArrays
using InteractiveUtils
using ChainRulesTestUtils

@time begin
@time @testset "SLArrays" begin include("slarrays.jl") end
@time @testset "LArrays" begin include("larrays.jl") end
@time @testset "DiffEq" begin include("diffeq.jl") end
@time @testset "ChainRules" begin include("chainrules.jl") end
end


text:
@testset "Generic elliptic curve" begin

  @testset "Constructors" begin
    @test_throws ErrorException EllipticCurve([1])
    @test_throws ErrorException EllipticCurve([1, 2, 3])
    @test_throws ErrorException EllipticCurve([1, 2, 3, 4, 5, 6])
    @test_throws ErrorException EllipticCurve([0, 0])
    @test_throws ErrorException EllipticCurve([0, 0, 0, 0, 0])

    E = @inferred EllipticCurve([1, 2], false)
    @test typeof(E) == EllCrv{fmpq}

    E = @inferred EllipticCurve([1, 2, 3, 4, 5])
    @test typeof(E) == EllCrv{fmpq}

    # this is Cremona: 11a2, lmfdb: 11.a1
    E11_a1 = @inferred EllipticCurve([0, -1, 1, -7820, -263580], false)

    # this is Cremona: 41a1, lmfdb: 41.a1
    E43_a1 = @inferred EllipticCurve([0, 1, 1, 0, 0])

    Qx, x = PolynomialRing(FlintQQ, "x")
    K, a = NumberField(x^2 - x - 1, "a")
    OK = maximal_order(K)

    E31_1_a1 = @inferred EllipticCurve([K(1), a + 1, a, a, K(0)])
    @test typeof(E31_1_a1) == EllCrv{nf_elem}

    # lmfdb: 116.1-a1
    E116_1_a1 = @inferred EllipticCurve([K(1), K(-1), a, -a, K(0)] )
    @test typeof(E31_1_a1) == EllCrv{nf_elem}

    # short example
    Eshort = @inferred EllipticCurve([4, 0])
    @test typeof(Eshort) == EllCrv{fmpq}
  end

  # The following curves will be used in later tests
  # Creation of these was tested in previous testset
  E11_a1 = EllipticCurve([0, -1, 1, -7820, -263580], false)

  E43_a1 = EllipticCurve([0, 1, 1, 0, 0])

  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - x - 1, "a")
  OK = maximal_order(K)

  E31_1_a1 = EllipticCurve([K(1), a + 1, a, a, K(0)])

  E116_1_a1 =EllipticCurve([K(1), K(-1), a, -a, K(0)] )

  Eshort = EllipticCurve([4, 0])

  @testset "Field access" begin
    @test base_field(E11_a1) == FlintQQ
    @test base_field(E43_a1) == FlintQQ
    @test base_field(E31_1_a1) == K
    @test base_field(E116_1_a1) == K
    @test base_field(Eshort) == FlintQQ
  end

  @testset "Weierstraß model computation" begin
    E = EllipticCurve([1,2,3,4,5])
    EE, f, g = @inferred short_weierstrass_model(E)
    @test isshort(EE)
    @test EE.coeff == [fmpq(61, 16), fmpq(127, 32)]
    P = E([1, 2])
    @test P == g(f(P))

    R = ResidueRing(FlintZZ, 5)
    E = EllipticCurve(map(R, [1, 2, 3, 4, 5]))
    EE, f, g = @inferred short_weierstrass_model(E)
    @test isshort(EE)
    @test EE.coeff == [R(1), R(1)]
    P = rand(EE)
    @test P == f(g(P))
    # @inferred will break the tests
  end

  @testset "Point construction" begin
    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred isfinite(P)
    @test @inferred !isinfinite(P)

    P = @inferred E43_a1([-1, 0], false)
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred isfinite(P)
    @test @inferred !isinfinite(P)

    P = @inferred E43_a1([fmpz(-1), fmpz(0)])
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred isfinite(P)
    @test @inferred !isinfinite(P)

# the error is/was from doing QQ(K(0)) - which is possible now
#    @test_throws MethodError E43_a1([gen(K), gen(K)])

    @test_throws ErrorException E43_a1([2, 2])

    P = @inferred infinity(E43_a1)
    @test @inferred !isfinite(P)
    @test @inferred isinfinite(P)

    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])

    P = @inferred E116_1_a1([K(0), -K(a)])
    @test typeof(P) == EllCrvPt{nf_elem}
    @test parent(P) == E116_1_a1
    @test @inferred isfinite(P)
    @test @inferred !isinfinite(P)

    P = @inferred infinity(E116_1_a1)
    @test !isfinite(P)
    @test isinfinite(P)

    @test_throws ErrorException E116_1_a1([1, 1])

    P = @inferred Eshort([2, 4], false)
    @test @inferred isfinite(P)
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == Eshort
  end

  @testset "Discriminant" begin
    @test (2*a + 10)*OK == @inferred (disc(E116_1_a1)*OK)
    @test -43 == @inferred disc(E43_a1)
  end

  @testset "j-invariant" begin
    b = (fmpq(-215055, 58) * a - fmpq(65799, 29))
    @test  b == @inferred j_invariant(E116_1_a1)
    @test fmpq(-4096, 43) == @inferred j_invariant(E43_a1)
  end

  @testset "Point aritmetic" begin
    #addition
    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])
    O = infinity(E43_a1)

    @test P == @inferred P + O
    @test E43_a1([2, -4]) == P + P
    @test E43_a1([fmpq(-2, 9), fmpq(1, 27)]) == P + P + P

    P = @inferred E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([1, -1]) == @inferred P + P
    @test E116_1_a1([K(1), -a]) == @inferred P + P + P
    @test infinity(E116_1_a1) == @inferred P + P + P + P + P
    @test P == @inferred P + infinity(E116_1_a1)

    P = @inferred Eshort([2, 4])
    O = @inferred infinity(Eshort)
    @test_throws ErrorException Eshort([-1, -1])

    @test Eshort([0, 0]) == @inferred P + P
    @test P == @inferred O + P

    P = Eshort([2, 4])
    @test Eshort([2, -4]) == @inferred -P
    P = infinity(Eshort)
    @test P == @inferred -P

    P = E43_a1([-1, 0])
    @test E43_a1([-1, -1]) == @inferred -P
    P = infinity(E43_a1)
    @test P == @inferred -P

    # inversion
    P = @inferred E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([0, 0]) == @inferred -P
    P = infinity(E116_1_a1)
    @test P == @inferred -P

    # equality
    P1 = Eshort([2, 4])
    @test @inferred ==(P1, P1)
    P2 = infinity(Eshort)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    P1 = E43_a1([-1, 0])
    @test @inferred ==(P1, P1)
    P2 = infinity(E43_a1)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    P1 = E116_1_a1([K(0), -K(a)])
    @test @inferred ==(P1, P1)
    P2 = infinity(E116_1_a1)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    # scalar multiplication
    P1 = Eshort([2, 4])
    @test Eshort([0, 0]) == @inferred 2*P1
    @test infinity(Eshort) == @inferred 4*P1

    P1 = E43_a1([-1, 0])
    @test E43_a1([fmpq(11, 49), fmpq(-363, 343)]) == @inferred 4*P1

    P1 = E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([K(0), K(0)]) == @inferred 4*P1
    @test infinity(E116_1_a1) == @inferred 5*P1
  end
end


text:
using CmeSIM: Acquisitions, FluorophoreMarkers
using Test
using Images
using Unitful: μm

@testset "CmeSIM.jl" begin
    global unreconstructedimages = load("./assets/test_unreconstructed.tiff")
    global unreconstructedimages_rgb = load("./assets/test_unreconstructed_rgb.tiff")
    global reconstructedimages = load("./assets/test_reconstructed.tiff")
    global reconstructedimages_rgb = load("./assets/test_reconstructed_rgb.tiff")
    global NA = 1.49
    global M = 109
    global pixelsize = 10μm

    include("./acquisitions_test.jl")
    include("./track_test.jl")
end


text:
module TestExtras
using Compat
using Compat.Test
using CategoricalArrays

const ≅ = isequal

@testset "cut($(Union{Int, T})[...])" for T in (Union{}, Missing)
    x = @inferred cut(Vector{Union{Int, T}}([2, 3, 5]), [1, 3, 6])
    @test x == ["[1, 3)", "[3, 6)", "[3, 6)"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[1, 3)", "[3, 6)"]

    err = @test_throws ArgumentError cut(Vector{Union{T, Int}}([2, 3, 5]), [3, 6])
    if T === Missing
        @test err.value.msg == "value 2 (at index 1) does not fall inside the breaks: adapt them manually, or pass extend=true or allow_missing=true"
    else
        @test err.value.msg == "value 2 (at index 1) does not fall inside the breaks: adapt them manually, or pass extend=true"
    end

    err = @test_throws ArgumentError cut(Vector{Union{T, Int}}([2, 3, 5]), [2, 5])
    if T === Missing
        @test err.value.msg == "value 5 (at index 3) does not fall inside the breaks: adapt them manually, or pass extend=true or allow_missing=true"
    else
        @test err.value.msg == "value 5 (at index 3) does not fall inside the breaks: adapt them manually, or pass extend=true"
    end

    if T === Missing
        x = @inferred cut(Vector{Union{T, Int}}([2, 3, 5]), [2, 5], allow_missing=true)
        @test x ≅ ["[2, 5)", "[2, 5)", missing]
        @test isa(x, CategoricalVector{Union{String, T}})
        @test isordered(x)
        @test levels(x) == ["[2, 5)"]
    else
        err = @test_throws ArgumentError cut(Vector{Int}([2, 3, 5]), [2, 5], allow_missing=true)
        @test err.value.msg == "value 5 (at index 3) does not fall inside the breaks: adapt them manually, or pass extend=true"
    end

    x = @inferred cut(Vector{Union{T, Int}}([2, 3, 5]), [3, 6], extend=true)
    @test x == ["[2, 3)", "[3, 6]", "[3, 6]"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[2, 3)", "[3, 6]"]

    x = @inferred cut(Vector{Union{T, Int}}([2, 3, 5, 6]), [3, 6], extend=true)
    @test x == ["[2, 3)", "[3, 6]", "[3, 6]", "[3, 6]"]
    @test isordered(x)
    @test levels(x) == ["[2, 3)", "[3, 6]"]

    x = @inferred cut(Vector{Union{T, Int}}([1, 2, 4]), [1, 3, 6])
    @test x == ["[1, 3)", "[1, 3)", "[3, 6)"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[1, 3)", "[3, 6)"]

    x = @inferred cut(Vector{Union{T, Int}}([1, 2, 4]), [3, 6], extend=true)
    @test x == ["[1, 3)", "[1, 3)", "[3, 6]"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[1, 3)", "[3, 6]"]

    x = @inferred cut(Vector{Union{T, Int}}([1, 2, 4]), [3], extend=true)
    @test x == ["[1, 3)", "[1, 3)", "[3, 4]"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[1, 3)", "[3, 4]"]

    x = @inferred cut(Vector{Union{T, Int}}([1, 5, 7]), [3, 6], extend=true)
    @test x == ["[1, 3)", "[3, 6)", "[6, 7]"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[1, 3)", "[3, 6)", "[6, 7]"]

    ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
    breaks = [18, 25, 35, 60, 100]
    x = @inferred cut(Vector{Union{T, Int}}(ages), breaks)
    @test x == ["[18, 25)", "[18, 25)", "[25, 35)", "[25, 35)", "[18, 25)", "[18, 25)",
                "[35, 60)", "[25, 35)", "[60, 100)", "[35, 60)", "[35, 60)", "[25, 35)"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[18, 25)", "[25, 35)", "[35, 60)", "[60, 100)"]

    breaks = [1, 6, 3] # Unsorted breaks
    labels = ["b", "a"] # Differs from lexical ordering
    x = @inferred cut(Vector{Union{T, Int}}([2, 3, 5]), breaks, labels=labels)
    @test x == ["b", "a", "a"]
    @test isa(x, CategoricalVector{Union{String, T}})
    @test isordered(x)
    @test breaks == [1, 6, 3] # Check that breaks are copied before sorting
    labels[1] = "c" # Check that labels are copied
    @test levels(x) == ["b", "a"]

    x = @inferred cut(Matrix{Union{Float64, T}}([-1.1 3.0; 1.456 10.394]), [-2.134, 3.0, 12.5])
    @test x == ["[-2.134, 3.0)" "[3.0, 12.5)"; "[-2.134, 3.0)" "[3.0, 12.5)"]
    @test isa(x, CategoricalMatrix{Union{String, T}})
    @test isordered(x)
    @test levels(x) == ["[-2.134, 3.0)", "[3.0, 12.5)"]
end

# TODO: test on arrays supporting missing values once a quantile() method is provided for them
@testset "cut([5, 4, 3, 2], 2)" begin
    x = @inferred cut([5, 4, 3, 2], 2)
    @test x == ["[3.5, 5.0]", "[3.5, 5.0]", "[2.0, 3.5)", "[2.0, 3.5)"]
    @test isa(x, CategoricalArray)
    @test isordered(x)
    @test levels(x) == ["[2.0, 3.5)", "[3.5, 5.0]"]
end

@testset "cut with formatter function" begin
  my_formatter(from, to, i; closed) = "$i: $from -- $to"

  x = 0.15:0.20:0.95
  p = [0, 0.4, 0.8, 1.0]

  @test cut(x, p, labels=my_formatter) ==
      ["1: 0.0 -- 0.4", "1: 0.0 -- 0.4", "2: 0.4 -- 0.8", "2: 0.4 -- 0.8", "3: 0.8 -- 1.0"]
end

end


text:
module TestSubDataTable
    using Base.Test
    using DataTables

    @testset "view -- DataTable" begin
        dt = DataTable(x = 1:10, y = 1.0:10.0)
        @test view(dt, 1) == head(dt, 1)
        @test view(dt, UInt(1)) == head(dt, 1)
        @test view(dt, BigInt(1)) == head(dt, 1)
        @test view(dt, 1:2) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8))) == head(dt, 2)
        @test view(dt, [1, 2]) == head(dt, 2)
        @test view(dt, 1, :x) == head(dt[[:x]], 1)
        @test view(dt, 1:2, :x) == head(dt[[:x]], 2)
        @test view(dt, vcat(trues(2), falses(8)), :x) == head(dt[[:x]], 2)
        @test view(dt, [1, 2], :x) == head(dt[[:x]], 2)
        @test view(dt, 1, 1) == head(dt[[:x]], 1)
        @test view(dt, 1:2, 1) == head(dt[[:x]], 2)
        @test view(dt, vcat(trues(2), falses(8)), 1) == head(dt[[:x]], 2)
        @test view(dt, [1, 2], 1) == head(dt[[:x]], 2)
        @test view(dt, 1, [:x, :y]) == head(dt, 1)
        @test view(dt, 1:2, [:x, :y]) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), [:x, :y]) == head(dt, 2)
        @test view(dt, [1, 2], [:x, :y]) == head(dt, 2)
        @test view(dt, 1, [1, 2]) == head(dt, 1)
        @test view(dt, 1:2, [1, 2]) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), [1, 2]) == head(dt, 2)
        @test view(dt, [1, 2], [1, 2]) == head(dt, 2)
        @test view(dt, 1, trues(2)) == head(dt, 1)
        @test view(dt, 1:2, trues(2)) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), trues(2)) == head(dt, 2)
        @test view(dt, [1, 2], trues(2)) == head(dt, 2)
        @test view(dt, Integer[1, 2]) == head(dt, 2)
        @test view(dt, UInt[1, 2]) == head(dt, 2)
        @test view(dt, BigInt[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{Int}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{Integer}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{UInt}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{BigInt}[1, 2]) == head(dt, 2)
        @test view(dt, NullableArray([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{Integer}([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{UInt}([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{BigInt}([1, 2])) == head(dt, 2)
        @test_throws NullException view(dt, [Nullable(), 1])
    end

    @testset "view -- SubDataTable" begin
        dt = view(DataTable(x = 1:10, y = 1.0:10.0), 1:10)
        @test view(dt, 1) == head(dt, 1)
        @test view(dt, UInt(1)) == head(dt, 1)
        @test view(dt, BigInt(1)) == head(dt, 1)
        @test view(dt, 1:2) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8))) == head(dt, 2)
        @test view(dt, [1, 2]) == head(dt, 2)
        @test view(dt, 1, :x) == head(dt[[:x]], 1)
        @test view(dt, 1:2, :x) == head(dt[[:x]], 2)
        @test view(dt, vcat(trues(2), falses(8)), :x) == head(dt[[:x]], 2)
        @test view(dt, [1, 2], :x) == head(dt[[:x]], 2)
        @test view(dt, 1, 1) == head(dt[[:x]], 1)
        @test view(dt, 1:2, 1) == head(dt[[:x]], 2)
        @test view(dt, vcat(trues(2), falses(8)), 1) == head(dt[[:x]], 2)
        @test view(dt, [1, 2], 1) == head(dt[[:x]], 2)
        @test view(dt, 1, [:x, :y]) == head(dt, 1)
        @test view(dt, 1:2, [:x, :y]) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), [:x, :y]) == head(dt, 2)
        @test view(dt, [1, 2], [:x, :y]) == head(dt, 2)
        @test view(dt, 1, [1, 2]) == head(dt, 1)
        @test view(dt, 1:2, [1, 2]) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), [1, 2]) == head(dt, 2)
        @test view(dt, [1, 2], [1, 2]) == head(dt, 2)
        @test view(dt, 1, trues(2)) == head(dt, 1)
        @test view(dt, 1:2, trues(2)) == head(dt, 2)
        @test view(dt, vcat(trues(2), falses(8)), trues(2)) == head(dt, 2)
        @test view(dt, [1, 2], trues(2)) == head(dt, 2)
        @test view(dt, Integer[1, 2]) == head(dt, 2)
        @test view(dt, UInt[1, 2]) == head(dt, 2)
        @test view(dt, BigInt[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{Int}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{Integer}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{UInt}[1, 2]) == head(dt, 2)
        @test view(dt, Nullable{BigInt}[1, 2]) == head(dt, 2)
        @test view(dt, NullableArray([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{Integer}([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{UInt}([1, 2])) == head(dt, 2)
        @test view(dt, NullableArray{BigInt}([1, 2])) == head(dt, 2)
        @test_throws NullException view(dt, [Nullable(), 1])
    end
end


text:
using StaticArrays, LinearAlgebra, BenchmarkTools
import TrajectoryOptimization: dynamics, jacobian, discrete_dynamics, discrete_jacobian
import TrajectoryOptimization: RigidBody

# Get Quadrotor Model
model = Dynamics.Quadrotor2{UnitQuaternion{Float64,VectorPart}}()
n,m = size(model)
N = 101  # knot points

# Build trajectory
X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
dt = fill(0.1, N)
Z = Traj(X,U,dt)  # vector of KnotPoints

# Allocate storage arrays
fVal = [@SVector zeros(n) for k = 1:N-1]
∇f = [@SMatrix zeros(n,n+m) for k = 1:N-1]
∇fd = [@SMatrix zeros(n,n+m+1) for k = 1:N-1]  # discrte Jacobian also takes Jacobian wrt time

# Continuous Dynamics
function dynamics(fVal, model::AbstractModel, Z::Traj)
    @inbounds for k in eachindex(fVal)
        # Most straightforward, passes through 2 "conversion functions":
        #  dynamics(model, z::KnotPoint) = dynamics(model, state(z), control(z), z.t)
        #  dynamics(model, x, u, t) = dynamics(model, x, u)
        #       The first extracts the state, control, and time from KnotPoint
        #       The second will drop dependence on time, by default
        fVal[k] = dynamics(model, Z[k])

        # ALTERNATE: call the dynamics function directly, by-passing the above pipeline
        #   dynamics(model, state(Z[k]), control(Z[k]))
        #
        # This function is defined in the Dynamics module
    end
end

# Continuous Jacobian
function jacobian(∇f, model::AbstractModel, Z::Traj)
    @inbounds for k in eachindex(∇f)
        # This function is defined in static_model.jl::93
        #   It's best to pass in the KnotPoint directly, since having the state and control
        #   stacked together avoids some allocations when using ForwardDiff
        ∇f[k] = jacobian(model, Z[k])
    end
end

# Discrete Dynamics
function discrete_dynamics(fVal, model::AbstractModel, Z::Traj)
    @inbounds for k in eachindex(fVal)
        # Most straightforward, passes through 1 "conversion function":
        #   discrete_dynamics(::Type{Q}, model, z::KnotPoint) = discrete_dynamics(Q, model, state(z), control(z), z.t, z.dt)
        fVal[k] = discrete_dynamics(RK3, model, Z[k])

        # ALTERNATE: bypass the extra function by making the call directly to the RK3 function
        #   defined in integration.jl::5
        #   discrete_dynamics(RK3, model, state(z), control(z), z.t, z.dt)
        end
    end
end

# Discrete Dynamics jacobian
function discrete_jacobian(∇fd, model::AbstractModel, Z::Traj)
    @inbounds for k in eachindex(∇fd)
        # This function is defined in static_model.jl::138
        ∇fd[k] = discrete_jacobian(RK3, model, Z[k])
    end
end

# Test functions
dynamics(fVal, model, Z)
model
get_orientation(model::RigidBody{UnitQuaternion{T,D}}, x::SVector{N,T2}) where {T,D,N,T2} =
    UnitQuaternion{T2,D}(x[4],x[5],x[6],x[7])
x,u = rand(model)
normalize(get_orientation(model,x))
jacobian(∇f, model, Z)
discrete_dynamics(fVal, model, Z)
discrete_jacobian(∇fd, model, Z)

# Benchmark functions (should all have 0 allocations)
@btime dynamics($fVal, $model, $Z)
@btime jacobian($∇f, $model, $Z)
@btime discrete_dynamics($fVal, $model, $Z)
@btime discrete_jacobian($∇fd, $model, $Z)


text:
function update_eps!(i::Int, s::State, c::Constants, d::Data)
  a, b = params(c.eps_prior[i])
  num0 = sum(s.lam[i] .== 0)
  a += num0
  b += (d.N[i] - num0)
  s.eps[i] = rand(Beta(a, b))
end

function update_eps!(s::State, c::Constants, d::Data)
  for i in 1:d.I
    update_eps!(i, s, c, d)
  end
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
using Test

using EquationsOfState.FiniteStrains

@testset "Test data from Pymatgen" begin
    volumes = [
        25.987454833,
        26.9045702104,
        27.8430241908,
        28.8029649591,
        29.7848370694,
        30.7887887064,
        31.814968055,
        32.8638196693,
        33.9353435494,
        35.0299842495,
        36.1477417695,
        37.2892088485,
        38.4543854865,
        39.6437162376,
        40.857201102,
        42.095136449,
        43.3579668329,
        44.6456922537,
        45.9587572656,
        47.2973100535,
        48.6614988019,
        50.0517680652,
        51.4682660281,
        52.9112890601,
        54.3808371612,
        55.8775030703,
        57.4014349722,
        58.9526328669,
    ]
    v0 = 40.98926572870838
    @test map(v -> getstrain(EulerianStrain(), v0, v), volumes) ≈ [
        0.17749734613767998,
        0.16201229839015363,
        0.14705196468522874,
        0.13259433633439055,
        0.11861428460608225,
        0.10509239492761757,
        0.09201008130906252,
        0.07934607736917354,
        0.06708557020358508,
        0.0552093893237503,
        0.043704151717372075,
        0.032551138219835574,
        0.021738453213956288,
        0.011250631487461304,
        0.0010768691828133559,
        -0.00879571243518995,
        -0.01838039465214164,
        -0.02768646727024454,
        -0.03672601609939291,
        -0.04550846145018128,
        -0.05404287277132974,
        -0.06233970600633609,
        -0.07040712280723271,
        -0.07825377702272912,
        -0.08588638996615916,
        -0.09331442415000685,
        -0.10054463011940301,
        -0.1075828678792159,
    ]
    @test [getstrain(LagrangianStrain(), v0, v) for v in volumes] ≈ [
        -0.13099486451833953,
        -0.12236351105872695,
        -0.11363226812607319,
        -0.10480202613156259,
        -0.09587095509895383,
        -0.08683995684674645,
        -0.07770989398154193,
        -0.06847899767396914,
        -0.059149424468252,
        -0.04971943052962041,
        -0.04019111457888069,
        -0.030561514081675456,
        -0.020832711371037982,
        -0.011003048988641961,
        -0.0010745548727579268,
        0.008953212194864557,
        0.0190818588446644,
        0.029309415622954416,
        0.039637468728734304,
        0.050065246093897176,
        0.06059200478213067,
        0.0712192845248587,
        0.08194633401200402,
        0.09277353626347284,
        0.10369906697722475,
        0.11472551485871496,
        0.125852144820906,
        0.13707718021610538,
    ]
    @test [getstrain(NaturalStrain(), v0, v) for v in volumes] ≈ [
        -0.15189876859201823,
        -0.14033801748544264,
        -0.1289092546069999,
        -0.1176106297063706,
        -0.10643692669611844,
        -0.09538653341846574,
        -0.08445778277754971,
        -0.07364595812672871,
        -0.06295105569265277,
        -0.05236861111679098,
        -0.04189858087249365,
        -0.03153541438302606,
        -0.02127915820568489,
        -0.011125922058677679,
        -0.0010757111979426712,
        0.008873996435816848,
        0.018726748238692272,
        0.028482534552264627,
        0.0381447307354657,
        0.04771440101503085,
        0.05719263915631435,
        0.06658253672817539,
        0.07588506772815591,
        0.0851021663162455,
        0.09423387093583208,
        0.10328387622610548,
        0.11225302719931324,
        0.1211413559225737,
    ]
    @test [getstrain(InfinitesimalStrain(), v0, v) for v in volumes] ≈ [
        -0.1640423928171002,
        -0.1506626772344306,
        -0.13758688871244362,
        -0.12480606002491879,
        -0.1123077673073063,
        -0.10008399218206754,
        -0.0881269055666829,
        -0.07642563827621052,
        -0.06497471350599215,
        -0.05376410009427657,
        -0.04278871466598844,
        -0.0320379239348092,
        -0.021507173948334346,
        -0.011188045308548222,
        -0.0010762899827498895,
        0.00883473873948748,
        0.01855249213433896,
        0.02808073099690478,
        0.03742638317829727,
        0.046593960004638224,
        0.05558787891231487,
        0.06441430751249311,
        0.07307726622682587,
        0.08158155182153393,
        0.08993010154841308,
        0.09812908257335062,
        0.10618193139700183,
        0.11409127770318894,
    ]
end


text:
using AutoDrivers.GaussianMixtureRegressionDrivers
using Distributions

include("test_GMR.jl")


text:
using Test, CSV, DelimitedFiles,
ModelSelection.CrossValidation,
ModelSelection.AllSubsetRegression,
ModelSelection.PreliminarySelection,
ModelSelection.Preprocessing

data_fat = CSV.read(DATABASE_FAT)

@testset "Kfold cross validation" begin
    @testset "Fat dataset" begin

        #funcgen("y *", data; ttest=true, cross=(k=3))

        dataorig = Preprocessing.input("y *", data_fat)

        datalassogen, vars = PreliminarySelection.lasso(dataorig)

        AllSubsetRegression.ols!(datalassogen, ttest=true)

        info = CrossValidation.kfoldcrossvalidation(datalassogen, dataorig, 3, 0.03)

        #exportar a latex, con info y best model.

        @test 1 == 1
    end
end


text:
using DecFP
using DecFP: Dec128

Base.convert(::Type{Dec128}, x::Int64) = convert(Dec128, convert(Float64, x))

include("../matplotlib_settings.jl")

using GeometricExamples

include("guiding_center_4d_settings_vprk_quad.jl")

tableau_list = get_tableau_list_vprk_projection(T=Dec128)

equ = problem.guiding_center_4d_iode_dec128()

include("guiding_center_4d.jl")

run_tableau_list(tableau_list, equ, Dec128(Δt), ntime)

generate_html(tableau_list, figures, prefix, fig_suff, fig_dir, out_dir, html_file;
              title=coll_title, params=get_params(Δt, ntime, problem))


text:
using Symplectic
a1, a2, a3 = [1 0; 1 1], [2 3; 4 5], [7 8;9 10]
b1, b2, b3, b4 = [6 3; 9 4], [1 7; 3 2], [4 3; 5 2], [6 1; 8 3]
sequence = [a1, a2, a3]
interspersion = [b1, b2, b3, b4]
ms = MatrixSequence( [a1, a2, a3] )
ms2 = MatrixSequence(sequence, interspersion)
replaceSequence(ms2, ms, 2)


text:
abstract type AbstractVarianceAccumulator{T} end

"""
    isempty(V::AbstractVarianceAccumulator)

Returns true if the given variance accumulator is empty.
"""
function Base.isempty(::AbstractVarianceAccumulator) end

"""
    empty!(V::AbstractVarianceAccumulator)

Clear the given variance accumulator.
"""
function Base.empty!(::AbstractVarianceAccumulator) end

Base.:(==)(a::AbstractVarianceAccumulator, b::AbstractVarianceAccumulator) = false
Base.isapprox(a::AbstractVarianceAccumulator, b::AbstractVarianceAccumulator; kwargs...) = false

"""
    mean(V::AbstractVarianceAccumulator)

Calculates the mean of a given variance accumulator.
"""
function mean(V::AbstractVarianceAccumulator) end

"""
    varN(V::AbstractVarianceAccumulator)

Calculates the variance/N of a given variance accumulator.
"""
varN(V::AbstractVarianceAccumulator{T}) where T = var(V) / V.count

"""
    var(V::AbstractVarianceAccumulator)

Calculates the variance of a given variance accumulator.
"""
function var(::AbstractVarianceAccumulator) end

"""
    push!(V::AbstractVarianceAccumulator{T}, value::T)

Pushes a new value into the variance accumulator.
"""
function Base.push!(::AbstractVarianceAccumulator{T}, value::T) where T end

@inline _prod(x::T, y::T) where T = x * y
@inline _prod(x::Complex, y::Complex) = Complex(real(x) * real(y), imag(x) * imag(y))
@inline _prod(x::AbstractArray, y::AbstractArray) = _prod.(x, y)

"""
    FastVariance <: AbstractVarianceAccumulator

The `FastVariance` accumulator keeps track of `∑x` and `∑x²` where `x` refers to
the pushed values. It is significantly faster than the `Variance` accumulator,
but can become unstable for some inputs. (For example when the mean is much
large than the variance of the pushed values.)

See also: [`AbstractVarianceAccumulator`](@ref), [`Variance`](@ref)
"""
mutable struct FastVariance{T} <: AbstractVarianceAccumulator{T}
    x_sum::T
    x2_sum::T
    count::Int
end
function FastVariance{T}(T0) where T
    FastVariance{T}(deepcopy(T0), deepcopy(T0), zero(Int))
end
FastVariance{T}() where T = FastVariance(zero(T), zero(T), zero(Int))
FastVariance() = FastVariance{Float64}()

@inline Base.isempty(V::FastVariance) = iszero(V.count)
function Base.empty!(V::FastVariance)
    V.x_sum, V.x2_sum = zero(V.x_sum), zero(V.x2_sum)
    V.count = zero(Int)
    return V
end
function Base.isapprox(a::FastVariance{T}, b::FastVariance{T}; kwargs...) where T
    return isapprox(a.x_sum, b.x_sum; kwargs...) &&
           isapprox(a.x2_sum, b.x2_sum; kwargs...) &&
           isapprox(a.count, b.count; kwargs...)
end
function Base.:(==)(a::FastVariance{T}, b::FastVariance{T}) where T
    return (a.x_sum == b.x_sum) && (a.x2_sum == b.x2_sum) && (a.count == b.count)
end
function Base.copy(V::FastVariance{T}) where T
    return FastVariance{T}(copy(V.x_sum), copy(V.x2_sum), copy(V.count))
end

mean(V::FastVariance) = isempty(V) ? zero(V.x_sum) : V.x_sum / V.count

function var(V::FastVariance{T}) where {T <: Real}
    n, X, X2 = V.count, V.x_sum, V.x2_sum
    return X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(V::FastVariance{T}) where {T <: Complex}
    n, X, X2 = V.count, V.x_sum, V.x2_sum
    return (real(X2) + imag(X2)) / (n - 1) - abs2(X) / (n*(n - 1))
end

function var(V::FastVariance{<: AbstractArray{T, D}}) where {D, T <: Real}
    n, X, X2 = V.count, V.x_sum, V.x2_sum
    return @. X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(V::FastVariance{<: AbstractArray{T, D}}) where {D, T <: Complex}
    n, X, X2 = V.count, V.x_sum, V.x2_sum
    return @. (real(X2) + imag(X2)) / (n - 1) - abs2(X) / (n*(n - 1))
end

function _push!(V::FastVariance{T}, value::T) where T
    V.x_sum += value
    V.x2_sum += _prod(value, value)
    V.count += 1
    return V
end

function _push!(V::FastVariance{T}, value::T) where T <: AbstractArray
    V.x_sum .+= value
    @. V.x2_sum += _prod(value, value)
    V.count += 1
    return V
end

"""
    Variance <: AbstractVarianceAccumulator

The `Variance` accumulator uses Welfords algorithm to keep track of pushed
values and to compute variances and errors. It is less error prone than
`FastVariance` but also slower.

See also: [`AbstractVarianceAccumulator`](@ref), [`FastVariance`](@ref)
"""
mutable struct Variance{T} <: AbstractVarianceAccumulator{T}
    δ::T
    m1::T
    m2::T
    count::Int
end
function Variance{T}(T0) where {T}
    Variance{T}(deepcopy(T0), deepcopy(T0), deepcopy(T0), zero(Int))
end
Variance{T}() where T = Variance{T}(zero(T), zero(T), zero(T), zero(Int))
Variance() = Variance{Float64}()

@inline Base.isempty(V::Variance) = iszero(V.count)
function Base.empty!(V::Variance)
    V.m1, V.m2 = zero(V.m1), zero(V.m2)
    V.count = zero(Int)
    return V
end

function Base.isapprox(a::Variance{T}, b::Variance{T}; kwargs...) where T
    return isapprox(a.m1, b.m1; kwargs...) &&
           isapprox(a.m2, b.m2; kwargs...) &&
           isapprox(a.count, b.count; kwargs...)
end
function Base.:(==)(a::Variance{T}, b::Variance{T}) where T
    return (a.m1 == b.m1) && (a.m2 == b.m2) && (a.count == b.count)
end
function Base.copy(V::Variance{T}) where T
    return Variance{T}(copy(V.δ), copy(V.m1), copy(V.m2), copy(V.count))
end

mean(V::Variance) = V.m1

var(V::Variance{T}) where {T <: Real} = V.m2 / (V.count - 1)
var(V::Variance{T}) where {T <: Complex} = (real(V.m2) + imag(V.m2)) / (V.count - 1)
var(V::Variance{<: AbstractArray{T, D}}) where {D, T <: Real} =
    @. V.m2 / (V.count - 1)
var(V::Variance{<: AbstractArray{T, D}}) where {D, T <: Complex} =
    @. (real(V.m2) + imag(V.m2)) / (V.count - 1)

function _push!(V::Variance{T}, value::S) where {S, T}
    δ = value - mean(V)
    V.count += 1
    V.m1 += δ / V.count
    V.m2 += _prod(δ, value - mean(V))
    return V
end

function _push!(V::Variance{T}, value::S) where {S, T <: AbstractArray}
    @. V.δ = value - V.m1
    V.count += 1
    @. V.m1 += V.δ / V.count
    @. V.m2 += _prod(V.δ, value - V.m1)
    return V
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
input = readlines("day4.txt")

draws = split(input[1], ",") .|> x -> parse(Int64,x)

# Create a 3D array with the cards
A = zeros(Int64,5, 5, Int64((size(input,1)-1)/6))
i,j=1,1
for line in input[3:end]
    if line == ""
        i+=1
        j=1
    else
        A[j,:,i] = split(line).|> x -> parse(Int64,x)
        j+=1
    end
end

# Run the simulation
i=1
stop = 0
check = 0
while stop == 0
    check = draws[i]
    A[findall(x->x==check, A)] .= -1

    keep_ind = 0
    keep_type = 0

    for card in 1:Int64((size(input,1)-1)/6)
        for ind in 1:5

            keep_ind = ind
            if all(y->y==-1, A[:,ind,card])
                keep_type = 2
                stop = 1
                keep_card = card
                break
            elseif all(y->y==-1, A[ind,:,card])
                keep_type = 1
                stop = 1
                keep_card = card
                break
            end
        end
    end
    i+=1
end

winner = A[:,:,keep_card]
winner[findall(x->x==-1, winner)] .= 0

sum(winner) * check


text:
function nlp_cvx_001_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear objective
    # - quadratic objective
    # - linear constraints forming a closed set
    # Variants
    #   010 - binding constraints
    #   011 - non-binding constraints

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @objective(model, Min, x)
    @constraint(model, x + y <= 5)
    @constraint(model, 2 * x - y <= 3)
    @constraint(model, 3 * x + 9 * y >= -10)
    @constraint(model, 10 * x - y >= -20)
    @constraint(model, -x + 2 * y <= 8)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -2.0430107680954848, tol = objective_tol)
    return check_solution(
        [x, y],
        [-2.0430107680954848, -0.4301075068564087],
        tol = primal_tol,
    )
end


text:
module InplaceOps

using LinearAlgebra

export @!

macro !(ex)
    if ex.head == :(=)
        C, ex2 = ex.args
        if ex2.head == :call
            op = ex2.args[1]
            length(ex2.args) == 3 || error("@! macro only supports 2-argument functions")
            A = ex2.args[2]
            B = ex2.args[3]
        else
            error("Invalid use of @! macro")
        end
    elseif ex.head in (:*=, :/=)
        op = Symbol(String(ex.head)[1:1])
        C = A = ex.args[1]
        B = ex.args[2]
    else
        error("@! does not support this expression")
    end
    C isa Symbol || error("Invalid use of @! macro")

    if op == :/
        if C == B
            error("@! macro cannot reuse array")
        elseif C == A
            return :($(esc(C)) = rdiv!($(esc(A)), $(esc(B))))
        else
            return :($(esc(C)) = rdiv!($(esc(C)), $(esc(A)), $(esc(B))))
        end
    elseif op == :\
        if C == A
            error("@! macro cannot reuse array")
        elseif C == B
            return :($(esc(C)) = ldiv!($(esc(A)), $(esc(B))))
        else
            return :($(esc(C)) = ldiv!($(esc(C)), $(esc(A)), $(esc(B))))
        end
    elseif op == :*
        if C == A
            if C == B
                error("@! macro cannot reuse array")
            end
            return :($(esc(C)) = rmul!($(esc(A)), $(esc(B))))
        elseif C == B
            return :($(esc(C)) = lmul!($(esc(A)), $(esc(B))))
        else
            return :($(esc(C)) = mul!($(esc(C)), $(esc(A)), $(esc(B))))
        end
    end
end

end  # module


text:
# Cardinality constrained least-squares example (nonconvex)
include("regressor_selection.jl")

function gen_sparse_reg_problem(m, n, p = 100/n)
    # generate sparse solution vector
    x = sprandn(n,1,p)

    # generate random data matrix
    A = randn(m,n)

    # normalize columns of A
    A = A ./ sqrt(sum(A.^2,1))

    # generate measurement b with noise
    b = A*x + sqrt(0.001)*randn(m,1)

    return A, x, b
end

# Ensure compiled
A, xtrue, b = gen_sparse_reg_problem(15, 50, 0.05);
xhat = regressor_selection_admm(A, b, iceil(50*0.05), 1.0; MaxIterations = 10, Quiet = true);

include("eval_sparse_reg.jl")

srand(2)
m = 1500       # number of examples
n = 5000       # number of features
p = 100/n      # sparsity density

A, xtrue, b = gen_sparse_reg_problem(m, n, p);
supportlen = length(find(xtrue))
#A, xtrue, b = map(fn -> readcsv(fn * ".csv"), ["A", "x", "b"]);
@time xhat = regressor_selection_admm(A, b, p*n, 1.0; MaxIterations = 1000, Quiet = false);
eval_sparse_reg(A, b, xtrue, xhat);
@time xhat2 = regressor_selection_admm(A, b, supportlen, 1.0; MaxIterations = 1000, Quiet = false);
eval_sparse_reg(A, b, xtrue, xhat2);

slens = collect(10 * (1:iceil(2*supportlen/10)));
numslens = length(slens)
res = zeros(5,numslens)
for i in 1:numslens
    xhat = regressor_selection_admm(A, b, slens[i], 1.1; MaxIterations = 100);
    res[:,i] = map(float, eval_sparse_reg(A, b, xtrue, xhat))
end


text:
using Base.Test

using CompScienceMeshes
using BEAST

const e0 = Point(0.0,0.0,0.0)
const e1 = Point(1.0,0.0,0.0)
const e2 = Point(0.0,1.0,0.0)
const e3 = Point(0.0,0.0,1.0)

p = patch(
    [
        Point(0.0,0.0,0.0),
        Point(1.0,0.0,0.0)
    ],
    Val{1}
)

q = patch(
    [
        Point(0.0,0.0,0.0),
        Point(0.5,0.0,0.0)
    ],
    Val{1}
)

f = LagrangeRefSpace{Float64,1}()
x = neighborhood(p, [0.0])
v = f(x)
@test v == Vec(0.0, 1.0)

Q = restrict(f, p, q)
@test Q == [
    1.0 0.5
    0.0 0.5
]

# Test restriction of RT elements
ni, no = 6, 7;
ui = transpose([triangleGaussA[ni] triangleGaussB[ni] ]);
uo = transpose([triangleGaussA[no] triangleGaussB[no] ]);
wi = triangleGaussW[ni];
wo = triangleGaussW[no];
universe = Universe(1.0, ui, wi, uo, wo);

p = patch([2*e0,2*e1,2*e2], Val{2})
x = neighborhood(p,[0.5, 0.5])
ϕ = RTRefSpace{Float64}()
v = ϕ(x)

Q = restrict(ϕ, p, p)
@test Q == eye(3)

q = patch([e0+e1, 2*e1, e1+e2], Val{2})
Q = restrict(ϕ, p, q)
@test Q == [
    2 -1 0
    0  1 0
    0 -1 2] // 4


text:
ops = open("boot_code.txt") do f
    map(m -> (m[1], parse(Int, m[2])), (match.(r"(\w{3}) ([+-]\d+)", readlines(f))))
end

function mutate(ops)
    mutations = []

    for (i, op) in enumerate(ops)
        ins, val = op
        if ins == "nop"
            mut = copy(ops)
            mut[i] = ("jmp", val)
            push!(mutations, mut)
        elseif ins == "jmp"
            mut = copy(ops)
            mut[i] = ("nop", val)
            push!(mutations, mut)
        end
    end

    return mutations
end

function run(ops)
    acc = 0
    i = 1
    visited = []

    while length(visited) == length(unique(visited))
        op, val = ops[i]

    #    println(length(visited)+1, " ", op, " ", val)
        if op == "acc"
            acc += val
            i += 1
        elseif op == "jmp"
            i += val
        elseif op == "nop"
            i += 1
        end

        if i > length(ops)
            break
        end

        push!(visited, i)
    end

    if i >= length(ops)
        eof = true
    else
        eof = false
    end

    return (acc, eof)
end

println(run(ops)[1])

for m in mutate(ops)
    acc, eof = run(m)
    if eof
        println(acc)
        break
    end
end


text:
using BinaryBuilder, Pkg

# Read in input `.json` file
json = String(read(ARGS[1]))
buff = IOBuffer(strip(json))
objs = []
while !eof(buff)
    push!(objs, BinaryBuilder.JSON.parse(buff))
end
# Merge the multiple outputs into one
merged = BinaryBuilder.merge_json_objects(objs)
BinaryBuilder.cleanup_merged_object!(merged)
# Determine build version
name = merged["name"]
jll_name = string(name, "_jll")
repo_name = string("JuliaBinaryWrappers/", jll_name, ".jl")

@info "Initializing $(repo_name)..."

BinaryBuilder.init_jll_package(name, joinpath(Pkg.devdir(), jll_name), repo_name)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil