text:
include("usingDG.jl")
include("EulerEqs.jl")
# Example 1 , sine initial condition
function bdr(xᵧ,yᵧ,time_now,nx,ny,ρ⁻,m⁻,n⁻,E⁻)
    ρ⁺=ρ⁻
    m⁺=m⁻
    n⁺=n⁻
    E⁺=E⁻
    return ρ⁺,m⁺,n⁺,E⁺
end
γ=1.4
ρ₀(x::Real,y::Real)=1+0.5sin(x+y); ρ₀(x::AbstractArray,y::AbstractArray)=ρ₀.(x,y);
m₀(x::Real,y::Real)=1+0.5sin(x+y); m₀(x::AbstractArray,y::AbstractArray)=m₀.(x,y);
n₀(x::Real,y::Real)=1+0.5sin(x+y); n₀(x::AbstractArray,y::AbstractArray)=n₀.(x,y);
E₀(x::Real,y::Real)=(1+0.5sin(x+y))+1/(γ-1); E₀(x::AbstractArray,y::AbstractArray)=E₀.(x,y);
ρ(x,y,t)=ρ₀(x-t,y-t)
CFL=0.1
P_order=2
T=1.0
α=1.1
xmin=0.0
ymin=0.0
xmax=2π
ymax=2π

Cell_N=60
Cell_M=60
mesh=Grid2DRect(xmin,xmax,ymin,ymax,Cell_N,Cell_M;Periodic=true)
mesh,Coef_ρ,Coef_m,Coef_n,Coef_E,time_now,num_iter,basisP = main(mesh,ρ₀,m₀,n₀,E₀,γ;
            CFL=CFL,P_order=P_order,T=T,α=α,bdr=bdr)

L2err=errLp(ρ,mesh,Coef_ρ,time_now;p=2)

# Example 2 , isentropic vortex advection problem
γ=1.4
ρ₀(x::Real,y::Real)=TT(x,y)^(1/(γ-1)); ρ₀(x::AbstractArray,y::AbstractArray)=ρ₀.(x,y);
m₀(x::Real,y::Real)=ρ₀(x,y)*(1-(y-5)*phi(x,y)); m₀(x::AbstractArray,y::AbstractArray)=m₀.(x,y);
n₀(x::Real,y::Real)=ρ₀(x,y)*(1+(x-5)*phi(x,y)); n₀(x::AbstractArray,y::AbstractArray)=n₀.(x,y);
E₀(x::Real,y::Real)=1/2*m₀(x,y)^2/ρ₀(x,y)+1/2*n₀(x,y)^2/ρ₀(x,y)+TT(x,y)^(γ/(γ-1))/(γ-1); E₀(x::AbstractArray,y::AbstractArray)=E₀.(x,y);
TT(x,y)=1-(γ-1)/(2γ)*phi(x,y)^2
phi(x,y)=5/(2π)*exp(0.5*(1-(x-5)^2-(y-5)^2))
ρ(x,y,t)=ρ₀(x-t,y-t)
function bdr(xᵧ,yᵧ,time_now,nx,ny,ρ⁻,m⁻,n⁻,E⁻)
    ρ⁺=ρ₀(xᵧ-time_now,yᵧ-time_now)
    m⁺=m₀(xᵧ-time_now,yᵧ-time_now)
    n⁺=n₀(xᵧ-time_now,yᵧ-time_now)
    E⁺=E₀(xᵧ-time_now,yᵧ-time_now)
    return ρ⁺,m⁺,n⁺,E⁺
end
CFL=0.1
P_order=2
T=1.0
α=1.1
xmin=0.0
ymin=0.0
xmax=10.0
ymax=10.0

Cell_N=60
Cell_M=60
mesh=Grid2DRect(xmin,xmax,ymin,ymax,Cell_N,Cell_M;Periodic=false)
mesh,Coef_ρ,Coef_m,Coef_n,Coef_E,time_now,num_iter,basisP = main(mesh,ρ₀,m₀,n₀,E₀,γ;
            CFL=CFL,P_order=P_order,T=T,α=α,bdr=bdr)

L2err=errLp(ρ,mesh,Coef_ρ,time_now;p=2)

# Example 3 , Burgers' like initial condition
function bdr(xᵧ,yᵧ,time_now,nx,ny,ρ⁻,m⁻,n⁻,E⁻)
    ρ⁺=ρ⁻
    m⁺=m⁻
    n⁺=n⁻
    E⁺=E⁻
    return ρ⁺,m⁺,n⁺,E⁺
end
γ=3
ρ₀(x::Real,y::Real)=(1+0.2sin((x+y)/2))/sqrt(6); ρ₀(x::AbstractArray,y::AbstractArray)=ρ₀.(x,y);
m₀(x::Real,y::Real)=ρ₀(x,y)*sqrt(γ/2)*ρ₀(x,y); m₀(x::AbstractArray,y::AbstractArray)=m₀.(x,y);
n₀(x::Real,y::Real)=ρ₀(x,y)*sqrt(γ/2)*ρ₀(x,y); n₀(x::AbstractArray,y::AbstractArray)=n₀.(x,y);
E₀(x::Real,y::Real)=1/2*m₀(x,y)^2/ρ₀(x,y)+1/2*n₀(x,y)^2/ρ₀(x,y)+ρ₀(x,y)^γ/(γ-1); E₀(x::AbstractArray,y::AbstractArray)=E₀.(x,y);
ρ(x,y,t)=Burgers(x/2,y/2,t/2;α=0.2,β=1.0)/sqrt(6)
CFL=0.1
P_order=2
T=0.5
α=1.1
xmin=0.0
ymin=0.0
xmax=4π
ymax=4π

Cell_N=60
Cell_M=60
mesh=Grid2DRect(xmin,xmax,ymin,ymax,Cell_N,Cell_M;Periodic=true)
mesh,Coef_ρ,Coef_m,Coef_n,Coef_E,time_now,num_iter,basisP = main(mesh,ρ₀,m₀,n₀,E₀,γ;
            CFL=CFL,P_order=P_order,T=T,α=α,bdr=bdr)

L2err=errLp(ρ,mesh,Coef_ρ,time_now;p=2)

# Exact Riemann solver
x=collect(range(-1,1,length=500))
t=0.15
left_rho=1.
left_u=-2.
left_p=0.4
left_gamma=1.4
right_rho=1.
right_u=2.
right_p=0.4
right_gamma=1.4
profiles=RiemannSolver(x, t,left_rho,left_u,left_p,left_gamma,right_rho,right_u,right_p,right_gamma)

density = [status.rho for status in profiles]
velocity = [status.u for status in profiles]
pressure = [status.p for status in profiles]
plot(x,density)
#=
file = matopen("Riemann.mat", "w")
write(file, "density", density)
write(file, "velocity", velocity)
write(file, "pressure", pressure)
write(file, "x", x)
close(file)
=#

#=
#  2D Double Rarefaction Riemann Problem
function bdr(x,y,t,nx,ny,ρ⁻,m⁻,n⁻,E⁻)
    if nx>0.5 || nx<-0.5
        if x<0 # left
            m⁺=-2.0
        else # right
            m⁺=2.0
        end
        ρ⁺=1.0
        n⁺=0.0
        E⁺=3.0
    else
        ρ⁺=ρ⁻
        m⁺=m⁻
        n⁺=-n⁻
        E⁺=E⁻
    end
    return ρ⁺,m⁺,n⁺,E⁺
end
γ=1.4
ρ₀(x,y)=1.0
m₀(x,y)=(x<0 ? -2.0 : 2.0)
n₀(x,y)=0.0
E₀(x,y)=3.0
ρ₀(x::AbstractArray,y::AbstractArray)=ρ₀.(x,y)
m₀(x::AbstractArray,y::AbstractArray)=m₀.(x,y)
n₀(x::AbstractArray,y::AbstractArray)=n₀.(x,y)
E₀(x::AbstractArray,y::AbstractArray)=E₀.(x,y)
xmin = -1.
xmax = 1.
ymin = -1.
ymax = 1.
Cell_N=80
Cell_M=80
mesh=Grid2DRect(xmin,xmax,ymin,ymax,Cell_N,Cell_M;Periodic=false)

mesh,Coef_ρ,Coef_m,Coef_n,Coef_E,time_now,num_iter,basisP = main(mesh,ρ₀,m₀,n₀,E₀,γ;
            CFL=0.01,P_order=2,T=0.15,α=1.1,bdr=bdr)

=#


text:
using VisualParticipationAnalytics
using Clustering
using Languages
using SparseArrays
using SQLite
using Tables
using Test
using TextAnalysis

@testset "VisualParticipationAnalytics" begin
    include("io.jl")
    include("preprocessing.jl")

    @testset "similarities" begin
        doc1 = StringDocument("a b")
        doc2 = StringDocument("b c")

        crps = Corpus([doc1, doc1])
        sims = similarities(crps)
        @test sims[1] ≈ sims[4] ≈ 1.0 && isnan(sims[2]) && isnan(sims[3])
        @test similarities(crps, false) ≈ [1.0 1.0; 1.0 1.0]

        crps = Corpus([doc1, doc2])
        @test similarities(crps) ≈ [1.0 0.0; 0.0 1.0]
        @test similarities(crps, false) ≈ [1.0 0.5; 0.5 1.0]
    end

    @testset "clustering" begin
        doc1 = StringDocument("a")
        doc2 = StringDocument("b")
        doc3 = StringDocument("c")
        crps = Corpus([doc1, doc2, doc3])

        c = clustering(crps, 2)
        @test typeof(c) <: ClusteringResult && nclusters(c) == 2

        c = clustering(crps, 2, false)
        @test typeof(c) <: ClusteringResult && nclusters(c) == 2
    end

    @testset "topicmodel" begin
        doc1 = StringDocument("a")
        doc2 = StringDocument("b")
        crps = Corpus([doc1, doc2])
        topicword, topicdoc = topicmodel(crps, 2, 100, 0.1, 0.1)
        @test typeof(topicword) == SparseMatrixCSC{Float64,Int}
        @test typeof(topicdoc) == Array{Float64,2}
    end

    @testset "topkwords" begin
        doc1 = StringDocument("programming in julia")
        doc2 = StringDocument("python programming")
        doc3 = StringDocument("julia and ada")
        crps = Corpus([doc1, doc2, doc3])
        topicword, topicdoc = topicmodel(crps, 2, 100, 0.1, 0.1)
        words = topkwords(topicword, 1, crps, 2)
        @test typeof(words) == Array{Tuple{String,Float64},1}
        @test length(words) == 2
        @test words[1][2] >= words[2][2]
    end
end

ext: jl


text:
# ---
# title: Density plots
# cover: assets/density_plots.png
# description: Visualizing kernel density estimation of data
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
using AlgebraOfGraphics: density
set_aog_theme!() #src

df = (x=randn(1000), c=rand(["a", "b"], 1000))
plt = data(df) * mapping(:x, color=:c) * density(bandwidth=0.5)
fg = draw(plt)

#

df = (x=randn(1000), c=rand(["a", "b"], 1000))
plt = data(df) * mapping(:x, color=:c) * density(bandwidth=0.5) * visual(orientation=:vertical)
"Not yet supported" # hide

# save cover image #src
mkpath("assets") #src
save("assets/density_plots.png", fg) #src


text:
using ReverseDiff: InstructionTape, GradientConfig, JacobianConfig, HessianConfig,
                   value, deriv, tape, valtype,
                   derivtype, track, track!

using Random

const COMPILED_TAPE_LIMIT = 5000

# These functions correctly emit NaNs for certain arguments, but ReverseDiff's test
# machinery is currently too dumb to handle them properly.
const SKIPPED_BINARY_SCALAR_TESTS = Symbol[:hankelh1, :hankelh1x, :hankelh2, :hankelh2x,
                                           :pow, :besselj, :besseli, :bessely, :besselk,
                                           :polygamma, :ldexp]

# make RNG deterministic, and thus make result inaccuracies
# deterministic so we don't have to retune EPS for arbitrary inputs
Random.seed!(1)

test_println(kind, f, pad = "  ") = println(pad, "testing $(kind): `$(f)`...")

@inline test_approx(A, B, _atol = 1e-5) = @test isapprox(A, B, atol = _atol)

tracked_is(a, b) = value(a) === value(b) && deriv(a) === deriv(b) && tape(a) === tape(b)
tracked_is(a::AbstractArray, b::AbstractArray) = all(map(tracked_is, a, b))
tracked_is(a::Tuple, b::Tuple) = all(map(tracked_is, a, b))

# ensure that input is in domain of function
# here `x` is a scalar or array generated with `rand(dims...)`, i.e., values of `x`
# are between 0 and 1
function modify_input(f, x)
    return if in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth))
        x .+ one(eltype(x))
    elseif f === :log1mexp || f === :log2mexp
        x .- one(eltype(x))
    else
        x
    end
end

max_stars_repo_name: jrevels/ReverseDiffPrototype

max_issues_repo_name: jrevels/ReverseDiffPrototype

max_forks_repo_name: jrevels/ReverseDiffPrototype


text:
using Parameters
using Random
using BenchmarkTools

@with_kw mutable struct Config
  char = (skip='?',less='>',more='<',num='$',klass='!')
  str  = (skip="?",)
  some = (max=32,)
  div  = (few=1024,divs=16, cohen=0.3, trivial=1.05)
  seed = 1
end

the=Config()
Random.seed!(the.seed)

same(s) = s
int(x)  = floor(Int,x)
any(a)  = a[ int(length(a) * rand()) + 1]
few(a,n=the.divs.few) = length(a) < n  ? a : [any(a)  for _ in 1:n]

function say(i)
  s,pre="$(typeof(i)){",""
  for f in sort!([x for x in fieldnames(typeof(i)) if !("$x"[1] == '_')])
    g = getfield(i,f)
    s = s * pre * "$f=$g"
    pre=", "
  end
  print(s * "}")
end

@with_kw mutable struct Num
  pos=0; txt=""; w=1; n=0; lo=10^32; hi=-1*10^32; mu=0; m2=0; sd=nothing
end
@with_kw mutable struct Some
  pos=0; txt=""; w=1; n=0; _all=[]; max=the.some.max; stale=false
end
@with_kw mutable struct Sym
  os=0; txt=""; w=1; n=0; seen=Dict();  mode=nothing; ent=nothing
end

incs!(i,inits)      = begin [inc!(i,x) for x in inits]; i end
nump(s,c=the.char)  = c.less in s || c.more in s || c.num in s
goalp(s,c=the.char) = c.less in s || c.more in s || c.klass in s

function col(txt,pos)
   x = nump(txt) ? Num : Sym
   x(txt=txt, pos=pos, w= the.char.less in txt ? -1 : 1)
end

function all(i::Some)
  if i.stale
    i._all = sort(i._all)
    i.stale=false
  end
  return i._all
end

function inc!(i,x)
  if x != the.char.skip
    i.n += 1
    inc1!(i,x)
  end
  x
end

function inc1!(i::Some, x)
  m = length(i._all)
  if m < i.max
    i.stale=true
    push!(i._all,x)
  elseif rand() < m/i.n
    i.stale=true
    i._all[ int(m*rand()) + 1 ] = x end
end

function var(i::Num)
  if     i.m2 < 0  0
  elseif i.n  < 2  0
  else             (i.m2 / (i.n - 1 + 10^-32))^0.5
  end
end

function inc1!(i::Num,x)
  i.lo  = min(i.lo, x)
  i.hi  = max(i.hi, x)
  d     = x - i.mu
  i.mu += d / i.n
  i.m2 += d * (x - i.mu)
end

function div(lst, x, y)
  lst       = few(lst,the.div.few)
  lst       = sort([z for z in lst if x(z) != the.str.skip], by=x)
  sd(a,f=x) = var(incs!(Num(), [f(z) for z in a]))
  mid(a,f=x) = f( a[ int(length(a)/2) ] )
  eps = the.div.cohen * sd(lst,x)
  function chop(a)
    m = length(a)
    tmp, out, n = [], [], m / the.div.divs
    last=nothing
    for (i,one) in enumerate(a)
      if length(tmp)>=n && m - i > n && x(one) != x(last)
        if x(one) - x(tmp[1]) > eps
          push!(out, tmp)
          tmp = []
      end end
      push!(tmp,one)
      last = one
    end
    if length(tmp) > 0 push!(out,tmp) end
    out
  end
  function merge(a)
    tmp, out, j, m = [], [], 1, length(a)
    while j <= m
      one = a[j]
      if j < m
        two       = a[j+1]
        three     = [ one ; two ]
        n1, n2, n3= length(one), length(two), length(three)
        sd1,sd2,sd3= sd(one,y), sd(two,y), sd(three,y)
        sd12      = n1/n3*sd1 + n2/n3*sd2
        if abs(sd1 - sd2) < 0.01 || sd12*the.div.trivial >sd3
          one = three
          j += 1
        end
      end
      push!(tmp,one)
      j += 1
    end
    return length(tmp) < length(a) ? merge(tmp) : a
  end
  merge( chop(lst))
end

function main()
  two(r) = if r<0.2 0 elseif r<0.4 1 elseif r<0.6 0 elseif r<0.8 1 else 0 end
  one(r) = [int(100*r), two(r)]
  lst = sort([rand() for _ in 1:10^8])
  lst = [one(r) for r in lst]
  first(z)  = z[1]
  second(z) = z[2]
  println(100)
  for (i,one) in enumerate(div(lst,first,second))
    println(i," ",length(one)," ", one[1]," ",last(one))
  end
end

main1() = for _ in 1:10; print(any([1,2,3])) end

try
  println("\n--| ",int(time() % 1000)," |----------------")
  main()
catch e
   showerror(stdout, e, catch_backtrace()[1:10])
end

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
mutable struct Fact
    uin::Vector{Float64}
    uout::Vector{Ptr{Float64}}
    neigs::Vector{Int}
    w::Vector{Float64}
    b::Int
end

Fact() = Fact(Float64[], Ptr{Float64}[], Int[], Float64[], 1)

deg(f::Fact) = length(f.uin)

mutable struct FGBMatching
    N::Int
    γ::Float64
    fnodes::Vector{Fact}
    adjlist::Vector{Vector{Int}}
end

show(io, ::FGBMatching) = show(io, "FGBMatching(...)")

function FGBMatching(net::Network; γ=Inf)
    @assert has_eprop(net, "w")
    @assert has_vprop(net, "b")

    N = nv(net)
    fnodes = [Fact() for i=1:N]
    wmap = eprop(net, "w")
    bmap = vprop(net, "b")

    # Prune Graph
    w = [Float64[] for i=1:N]
    adjlist = [Int[] for i=1:N]
    for i=1:N
        for e in edges(net, i)
            wij = wmap[e]
            if wij < γ
                push!(w[i], wij)
                push!(adjlist[i], dst(e))
            end
        end
    end

    for (i, f) in enumerate(fnodes)
        resize!(f.uin, length(adjlist[i]))
        resize!(f.uout, length(adjlist[i]))
        resize!(f.w, length(adjlist[i]))
    end

    for (i, f) in enumerate(fnodes)
        f.b = bmap[i]
        for (ki, j) in enumerate(adjlist[i])
            f.uin[ki] = 0
            kj = findfirst(==(i), adjlist[j])
            fnodes[j].uout[kj] = getref(f.uin, ki)
            f.w[ki] = w[i][ki]
        end
    end

    FGBMatching(N, γ, fnodes, adjlist)
end

function initrand!(g::FGBMatching)
    for f in g.fnodes
        for k=1:deg(f)
            f.uin[k] = randn()
        end
    end
end

"""
    topk(x, k; rev=false)

Return the indexes of the largest `k` elements in `x`.
If `rev=true`, return the smallest elements.
"""
function topk(x, k; rev=false)
    is = partialsortperm(x, 1:k; rev=!rev)
    return is # xs[is]
end

function update!(f::Fact)
    @unpack w, uin, uout, b = f
    h = w .- uin
    is = topk(h, b+1, rev=true)
    m1 = h[is[b]]
    m2 = h[is[b+1]]

    Δ = abs(uout[is[b]][] - m2)
    for i=1:deg(f)
        uout[i][] = m1
    end
    for i in is[1:b]
        uout[i][] = m2
    end
    return Δ
end

findmatch(f::Fact) = findmatch(f, f.b)

function findmatch(f::Fact, k)
    h = f.w .- f.uin
    is = topk(h, k, rev=true)
    return is
end

function oneBPiter!(g::FGBMatching)
    Δ = 0.
    for a in randperm(g.N)
        d = update!(g.fnodes[a])
        Δ = max(Δ, d)
    end
    return Δ
end

function converge!(g::FGBMatching; maxiters=100, ϵ=1e-8, verbose=true)
    Eold = 0.
    tstop = 0
    it = 0
    while it < maxiters
        it += 1
        Δ = oneBPiter!(g)
        E, matchmap, nviolations = energy(g)
        verbose && @printf("it=%d  E=%.5f  nviolations=%d \tΔ=%f \n", it, E, nviolations, Δ)

        if abs(Eold - E) < ϵ && nviolations == 0
            tstop += 1
            if tstop == 10
                verbose && println("Found ground state")
                break
            end
        else
            tstop = 0
        end

        Eold = E
    end
    return Eold, it
end

function energy(g::FGBMatching)
    @unpack fnodes, N, adjlist = g
    E = 0.
    matchmap = Vector{Vector{Int}}(undef, N)
    for i=1:N
        f = fnodes[i]
        ks = findmatch(f)
        E += sum(f.w[ks])
        matchmap[i] = adjlist[i][ks]
    end
    E /= 2
    nviolations = 0
    for i=1:N
        for j in matchmap[i]
            nviolations += i ∉ matchmap[j]
        end
    end
    return E, matchmap, nviolations
end

"""
    run_bp(net::Network; [γ, maxiters, ϵ, seed, verbose])

Computes the minimum weight perfect b-matching on graph `net`.
The input `net` has to contain the edge property `w` and
the vertex property `b`.

Can optionally impose a cutoff `γ` on the edge costs `w`
in order to sparsify the graph and improve performance.

The output object contains the following fields:

- `energy`: the optimal cost.
- `match`: where `match[i]` is the list of neighbors of vertex
           `i` in the optimal matching.
- `iters`: number of iterations performed.
- `ok`: whether the algorithm was successful or not.
- `num_violations`: the number of inconsistencies in the solution.
- `bpgraph`: a type storing the BP messages.
"""
function run_bp(net::Network;
                γ = Inf,
                maxiters = 10000,
                ϵ = 1e-4,
                seed = -1,
                verbose=true)
    seed > 0 && Random.seed!(seed)
    g = FGBMatching(net; γ)
    initrand!(g)
    E, iters = converge!(g; maxiters, ϵ, verbose)
    E, matchmap, nviolations = energy(g)

    res = (energy = E,
           match = matchmap,
           bpgraph = g,
           num_violations =  nviolations,
           iters = iters,
           ok = (nviolations == 0),)

    return res
end


text:
# This is just codegen for making some convenient names for commonly used fonts.
# don't remmember if this is in use or not anymore...

cols = "black, red, green, yellow, blue, magenta, cyan, light_gray, default, dark_gray, light_red, light_green, light_yellow, light_blue, light_magenta, light_cyan, white"
fxs = "reset, bold, faint, italics, underline, blink, negative, conceal, strikethrough"

for (col, fx) in Iterators.product(split(cols, ", "), split(fxs, ", "))
    name = Symbol(col * "_" * fx)
    global colsym = Symbol(col)
    global fxsym = Symbol(fx)
    @eval $colsym = colsym
    @eval $fxsym = fxsym
    @eval $name() = Crayon(foreground=$colsym, $fxsym=true)
end


text:
P = [1	2	3	1e3	1e3	1e3	1e3;
3	1e3	2	1	1e3	1e3	1e3;
1e3	3	1e3	1e3	1	2	1e3;
1	1e3	1e3	1e3	3	1e3	2;
1	3	1e3	1e3	2	1e3	1e3;
2	1e3	3	1e3	1	1e3	1e3;
1e3	1e3	2	1	1e3	1e3	3;
1	3	1e3	1e3	2	1e3	1e3;
1	3	1e3	1e3	2	1e3	1e3;
1e3	1e3	2	3	1e3	1e3	1;
1	2	3	1e3	1e3	1e3	1e3;
1	2	3	1e3	1e3	1e3	1e3;
1e3	1e3	1	2	1e3	1e3	3;
1e3	1e3	3	2	1e3	1e3	1e3;
2	1	1e3	3	1e3	1e3	1e3;
1	2	1e3	3	1e3	1e3	1e3;
1e3	1	1e3	2	1e3	1e3	3;
1e3	3	1	1e3	1e3	1e3	2;
1e3	1e3	1e3	3	1	2	1e3;
1e3	2	1e3	3	1e3	1e3	1;
3	2	1e3	1	1e3	1e3	1e3;
3	2	1	1e3	1e3	1e3	1e3;
1	2	3	1e3	1e3	1e3	1e3;
1e3	1e3	2	3	1e3	1e3	1;
1e3	1	3	2	1e3	1e3	1e3;
1	1e3	2	1e3	3	1e3	1e3;
1	2	1e3	1e3	1e3	1e3	3;
1	2	1e3	3	1e3	1e3	1e3;
1	2	1e3	3	1e3	1e3	1e3;
2	1	3	1e3	1e3	1e3	1e3;
1e3	1e3	3	1e3	1e3	1	2;
3	2	1e3	1e3	1	1e3	1e3;
1	1e3	2	1e3	3	1e3	1e3;
1e3	2	1e3	1e3	1	1e3	3;
3	1e3	2	1	1e3	1e3	1e3;
1e3	1	3	1e3	1e3	1e3	2;
1e3	1e3	2	1	1e3	3	1e3;
1e3	3	1	2	1e3	1e3	1e3;
1	1e3	2	3	1e3	1e3	1e3;
1	1e3	3	2	1e3	1e3	1e3;
3	1	1e3	1e3	1e3	1e3	2;
1	1e3	2	1e3	3	1e3	1e3;
1	1e3	1e3	1e3	2	3	1e3;
3	1e3	2	1e3	1e3	1e3	1;
2	1	3	1e3	1e3	1e3	1e3;
1e3	3	2	1e3	1e3	1	1e3;
1e3	2	1e3	3	1	1e3	1e3;
1	3	1e3	1e3	1e3	1e3	2;
2	3	1e3	1	1e3	1e3	1e3;
1e3	1	1e3	2	1e3	1e3	3;]


text:
@doc raw"""
    ButcherTableau <: AbstractButcherTableau

A composite type for the Butcher tableau of an [`AbstractRungeKuttaSolver`](@ref).

# Constructors
```julia
ButcherTableau(A, b, c, s, p[, d, q])
ButcherTableau(tableau)
```

# Arguments
- `A :: AbstractMatrix`       : matrix of coefficients.
- `b :: AbstractVector`       : vector of weights.
- `c :: AbstractVector`       : vector of nodes.
- `s :: Integer`              : number of stages.
- `p :: Integer`              : order of accuracy.
- `d :: AbstractVector`       : embedding's vector of weights.
- `q :: Integer`              : embedding's order of accuracy.
- `tableau :: AbstractMatrix` : matrix from which all other fields are automatically constructed, with structure:
```math
\begin{array}{c|c}
    c & A \\
    \hline
    p & b \\
    q & d
\end{array}
```

# Functions
- [`show`](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
struct ButcherTableau{A_T, b_T, c_T, s_T, p_T, d_T, q_T} <: AbstractButcherTableau
    A::A_T
    b::b_T
    c::c_T
    s::s_T
    p::p_T
    d::d_T
    q::q_T
end

function ButcherTableau(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, s::Integer, p::Integer)
    d = nothing
    q = nothing
    return ButcherTableau(A, b, c, s, p, d, q)
end

function ButcherTableau(tableau::AbstractMatrix)
    nrows, ncols = size(tableau)
    s = ncols - 1
    A = tableau[1:s, 2:ncols]
    b = tableau[s+1, 2:ncols]
    c = tableau[1:s, 1]
    p = convert(Integer, tableau[s+1, 1])
    if nrows == ncols
        return ButcherTableau(A, b, c, s, p)
    else
        d = tableau[nrows, 2:ncols]
        q = convert(Integer, tableau[nrows, 1])
        return ButcherTableau(A, b, c, s, p, d, q)
    end
end

############################################################################################
#                                         PRINTING                                         #
############################################################################################

"""
    show(io::IO, tableau::AbstractButcherTableau)

prints a full description of `tableau` and its contents to a stream `io`.
"""
Base.show(io::IO, tableau::AbstractButcherTableau) = NSDEBase._show(io, tableau)

"""
    summary(io::IO, tableau::AbstractButcherTableau)

prints a brief description of `tableau` to a stream `io`.
"""
Base.summary(io::IO, tableau::AbstractButcherTableau) = NSDEBase._summary(io, tableau)


text:
module JuliennedArrays

import Base: axes, getindex, setindex!, size
using Base: @pure, tail

@inline is_in(needle::Needle, straw1::Needle, straws...) where {Needle} = True()
@inline is_in(needle, straw1, straws...) = is_in(needle, straws...)
@inline is_in(needle) = False()

@inline in_unrolled(straws, needle1, needles...) =
    is_in(needle1, straws...), in_unrolled(straws, needles...)...
@inline in_unrolled(straws) = ()

@pure as_vals(them::Int...) = map(Val, them)

abstract type TypedBool end
"""
    struct True

Typed `true`
"""
struct True <: TypedBool end
"""
    struct False

Typed `false`
"""
struct False <: TypedBool end

@inline untyped(::True) = true
@inline untyped(::False) = false

@inline not(::False) = True()
@inline not(::True) = False()

export True
export False

@inline getindex_unrolled(into::Tuple{}, switches::Tuple{}) = ()
@inline function getindex_unrolled(into, switches)
    next = getindex_unrolled(tail(into), tail(switches))
    if untyped(first(switches))
        (first(into), next...)
    else
        next
    end
end

@inline setindex_unrolled(old::Tuple{}, something, ::Tuple{}) = ()
@inline setindex_unrolled(old, new, switches) =
    if untyped(first(switches))
        first(new), setindex_unrolled(tail(old), tail(new), tail(switches))...
    else
        first(old), setindex_unrolled(tail(old), new, tail(switches))...
    end

struct Slices{Item,Dimensions,Whole,Alongs} <: AbstractArray{Item,Dimensions}
    whole::Whole
    alongs::Alongs
end
@inline Slices{Item,Dimensions}(
    whole::Whole,
    alongs::Alongs,
) where {Item,Dimensions,Whole,Alongs} = Slices{Item,Dimensions,Whole,Alongs}(whole, alongs)

@inline axes(slices::Slices) =
    getindex_unrolled(axes(slices.whole), map(not, slices.alongs))
@inline size(slices::Slices) = map(length, axes(slices))

@inline slice_index(slices, indices) =
    setindex_unrolled(axes(slices.whole), indices, map(not, slices.alongs))
@inline getindex(slices::Slices, indices::Int...) =
    view(slices.whole, slice_index(slices, indices)...)
@inline setindex!(slices::Slices, value, indices::Int...) =
    slices.whole[slice_index(slices, indices)...] = value

@inline axis_or_1(switch, axis) =
    if untyped(switch)
        axis
    else
        1
    end
"""
    Slices(whole, alongs::TypedBool...)

Slice `whole` into `view`s.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be replaced with `:` when slicing.

```jldoctest
julia> using JuliennedArrays

julia> whole = [1 2; 3 4];

julia> slices = Slices(whole, False(), True())
2-element Slices{SubArray{Int64,1,Array{Int64,2},Tuple{Int64,Base.OneTo{Int64}},true},1,Array{Int64,2},Tuple{False,True}}:
 [1, 2]
 [3, 4]

julia> slices[1] == whole[1, :]
true

julia> slices[1] = [2, 1];

julia> whole
2×2 Array{Int64,2}:
 2  1
 3  4

julia> larger = rand(5, 5, 5);

julia> larger_slices = Slices(larger, True(), False(), False());

julia> size(first(larger_slices))
(5,)
```
"""
Slices(whole::AbstractArray, alongs::TypedBool...) = Slices{
    typeof(@inbounds view(whole, map(axis_or_1, alongs, axes(whole))...)),
    length(getindex_unrolled(alongs, map(not, alongs))),
}(
    whole,
    alongs,
)

"""
    Slices(whole, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be replaced with `:` when slicing.

```jldoctest
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{Int64}, 2, 2, 2) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> s = Slices(input, 1, 3)
2-element Slices{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{True,False,True}}:
 [1 5; 2 6]
 [3 7; 4 8]

julia> map(sum, s)
2-element Array{Int64,1}:
 14
 22
```
"""
Slices(
    whole::AbstractArray{Item,NumberOfDimensions},
    alongs::Int...,
) where {Item,NumberOfDimensions} =
    Slices(whole, in_unrolled(as_vals(alongs...), ntuple(Val, NumberOfDimensions)...)...)
export Slices

struct Align{Item,Dimensions,Sliced,Alongs} <: AbstractArray{Item,Dimensions}
    slices::Sliced
    alongs::Alongs
end
@inline Align{Item,Dimensions}(
    slices::Sliced,
    alongs::Alongs,
) where {Item,Dimensions,Sliced,Alongs} =
    Align{Item,Dimensions,Sliced,Alongs}(slices, alongs)

@inline axes(aligned::Align) = setindex_unrolled(
    setindex_unrolled(aligned.alongs, axes(aligned.slices), map(not, aligned.alongs)),
    axes(first(aligned.slices)),
    aligned.alongs,
)
@inline size(aligned::Align) = map(length, axes(aligned))

@inline split_indices(aligned, indices) =
    getindex_unrolled(indices, map(not, aligned.alongs)),
    getindex_unrolled(indices, aligned.alongs)
@inline function getindex(aligned::Align, indices::Int...)
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...]
end
@inline function setindex!(aligned::Align, value, indices::Int...)
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...] = value
end

"""
    Align(slices, alongs::TypedBool...)

`Align` an array of arrays, all with the same size.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be taken up by the inner arrays. Inverse of [`Slices`](@ref).

```jldoctest
julia> using JuliennedArrays

julia> slices = [[1, 2], [3, 4]];

julia> aligned = Align(slices, False(), True())
2×2 Align{Int64,2,Array{Array{Int64,1},1},Tuple{False,True}}:
 1  2
 3  4

julia> aligned[1, :] == slices[1]
true

julia> aligned[1, 1] = 0;

julia> slices
2-element Array{Array{Int64,1},1}:
 [0, 2]
 [3, 4]
```
"""
@inline Align(
    slices::AbstractArray{<:AbstractArray{Item,InnerDimensions},OuterDimensions},
    alongs::TypedBool...,
) where {Item,InnerDimensions,OuterDimensions} =
    Align{Item,OuterDimensions + InnerDimensions}(slices, alongs)
export Align

"""
    Along(slices, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be taken up by the inner arrays.

```jldoctest
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{Int64}, 2, 2, 2) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> slices = collect(Slices(input, 1, 3))
2-element Array{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1}:
 [1 5; 2 6]
 [3 7; 4 8]

julia> Align(slices, 1, 3)
2×2×2 Align{Int64,3,Array{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1},Tuple{True,False,True}}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8
```
"""
Align(
    slices::AbstractArray{<:AbstractArray{Item,InnerDimensions},OuterDimensions},
    alongs::Int...,
) where {Item,InnerDimensions,OuterDimensions} = Align(
    slices,
    in_unrolled(
        as_vals(alongs...),
        ntuple(Val, InnerDimensions + OuterDimensions)...,
    )...,
)

end


text:
function createQuadraticModelT(qpdata; T = Float128, name = "qp_pb")
  return QuadraticModel(
    convert(Array{T}, qpdata.c),
    qpdata.qrows,
    qpdata.qcols,
    convert(Array{T}, qpdata.qvals),
    Arows = qpdata.arows,
    Acols = qpdata.acols,
    Avals = convert(Array{T}, qpdata.avals),
    lcon = convert(Array{T}, qpdata.lcon),
    ucon = convert(Array{T}, qpdata.ucon),
    lvar = convert(Array{T}, qpdata.lvar),
    uvar = convert(Array{T}, qpdata.uvar),
    c0 = T(qpdata.c0),
    x0 = zeros(T, length(qpdata.c)),
    name = name,
  )
end

@testset "multi_mode" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "Float16, Float32, Float128" begin
  qm128_1 = createQuadraticModelT(qps1, T = Float128)
  stats1 = ripqp(
    qm128_1,
    itol = InputTol(Float128, ϵ_rb32 = Float128(0.1), ϵ_rb64 = Float128(0.01)),
    iconf = InputConfig(mode = :multi, normalize_rtol = false),
    display = false,
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  for T ∈ [Float32, Float128]
    qmT_2 = createQuadraticModelT(qps2, T = T)
    stats2 = ripqp(qmT_2, display = false)
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
    @test stats2.status == :acceptable
  end

  qm16 = QuadraticModel(
    Float16.(c),
    Float16.(Q),
    A = Float16.(A),
    lcon = Float16.(b),
    ucon = Float16.(b),
    lvar = Float16.(l),
    uvar = Float16.(u),
    c0 = Float16(0.0),
    x0 = zeros(Float16, 3),
    name = "QM16",
  )
  stats_dense = ripqp(
    qm16,
    itol = InputTol(
      Float16,
      ϵ_pdd = Float16(1.0e-1),
      ϵ_rb = Float16(1.0e-1),
      ϵ_rc = Float16(1.0e-1),
    ),
    display = false,
  )
  @test isapprox(stats_dense.objective, 1.1249999990782493, atol = 1e-2)
  @test stats_dense.status == :acceptable
end


text:
using REALMSupport
using REALMSupport.RegisterCore: CenterIndexedArray
using Optim, Images, RegisterMismatch, FFTW
FFTW.set_num_threads(16)

filelist = readdir()
filelist = filter(x -> occursin("-fi.tif", x), filelist)
filelist = filter(x -> occursin("20220317_1um", x), filelist)

mm = CenterIndexedArray(zeros(241, 241, 241, size(filelist, 1)))
mm_r = []
img_org = load("20220317_1um_x-0_1_1-fi.tif")

for (i, j) in enumerate(axes(mm, 4))
    img = load(filelist[i])
    mm[:, :, :, j] = translate(img_org, img, (120, 120, 120))
    push!(mm_r, register_translate(img_org, img, (120, 120, 120)))
end

move = [
    20 20 40 40 0 0 -20 -20 -40 -40 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 20 20 40 40 0 -20 -20 -40 -40 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 20 40 40 0 0 -20 -20 -40 -40
]

#move = [20 20 20 40 40 40 0 0 0 -20 -20 -20 -40 -40 -40 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 20 20 40 40 40 0 0 0 -20 -20 -20 -40 -40 -40 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 20 20 40 40 40 0 0 0 -20 -20 -20 -40 -40 -40]
#move = [ 20 40 0 -20 -40 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 20 40 0 -20 -40 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 20 40 0 -20 -40]

matrix = Float64[-0.01 0 -1.84 0 1.56 0 1.198 0 -0.39]

result = translate_optim(mm, matrix, move; g_abstol = 1e-14)

mat = reshape(Optim.minimizer(result), 3, 3)
display(mat)
display((mat * move)')
display(mm_r)


text:
export RandomizedSystem, square_up

"""
    RandomizedSystem(F::Union{System,AbstractSystem}, k::Integer) <: AbstractSystem

Given a ``n × N`` system `F` with ``n > N`` this constructs the system
``\\mathfrak{R}(F; k)(x) = [I A]⋅F(x)`` where ``I`` is a ``k × k`` identity matrix and
``A`` is random complex ``k × n`` matrix. See Chapter 13.5 in [^SW05] for more details.

    RandomizedSystem(F::Union{System,AbstractSystem}, A::Matrix{ComplexF64})

Explicitly provide the used randomization matrix `A`.

[^SW05]: Sommese, A. J., & Wampler, C. W. (2005). The Numerical Solution of Systems of Polynomials Arising in Engineering and Science. World Scientific.
"""
struct RandomizedSystem{S<:AbstractSystem} <: AbstractSystem
    system::S
    A::Matrix{ComplexF64}
    u::Vector{ComplexF64}
    ū::Vector{ComplexDF64}
    U::Matrix{ComplexF64}
    taylor_U::Matrix{ComplexF64}
end

function RandomizedSystem(F::Union{AbstractSystem,System}, k::Integer)
    n, N = size(F)
    #TODO: Should the rows of [I A] be unitary?
    A = randn(ComplexF64, k, n - k)
    RandomizedSystem(F, A)
end
RandomizedSystem(F::System, A::Matrix{ComplexF64}) = RandomizedSystem(ModelKitSystem(F), A)
function RandomizedSystem(F::AbstractSystem, A::Matrix{ComplexF64})
    n, N = size(F)
    n > N || throw(ArgumentError("Then given system is not overdetermined."))
    u = zeros(ComplexF64, n)
    ū = zeros(ComplexDF64, n)
    U = zeros(ComplexF64, n, N)
    taylor_U = zeros(ComplexF64, n, 5)
    RandomizedSystem(F, A, u, ū, U, taylor_U)
end

function Base.show(io::IO, mime::MIME"text/plain", F::RandomizedSystem)
    println(io, typeof(F), ":")
    println(io, "A:")
    show(io, mime, F.A)
    println(io, "\n\nF:")
    show(io, mime, F.system)
end

Base.size(F::RandomizedSystem) = (size(F.A, 1), last(size(F.system)))
ModelKit.variables(F::RandomizedSystem) = variables(F.system)
ModelKit.parameters(F::RandomizedSystem) = parameters(F.system)
ModelKit.variable_groups(F::RandomizedSystem) = variable_groups(F.system)

"""
    square_up(F::Union{System, AbstractSystem})

Creates the [`RandomizedSystem`](@ref) ``\\mathfrak{R}(F(x); N)`` where ``N`` is the number
of variables of `F`.
"""
square_up(F::Union{AbstractSystem,System}) = RandomizedSystem(F, last(size(F)))

function randomize!(u, A, v::AbstractVector, ncols = 1)
    n, m = size(A, 1), length(v)

    @inbounds for i = 1:n
        u[i] = v[i]
    end
    @inbounds for j = 1:(m-n)
        vnj = v[n+j]
        for i = 1:n
            u[i] += A[i, j] * vnj
        end
    end
    u
end
function randomize!(U, A, V::AbstractMatrix, ncols = size(A, 1))
    n, m = size(A, 1), size(V, 1)

    for j = 1:ncols, i = 1:n
        U[i, j] = V[i, j]
    end
    for j = 1:ncols, k = 1:(m-n), i = 1:n
        U[i, j] += A[i, k] * V[n+k, j]
    end
    U
end

(F::RandomizedSystem)(x, p = nothing) = [LA.I F.A] * F.system(x, p)

function evaluate!(u, F::RandomizedSystem, x, p = nothing)
    if eltype(x) isa ComplexDF64
        evaluate!(F.ū, F.system, x, p)
        randomize!(u, F.A, F.ū)
    else
        evaluate!(F.u, F.system, x, p)
        randomize!(u, F.A, F.u)
    end
    u
end

function evaluate_and_jacobian!(u, U, F::RandomizedSystem, x, p = nothing)
    evaluate_and_jacobian!(F.u, F.U, F.system, x, p)
    randomize!(u, F.A, F.u)
    randomize!(U, F.A, F.U)
end

function taylor!(u::AbstractVector, v::Val, F::RandomizedSystem, tx, p = nothing)
    taylor!(F.u, v, F.system, tx, p)
    randomize!(u, F.A, F.u)
end
function taylor!(
    U::AbstractMatrix,
    v::Val{N},
    F::RandomizedSystem,
    tx::TaylorVector,
    p = nothing,
) where {N}
    taylor!(F.taylor_U, v, F.system, tx, p)
    randomize!(U, F.A, F.taylor_U, N + 1)
end


text:
using Colors, Gadfly, RDatasets

set_default_plot_size(10cm, 8cm)

iris = dataset("datasets", "iris")
palettef = Scale.lab_gradient(range(HSV(250,1,1), stop=HSV(0,1,1), length=100)...)
colkey = Guide.colorkey(title="Petal\nlength")

plot(iris, x=:SepalLength, y=:SepalWidth, color=:PetalLength,
    shape=:Species, Guide.shapekey(nothing), colkey,
    Scale.color_continuous(colormap=palettef, minvalue=0, maxvalue=7)
)


text:
using BandedMatrices
using LinearAlgebra

function diffs(a, b)
    n = length(a)
    d = Float64.(abs.(a - b))
    a1 = sum(d) / n
    r1 = sum(d ./ abs.(a)) / n
    a2 = sqrt(sum(d .* d) / n)
    r2 = sqrt(sum(d .* d ./ abs2.(a)) / n)
    a0 = maximum(d)
    r0 = maximum(d ./ abs.(a))
    a1, a2, a0, r1, r2, r0
end

function evaluate(n::Integer, T=Float64)

    M = BandedMatrix(0 => 6ones(n), 1 => -4ones(n-1), 2 => 1ones(n-2));
    MC = BandedMatrix(0 => 6ones(n), 1 => -4ones(n-1) .+ 0.1im, 2 => 1ones(n-2));
    M = T.(M)
    MC = Complex{T}.(MC)

    println("evaluation times with 5-banded test matrixes of dimension $n")
    println("LAPACK $T")
    Mb = @btime BandedMatrices.tridiagonalize!((Symmetric(copy($M))));

    println("generic $T")
    Mt = @btime tridiagonalize(Hermitian($M));

    println("generic with LinearAlgebra.givensAlgorithm $T")
    MtGA = @btime BandedMatrices.tridiagonalizeGA(Hermitian($M));

    println("generic Complex{$T}")
    Mtc = @btime tridiagonalize(Hermitian($MC));

    println("generic with LinearAlgebra.givensAlgorithm Complex{$T}")
    MtcGA = @btime BandedMatrices.tridiagonalizeGA(Hermitian($MC));

    println("generic BigFloat($(precision(BigFloat)))")
    Mbig = @time tridiagonalize(Hermitian(big.(M)));
    println("generic complex BigFloat($(precision(BigFloat)))")
    MCbig = @time tridiagonalize(Hermitian(big.(MC)));

    et = eigvals(Mt);
    etGA = eigvals(MtGA);
    etc = eigvals(Mtc);
    etcGA = eigvals(MtcGA);
    eb = eigvals(Mb);
    ebib = eigvals(Mbig);
    ebig = ebib;
    eCbig = eigvals(MCbig);

    println("absolute and relative errors")
    println("(mean / square-mean / max ) * (abs / rel)")

    println("LAPACK $T")
    println(diffs(eb, ebig))
    println("generic $T")
    println(diffs(et, ebig))
    println("generic with LinearAlgebra.givensAlgorithm $T")
    println(diffs(etGA, ebig))
    println("generic Complex{$T}")
    println(diffs(etc, eCbig))
    println("generic with LinearAlgebra.givensAlgorithm Complex{$T}")
    println(diffs(etcGA, eCbig))
end

evaluate( 1000);


text:
include("testga.jl")

test_serial()
test_parallel()


text:
module DistributedReinforcementLearning

using Distributed
using ReinforcementLearningBase
using ReinforcementLearningCore

include("actor_model.jl")
include("core.jl")
include("extensions.jl")

end


text:
@parameters t
D = Differential(t)
@variables x(t) # integrator

@named plf = plant_face_fluct()
#structural_simplify(plf)

function rep_plf(pl;name, simplify=true)
    @parameters t
    D = Differential(t)
        sts = @variables x(t)
    # integrate of i_L to check annual litter inputs
    tmp = compose(ODESystem([
        D(x) ~ pl.i_L
        ],t,sts,[];name), pl)
    simplify ? structural_simplify(tmp) : tmp
 end
 @named plf_rep = rep_plf(plf)
# @named plf_rep = embed_system(plf)
 equations(plf_rep)

p_plf = Dict(
    plf.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    plf.β_Ni0 => 25,
    plf.i_IN0 => 0,   ##<< input of mineral N,
    plf.t1 => 0.0,
    plf.t2 => 3.0,
    #
    plf.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    plf.β_Ni0 => 25,
    plf.i_IN0 => 0,   ##<< input of mineral N,
    plf.β_Pi0 => 25*20, ## leaf litter N:P ~20(massratio Kang10)
    plf.i_IP0 => 0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
    plf.s_EP0 => 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
)
u_plf = Dict(
    plf.Lagr => 0.2,
    x => 0.0,
)
#tspan = (-500.0, 120.0)
tspan = (-5.0, 5.0)
prob = ODEProblem(plf_rep, u_plf, tspan, p_plf)
#sol = solve(prob, saveat = 0:1:tspan[2])
#solver = CompositeAlgorithm((Tsit5(),Vern7()), (integrator) -> (Int(integrator.dt<0.001) + 1))
solver = Vern7() # implicit method
# solver = AutoVern9(Rodas4())
# solver = Vern9()
#using  deSolveDiffEq; solver = deSolveDiffEq.lsoda()

# solver = Tsit5() # takes forever and produces wrong solution
sol = solve(prob, solver,
    reltol = 1e-5,  # simulate litterfall with higher accuracy
    #progress = true, progress_steps = 1,
    #saveat = 0.005 # uncomment for proper plotting
);
# sol = sol_euler = solve(prob, Euler(),
#     dt = 0.005,
#     progress = true, progress_steps = 1,
#     #saveat = 0.005
# );

i_plot = () -> begin
    # using Plots
    ts = tspan
    Plots.plot(sol, tspan=ts, vars=[plf.i_L, plf.i_L_annual], xlab="time (yr)", ylab="litter input (g/m2/yr)")
    Plots.plot(sol, tspan=ts, vars=[plf.β_Ni], xlab="time (yr)", ylab="CN litter input (g/g)")
    Plots.plot(sol, tspan=ts, vars=[plf.inc_period], xlab="time (yr)", ylab="inc_period")
    series_sol(sol, )

    ts = (-2,5)
    ts = (-0.1,+0.1)
    ts = (0.7,1)
    ts = (-0.3,0)
    ts = (0,10) .+ tspan[1]
    ts = tspan
    Plots.plot(sol, tspan=ts, vars=[plf.i_Lagr, plf.dec_Lagr, plf.d_Lagr])

    x = (9.5/12:1/(12*10):11.5/12) .+ 1
    Plots.scatter(x, first.(sol.(x)))
    plot(sol, tspan=extrema(x), vars=[plf.i_Lagr, plf.dec_Lagr])

    plot(sol, vars=[plf.i_L / plf.β_Ni])
end

@testset "integration of litterfall over one year" begin
    isapprox(sol(0, idxs=x) - sol(-1, idxs=x), 400, atol=0.05)
    isapprox(sol(5, idxs=x) - sol(4, idxs=x), 480, atol=0.05)
end;

i_plot_integrated_annual = () -> begin
    Plots.plot(sol, vars=[x])
    ts = range(first(tspan)+1, last(tspan), length=200)
    xann = sol.(ts, idxs=x) .- sol.(ts .- 1.0, idxs=x)
    Plots.plot(ts, xann)
end

# dlit = shifloNormal(10/12, 11/12)
# #plot(dlit, xlim=(0,1))
# plot(dlit)
# ts = 0:1/24:3
# t = 10.5*1/12Sophienstraße 2, 07743 Jena
# ts = (9.5:0.05:11.5)/12
# plot(ts*12, iL.(ts,1))

# # check that integral matches annual litterfall
# using QuadGK
# lann = 3.2
# f1(x) = iL(x,lann)
# @test isapprox(first(quadgk(f1, 0, 1, rtol=1e-3)), lann, rtol=1e-2)

i_inspect_dlit = () -> begin
    autumn_start=8.5/12; autumn_end=11.5/12
    d_lit_agr = LocationScale(autumn_start, autumn_end - autumn_start,
        #     LogitNormal(0,sqrt(2)),
        fit_mode_flat(LogitNormal, 0.3; peakedness = 3)
        #    Uniform()
    )
    xs = 0.6:0.005:1
    Plots.plot(xs .* 12 .+ 1, pdf.(d_lit_agr, xs))
end


text:
# LOGICCIRCUITS LIBRARY ROOT

module LogicCircuits

# USE EXTERNAL MODULES

using Reexport

include("Utils/Utils.jl")
include("Data/Data.jl")

@reexport using .Data
@reexport using .Utils

# INCLUDE CHILD MODULES
include("Logical/Logical.jl")
include("IO/IO.jl")

# USE CHILD MODULES (in order to re-export some functions)
@reexport using .Logical
@reexport using .IO

end


text:
suite["flow"] = BenchmarkGroup()

edgs = [
  (1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
  (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
  (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
]

flow_graph = DiGraph(8)
capacity_matrix = zeros(Int,8,8)
for e in edgs
    u,v,f = e
    add_edge!(flow_graph,u,v)
    capacity_matrix[u,v] = f
end

suite["flow"]["push_relabel"] = @benchmarkable maximum_flow($flow_graph, 1, 8
                    , $capacity_matrix, algorithm=PushRelabelAlgorithm())
suite["flow"]["dinic"] = @benchmarkable maximum_flow($flow_graph, 1, 8
                    , $capacity_matrix, algorithm=DinicAlgorithm())


text:
function setup_warm_start()
    options = Dict(:verbose => false,
                   :eps_abs => 1e-08,
                   :eps_rel => 1e-08,
                   :polish => false,
                   :adaptive_rho => false,
                   :check_termination => 1)
    return options
end

tol = 1e-5

@testset "warm_start" begin

    @testset "warm_start_problem" begin
        srand(1)

        n = 100
        m = 200
        P = sprandn(n, n, 0.9)
        P = P' * P
        q = randn(n)
        A = sprandn(m, n, 0.9)
        u = rand(m) * 2
        l = -rand(m) * 2

        options = setup_warm_start()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        # Store optimal values
        x_opt = results.x
        y_opt = results.y
        tot_iter = results.info.iter

        # Warm start with zeros to check if number of iterations is the same
        OSQP.warm_start!(model, x=zeros(n), y=zeros(m))
        results = OSQP.solve!(model)
        @test results.info.iter == tot_iter

        # Warm start with optimal values and check that total number of iter < 10
        OSQP.warm_start!(model, x=x_opt, y=y_opt)
        results = OSQP.solve!(model)
        @test results.info.iter <= 10

    end

end


text:
module Function
using  ...PatternStructure.Trees
using  ...PatternStructure.Checks
using  ...PatternStructure.Special
using  ...PatternStructure.SlurpTypes
using  ...Analyzer.SlurpOptimizations
using  ...Helper
export analyze

#-----------------------------------------------------------------------------------
# Analysis state
#-----------------------------------------------------------------------------------

struct AnalysisState
  tree    :: PatternTree
  literal :: Bool
  mod     :: Module
  consts  :: Set{Symbol}
end
const newstate = AnalysisState
AnalysisState(state, tree::PatternTree) = newstate(tree, state.literal, state.mod, state.consts)
AnalysisState(state, literal::Bool)     = newstate(state.tree, literal, state.mod, state.consts)

#-----------------------------------------------------------------------------------
# analyze function: transform expressions into patterns
#-----------------------------------------------------------------------------------

function analyze(ex, mod=current_module())
  root   = PatternRoot()
  consts = Set{Symbol}()
  analyze!(ex, newstate(root, false, mod, consts))
  addbindings!(root)

  return root
end

#-----------------------------------------------------------------------------------
# analyze!: constructs patterns from expression and inserts them in a pattern tree
#-----------------------------------------------------------------------------------

function analyze!(ex::Symbol, state)
  is_binding_name(ex) && !state.literal ?
     newleaf!(newbinding(ex), state.tree) :
     newleaf!(equalsto(ex),   state.tree)
  #--
end

function analyze!(ex::Expr, state)
  is_special_expr(ex) && (return analyze_special!(ex, state))
  ex.head == :line && (return analyze!(LineNumberNode(0), state))

  step = getstep(ex.head)
  head = ex.head in [:kw, :(=)] ? :assign : ex.head
  node = newnode!(ExprHead(head), step, state.tree)
  args = step(ex)

  analyze_args!(args, node, state)
end

function analyze!(ex::QuoteNode, state)
  step = QuoteStep()
  node = newnode!(ExprHead(:quote), step, state.tree)
  analyze!(ex.value, newstate(state, node))
end

function analyze!(ex::LineNumberNode, state)
  analyze_special!(:(:P{$Helper.is_line_number}), state)
end

function analyze!(ex, state)
  newleaf!(equalsto(ex), state.tree)
end

#-----------------------------------------------------------------------------------
# analyze_args!: creates patterns from a vector of expressions
#-----------------------------------------------------------------------------------

function analyze_args!(args, node, state)
  nstate = newstate(state, node)
  for i in eachindex(args)
    analyze!(args[i], nstate)
    union!(node.bindings, bindings(node.children[i]))
  end
  optimize_slurps!(node, state)
end

#-----------------------------------------------------------------------------------
# analyze_special: patterns of the form :X{ys...} have special meaning
#-----------------------------------------------------------------------------------

function analyze_special!(ex, state)
  @assert length(ex.args) > 1 "Special syntax with no arguments is not supported."

  head = special_name(ex.args[1])
  args = ex.args[2:end]

  if head in [:?, :*]
    slurptype = head == :? ? GenericLazySlurp() : GenericGreedySlurp()

    node = newnode!(slurptype, SlurpStep(), state.tree)
    analyze_args!(args, node, state)

  elseif head == :binding
    @assert length(args) == 1    "the :binding pattern only accepts one argument"
    @assert isa(args[1], Symbol) "the :binding pattern only accepts a Symbol argument"
    newleaf!(newbinding(args[1]),  state.tree)

  elseif head == :autobinding
    @assert state.literal "Invalid :autobinding pattern. Only use it inside a :literal pattern."
    @assert length(args) == 1 "the :autobinding pattern only accepts one argument."
    analyze!(args[1], newstate(state, false))

  elseif head == :literal
    @assert !state.literal "Invalid :literal pattern. Don't use it inside another :literal pattern."
    @assert length(args) == 1 "the :literal pattern only accepts one argument."
    analyze!(args[1], newstate(state, true))

  elseif head == :predicate
    args = assertation_args(args)
    pred = eval(state.mod, args[2])
    gate = PatternGate(PredicateCheck(pred), depth(state.tree))
    insert!(state.tree, gate)
    analyze!(args[1], newstate(state, gate))

  elseif head == :type
    args = assertation_args(args)
    typ  = eval(state.mod, args[2])
    gate = PatternGate(TypeCheck{typ}(), depth(state.tree))
    insert!(state.tree, gate)
    analyze!(args[1], newstate(state, gate))

  elseif head == :equals
    args = assertation_args(args)
    val  = eval(state.mod, args[2])
    gate = PatternGate(EqualityCheck(val), depth(state.tree))
    insert!(state.tree, gate)
    analyze!(args[1], newstate(state, gate))

  elseif head == :iterable
    node = newnode!(ExprHead(:iterable), IterStep(), state.tree)
    analyze_args!(args, node, state)

  elseif head == :raw
    analyze!(Expr(QuoteStep()(args[1])[1], args[2:end]...), state)

  end
end

#-----------------------------------------------------------------------------------
# Find and add bindings
#-----------------------------------------------------------------------------------

function addbindings!(tree::PatternGate)
  addbindings!(tree.child)
  union!(tree.bindings, bindings(tree.child))
end

function addbindings!(tree::PatternNode)
  for child in tree.children
    addbindings!(child)
    union!(tree.bindings, bindings(child))
  end
end

function addbindings!(tree::PatternRoot)
  addbindings!(tree.child)
end

function addbindings!(tree::PatternLeaf)
end

#-----------------------------------------------------------------------------------
# Utility functions
#-----------------------------------------------------------------------------------

function getstep(head)
  head == :quote ? QuoteStep() :
  head == :block ? BlockStep() :
                  ArgsStep()
end

function is_macro_name(x::Symbol)
  string(x)[1] == '@'
end

function newbinding(ex)
  name = is_macro_name(ex) ? Symbol(string(ex)[2:end]) : ex
  Binding(name)
end

equalsto(x) = EqualityCheck(x)

assertation_args(args) =
  length(args) == 1 ? [gensym("x"), args[1]] : args

#-----------------------------------------------------------------------------------

end


text:
using NonUniformRandomVariateGeneration
using Distributions
using BenchmarkTools
using Random
using StatsBase

function benchNative1(n::Int64, ps::Vector{Float64}, N::Int64)
  v::Vector{Int64} = Vector{Int64}(undef, n)
  for i in 1:N
    sampleCategoricalSorted!(v, ps)
  end
end

function benchNative2(n::Int64, ps::Vector{Float64}, N::Int64)
  v::Vector{Int64} = Vector{Int64}(undef, n)
  scratch1::Vector{Float64} = Vector{Float64}(undef, n)
  scratch2::Vector{Float64} = Vector{Float64}(undef, length(ps))
  for i in 1:N
    sampleCategoricalSorted!(v, ps, scratch1, scratch2)
  end
  return v
end

function benchDistributions(n::Int64, ps::Vector{Float64}, N::Int64)
  multinom = Distributions.Multinomial(n, ps)
  v::Vector{Int64} = Vector{Int64}(undef, n)
  a = Vector{Int64}(undef, length(ps))
  for i in 1:length(ps)
    a[i] = i
  end
  ws = StatsBase.weights(ps)
  for i in 1:N
    Distributions.sample!(a, ws, v)
  end
  return v
end

n = 10000
m = 1000
ps = exp.(randexp(m))
ps ./= sum(ps)
N = 1024

@btime benchNative1(n, ps, N)
@btime benchNative2(n, ps, N)
@btime benchDistributions(n, ps, N)


text:
using LinearAlgebra
function mwe()
  local Q::Matrix{Float64}

  #Q = Matrix{Float64}(undef, 800, 800)
  X = rand(200,13)
  f = qr(X)
  @time Q = Matrix(f.Q[:,:])

  X = rand(400,13)
  f = qr(X)
  @time Q = Matrix(f.Q[:,:])

  X = rand(600,13)
  f = qr(X)
  @time Q = Matrix(f.Q[:,:])

  X = rand(800,13)
  f = qr(X)
  @time Q = Matrix(f.Q[:,:])
end

@time mwe()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module GradientTest

import Calculus

using Base.Test
using ForwardDiff
using ForwardDiff: Dual, Tag
using StaticArrays

include(joinpath(dirname(@__FILE__), "utils.jl"))

##################
# hardcoded test #
##################

f = DiffBase.rosenbrock_1
x = [0.1, 0.2, 0.3]
v = f(x)
g = [-9.4, 15.6, 52.0]

for c in (1, 2, 3), tag in (nothing, f)
    println("  ...running hardcoded test with chunk size = $c and tag = $tag")
    cfg = ForwardDiff.GradientConfig(tag, x, ForwardDiff.Chunk{c}())

    @test eltype(cfg) == Dual{typeof(Tag(typeof(tag), eltype(x))), eltype(x), c}

    @test isapprox(g, ForwardDiff.gradient(f, x, cfg))
    @test isapprox(g, ForwardDiff.gradient(f, x))

    out = similar(x)
    ForwardDiff.gradient!(out, f, x, cfg)
    @test isapprox(out, g)

    out = similar(x)
    ForwardDiff.gradient!(out, f, x)
    @test isapprox(out, g)

    out = DiffBase.GradientResult(x)
    ForwardDiff.gradient!(out, f, x, cfg)
    @test isapprox(DiffBase.value(out), v)
    @test isapprox(DiffBase.gradient(out), g)

    out = DiffBase.GradientResult(x)
    ForwardDiff.gradient!(out, f, x)
    @test isapprox(DiffBase.value(out), v)
end

########################
# test vs. Calculus.jl #
########################

for f in DiffBase.VECTOR_TO_NUMBER_FUNCS
    v = f(X)
    g = ForwardDiff.gradient(f, X)
    @test isapprox(g, Calculus.gradient(f, X), atol=FINITEDIFF_ERROR)
    for c in CHUNK_SIZES, tag in (nothing, f)
        println("  ...testing $f with chunk size = $c and tag = $tag")
        cfg = ForwardDiff.GradientConfig(tag, X, ForwardDiff.Chunk{c}())

        out = ForwardDiff.gradient(f, X, cfg)
        @test isapprox(out, g)

        out = similar(X)
        ForwardDiff.gradient!(out, f, X, cfg)
        @test isapprox(out, g)

        out = DiffBase.GradientResult(X)
        ForwardDiff.gradient!(out, f, X, cfg)
        @test isapprox(DiffBase.value(out), v)
        @test isapprox(DiffBase.gradient(out), g)
    end
end

##########################################
# test specialized StaticArray codepaths #
##########################################

x = rand(3, 3)
sx = StaticArrays.SArray{Tuple{3,3}}(x)
out = similar(x)
actual = ForwardDiff.gradient(prod, x)

@test ForwardDiff.gradient(prod, sx) == actual

ForwardDiff.gradient!(out, prod, sx)

@test out == actual

result = DiffBase.GradientResult(x)
sresult = DiffBase.GradientResult(sx)

ForwardDiff.gradient!(result, prod, x)
ForwardDiff.gradient!(sresult, prod, sx)

@test DiffBase.value(sresult) == DiffBase.value(result)
@test DiffBase.gradient(sresult) == DiffBase.gradient(result)

end # module


text:
function f(path)
    open(path) do f
        N = 0
        corrects = 0
        for line in readlines(f)
            tokens = split(line)
            if tokens[4] == tokens[5] && tokens[6] != -1
                corrects += 1
            end
            N += 1
        end

        println("Result for $path: ", corrects/N)
    end
end

f("92.jl.85.out")
f("92.jl.90.out")

#=
julia> @time include("./93.jl")
Result for 92.jl.85.out: 0.15019762845849802
Result for 92.jl.90.out: 0.2134387351778656
=#


text:
"""
$(TYPEDSIGNATURES)

Receives a tuple of arrays representing categorical factor levels, and returns a
`Base.Iterators.ProductIterator`.   This allows  full  factorial  designs to  be
arbitrarily large and  only be computed as needed.  To  compute an explicit full
factorial design, use [`explicit_fullfactorial`](@ref).

```jldoctest
julia> fullfactorial(Tuple([-1, 1] for i = 1:10))
Base.Iterators.ProductIterator{NTuple{10,Array{Int64,1}}}(([-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1]))

```
"""
function fullfactorial(factors::Tuple)
    Base.Iterators.product(factors...)
end

"""
$(TYPEDSIGNATURES)

Receives  a  `Base.Iterators.ProductIterator`  and  computes  an  explicit  full
factorial design. The generated array is exponentially large.

```jldoctest
julia> explicit_fullfactorial(fullfactorial(([-1, 1], [:a, :b, :c])))
6×2 Array{Any,2}:
 -1  :a
  1  :a
 -1  :b
  1  :b
 -1  :c
  1  :c

```
"""
function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator)
    permutedims(foldl((x, y) -> hcat(collect(x), collect(y)), iterator))
end

"""
$(TYPEDSIGNATURES)

Receives a tuple of arrays  representing categorical factor levels, and computes
an explicit full factorial design. The generated array is exponentially large.

```jldoctest
julia> explicit_fullfactorial(([-1, 1], [:a, :b, :c], [1, 2]))
12×3 Array{Any,2}:
 -1  :a  1
  1  :a  1
 -1  :b  1
  1  :b  1
 -1  :c  1
  1  :c  1
 -1  :a  2
  1  :a  2
 -1  :b  2
  1  :b  2
 -1  :c  2
  1  :c  2

```
"""
function explicit_fullfactorial(factors::Tuple)
    explicit_fullfactorial(fullfactorial(factors))
end


text:
# Julia wrapper for header: /usr/local/include/gdalwarper.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

"""
    GDALWarpNoDataMasker(void * pMaskFuncArg,
                         int nBandCount,
                         GDALDataType eType,
                         int,
                         int,
                         int nXSize,
                         int nYSize,
                         GByte ** ppImageData,
                         int bMaskIsFloat,
                         void * pValidityMask,
                         int * pbOutAllValid) -> CPLErr
"""
function GDALWarpNoDataMasker(pMaskFuncArg, nBandCount::Cint, eType::GDALDataType, nXOff::Cint, nYOff::Cint, nXSize::Cint, nYSize::Cint, papabyImageData, bMaskIsFloat::Cint, pValidityMask, pbOutAllValid)
    ccall((:GDALWarpNoDataMasker, libgdal), CPLErr, (Ptr{Void}, Cint, GDALDataType, Cint, Cint, Cint, Cint, Ptr{Ptr{GByte}}, Cint, Ptr{Void}, Ptr{Cint}), pMaskFuncArg, nBandCount, eType, nXOff, nYOff, nXSize, nYSize, papabyImageData, bMaskIsFloat, pValidityMask, pbOutAllValid)
end

"""
    GDALWarpDstAlphaMasker(void * pMaskFuncArg,
                           int nBandCount,
                           GDALDataType,
                           int nXOff,
                           int nYOff,
                           int nXSize,
                           int nYSize,
                           GByte **,
                           int bMaskIsFloat,
                           void * pValidityMask) -> CPLErr
"""
function GDALWarpDstAlphaMasker(pMaskFuncArg, nBandCount::Cint, eType::GDALDataType, nXOff::Cint, nYOff::Cint, nXSize::Cint, nYSize::Cint, arg1, bMaskIsFloat::Cint, pValidityMask)
    ccall((:GDALWarpDstAlphaMasker, libgdal), CPLErr, (Ptr{Void}, Cint, GDALDataType, Cint, Cint, Cint, Cint, Ptr{Ptr{GByte}}, Cint, Ptr{Void}), pMaskFuncArg, nBandCount, eType, nXOff, nYOff, nXSize, nYSize, arg1, bMaskIsFloat, pValidityMask)
end

"""
    GDALWarpSrcAlphaMasker(void * pMaskFuncArg,
                           int,
                           GDALDataType,
                           int nXOff,
                           int nYOff,
                           int nXSize,
                           int nYSize,
                           GByte **,
                           int bMaskIsFloat,
                           void * pValidityMask,
                           int * pbOutAllOpaque) -> CPLErr
"""
function GDALWarpSrcAlphaMasker(pMaskFuncArg, nBandCount::Cint, eType::GDALDataType, nXOff::Cint, nYOff::Cint, nXSize::Cint, nYSize::Cint, arg1, bMaskIsFloat::Cint, pValidityMask, pbOutAllOpaque)
    ccall((:GDALWarpSrcAlphaMasker, libgdal), CPLErr, (Ptr{Void}, Cint, GDALDataType, Cint, Cint, Cint, Cint, Ptr{Ptr{GByte}}, Cint, Ptr{Void}, Ptr{Cint}), pMaskFuncArg, nBandCount, eType, nXOff, nYOff, nXSize, nYSize, arg1, bMaskIsFloat, pValidityMask, pbOutAllOpaque)
end

"""
    GDALWarpSrcMaskMasker(void * pMaskFuncArg,
                          int,
                          GDALDataType,
                          int nXOff,
                          int nYOff,
                          int nXSize,
                          int nYSize,
                          GByte **,
                          int bMaskIsFloat,
                          void * pValidityMask) -> CPLErr
"""
function GDALWarpSrcMaskMasker(pMaskFuncArg, nBandCount::Cint, eType::GDALDataType, nXOff::Cint, nYOff::Cint, nXSize::Cint, nYSize::Cint, arg1, bMaskIsFloat::Cint, pValidityMask)
    ccall((:GDALWarpSrcMaskMasker, libgdal), CPLErr, (Ptr{Void}, Cint, GDALDataType, Cint, Cint, Cint, Cint, Ptr{Ptr{GByte}}, Cint, Ptr{Void}), pMaskFuncArg, nBandCount, eType, nXOff, nYOff, nXSize, nYSize, arg1, bMaskIsFloat, pValidityMask)
end

"""
    GDALWarpCutlineMasker(void * pMaskFuncArg,
                          int,
                          GDALDataType,
                          int nXOff,
                          int nYOff,
                          int nXSize,
                          int nYSize,
                          GByte **,
                          int bMaskIsFloat,
                          void * pValidityMask) -> CPLErr
"""
function GDALWarpCutlineMasker(pMaskFuncArg, nBandCount::Cint, eType::GDALDataType, nXOff::Cint, nYOff::Cint, nXSize::Cint, nYSize::Cint, arg1, bMaskIsFloat::Cint, pValidityMask)
    ccall((:GDALWarpCutlineMasker, libgdal), CPLErr, (Ptr{Void}, Cint, GDALDataType, Cint, Cint, Cint, Cint, Ptr{Ptr{GByte}}, Cint, Ptr{Void}), pMaskFuncArg, nBandCount, eType, nXOff, nYOff, nXSize, nYSize, arg1, bMaskIsFloat, pValidityMask)
end

"""
    GDALCreateWarpOptions() -> GDALWarpOptions *

Create a warp options structure.
"""
function GDALCreateWarpOptions()
    ccall((:GDALCreateWarpOptions, libgdal), Ptr{GDALWarpOptions}, ())
end

"""
    GDALDestroyWarpOptions(GDALWarpOptions * psOptions) -> void

Destroy a warp options structure.
"""
function GDALDestroyWarpOptions(arg1)
    ccall((:GDALDestroyWarpOptions, libgdal), Void, (Ptr{GDALWarpOptions},), arg1)
end

"""
    GDALCloneWarpOptions(const GDALWarpOptions * psSrcOptions) -> GDALWarpOptions *

Clone a warp options structure.
"""
function GDALCloneWarpOptions(arg1)
    ccall((:GDALCloneWarpOptions, libgdal), Ptr{GDALWarpOptions}, (Ptr{GDALWarpOptions},), arg1)
end

"""
    GDALSerializeWarpOptions(const GDALWarpOptions * psWO) -> CPLXMLNode *
"""
function GDALSerializeWarpOptions(arg1)
    ccall((:GDALSerializeWarpOptions, libgdal), Ptr{CPLXMLNode}, (Ptr{GDALWarpOptions},), arg1)
end

"""
    GDALDeserializeWarpOptions(CPLXMLNode * psTree) -> GDALWarpOptions *
"""
function GDALDeserializeWarpOptions(arg1)
    ccall((:GDALDeserializeWarpOptions, libgdal), Ptr{GDALWarpOptions}, (Ptr{CPLXMLNode},), arg1)
end

"""
    GDALReprojectImage(GDALDatasetH hSrcDS,
                       const char * pszSrcWKT,
                       GDALDatasetH hDstDS,
                       const char * pszDstWKT,
                       GDALResampleAlg eResampleAlg,
                       double dfWarpMemoryLimit,
                       double dfMaxError,
                       GDALProgressFunc pfnProgress,
                       void * pProgressArg,
                       GDALWarpOptions * psOptions) -> CPLErr

Reproject image.

### Parameters
* **hSrcDS**: the source image file.
* **pszSrcWKT**: the source projection. If NULL the source projection is read from from hSrcDS.
* **hDstDS**: the destination image file.
* **pszDstWKT**: the destination projection. If NULL the destination projection will be read from hDstDS.
* **eResampleAlg**: the type of resampling to use.
* **dfWarpMemoryLimit**: the amount of memory (in bytes) that the warp API is allowed to use for caching. This is in addition to the memory already allocated to the GDAL caching (as per GDALSetCacheMax()). May be 0.0 to use default memory settings.
* **dfMaxError**: maximum error measured in input pixels that is allowed in approximating the transformation (0.0 for exact calculations).
* **pfnProgress**: a GDALProgressFunc() compatible callback function for reporting progress or NULL.
* **pProgressArg**: argument to be passed to pfnProgress. May be NULL.
* **psOptions**: warp options, normally NULL.

### Returns
CE_None on success or CE_Failure if something goes wrong.
"""
function GDALReprojectImage(hSrcDS::GDALDatasetH, pszSrcWKT, hDstDS::GDALDatasetH, pszDstWKT, eResampleAlg::GDALResampleAlg, dfWarpMemoryLimit::Cdouble, dfMaxError::Cdouble, pfnProgress::GDALProgressFunc, pProgressArg, psOptions)
    ccall((:GDALReprojectImage, libgdal), CPLErr, (GDALDatasetH, Cstring, GDALDatasetH, Cstring, GDALResampleAlg, Cdouble, Cdouble, GDALProgressFunc, Ptr{Void}, Ptr{GDALWarpOptions}), hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, eResampleAlg, dfWarpMemoryLimit, dfMaxError, pfnProgress, pProgressArg, psOptions)
end

"""
    GDALCreateAndReprojectImage(GDALDatasetH hSrcDS,
                                const char * pszSrcWKT,
                                const char * pszDstFilename,
                                const char * pszDstWKT,
                                GDALDriverH hDstDriver,
                                char ** papszCreateOptions,
                                GDALResampleAlg eResampleAlg,
                                double dfWarpMemoryLimit,
                                double dfMaxError,
                                GDALProgressFunc pfnProgress,
                                void * pProgressArg,
                                GDALWarpOptions * psOptions) -> CPLErr

Reproject an image and create the target reprojected image.
"""
function GDALCreateAndReprojectImage(hSrcDS::GDALDatasetH, pszSrcWKT, pszDstFilename, pszDstWKT, hDstDriver::GDALDriverH, papszCreateOptions, eResampleAlg::GDALResampleAlg, dfWarpMemoryLimit::Cdouble, dfMaxError::Cdouble, pfnProgress::GDALProgressFunc, pProgressArg, psOptions)
    ccall((:GDALCreateAndReprojectImage, libgdal), CPLErr, (GDALDatasetH, Cstring, Cstring, Cstring, GDALDriverH, Ptr{Cstring}, GDALResampleAlg, Cdouble, Cdouble, GDALProgressFunc, Ptr{Void}, Ptr{GDALWarpOptions}), hSrcDS, pszSrcWKT, pszDstFilename, pszDstWKT, hDstDriver, papszCreateOptions, eResampleAlg, dfWarpMemoryLimit, dfMaxError, pfnProgress, pProgressArg, psOptions)
end

"""
    GDALAutoCreateWarpedVRT(GDALDatasetH hSrcDS,
                            const char * pszSrcWKT,
                            const char * pszDstWKT,
                            GDALResampleAlg eResampleAlg,
                            double dfMaxError,
                            const GDALWarpOptions * psOptions) -> GDALDatasetH

Create virtual warped dataset automatically.

### Parameters
* **hSrcDS**: The source dataset.
* **pszSrcWKT**: The coordinate system of the source image. If NULL, it will be read from the source image.
* **pszDstWKT**: The coordinate system to convert to. If NULL no change of coordinate system will take place.
* **eResampleAlg**: One of GRA_NearestNeighbour, GRA_Bilinear, GRA_Cubic, GRA_CubicSpline, GRA_Lanczos, GRA_Average or GRA_Mode. Controls the sampling method used.
* **dfMaxError**: Maximum error measured in input pixels that is allowed in approximating the transformation (0.0 for exact calculations).
* **psOptionsIn**: Additional warp options, normally NULL.

### Returns
NULL on failure, or a new virtual dataset handle on success.
"""
function GDALAutoCreateWarpedVRT(hSrcDS::GDALDatasetH, pszSrcWKT, pszDstWKT, eResampleAlg::GDALResampleAlg, dfMaxError::Cdouble, psOptions)
    ccall((:GDALAutoCreateWarpedVRT, libgdal), GDALDatasetH, (GDALDatasetH, Cstring, Cstring, GDALResampleAlg, Cdouble, Ptr{GDALWarpOptions}), hSrcDS, pszSrcWKT, pszDstWKT, eResampleAlg, dfMaxError, psOptions)
end

"""
    GDALCreateWarpedVRT(GDALDatasetH hSrcDS,
                        int nPixels,
                        int nLines,
                        double * padfGeoTransform,
                        GDALWarpOptions * psOptions) -> GDALDatasetH

Create virtual warped dataset.

### Parameters
* **hSrcDS**: The source dataset.
* **nPixels**: Width of the virtual warped dataset to create
* **nLines**: Height of the virtual warped dataset to create
* **padfGeoTransform**: Geotransform matrix of the virtual warped dataset to create
* **psOptions**: Warp options. Must be different from NULL.

### Returns
NULL on failure, or a new virtual dataset handle on success.
"""
function GDALCreateWarpedVRT(hSrcDS::GDALDatasetH, nPixels::Cint, nLines::Cint, padfGeoTransform, psOptions)
    ccall((:GDALCreateWarpedVRT, libgdal), GDALDatasetH, (GDALDatasetH, Cint, Cint, Ptr{Cdouble}, Ptr{GDALWarpOptions}), hSrcDS, nPixels, nLines, padfGeoTransform, psOptions)
end

"""
    GDALInitializeWarpedVRT(GDALDatasetH hDS,
                            GDALWarpOptions * psWO) -> CPLErr

Set warp info on virtual warped dataset.

### Parameters
* **hDS**: dataset previously created with the VRT driver, and a SUBCLASS of "VRTWarpedDataset".
* **psWO**: the warp options to apply. Note that ownership of the transformation information is taken over by the function though everything else remains the property of the caller.

### Returns
CE_None on success or CE_Failure if an error occurs.
"""
function GDALInitializeWarpedVRT(hDS::GDALDatasetH, psWO)
    ccall((:GDALInitializeWarpedVRT, libgdal), CPLErr, (GDALDatasetH, Ptr{GDALWarpOptions}), hDS, psWO)
end

"""
    GDALCreateWarpOperation(const GDALWarpOptions *) -> GDALWarpOperationH
"""
function GDALCreateWarpOperation(arg1)
    ccall((:GDALCreateWarpOperation, libgdal), GDALWarpOperationH, (Ptr{GDALWarpOptions},), arg1)
end

"""
    GDALDestroyWarpOperation(GDALWarpOperationH) -> void
"""
function GDALDestroyWarpOperation(arg1::GDALWarpOperationH)
    ccall((:GDALDestroyWarpOperation, libgdal), Void, (GDALWarpOperationH,), arg1)
end

"""
    GDALChunkAndWarpImage(GDALWarpOperationH,
                          int,
                          int,
                          int,
                          int) -> CPLErr
"""
function GDALChunkAndWarpImage(arg1::GDALWarpOperationH, arg2::Cint, arg3::Cint, arg4::Cint, arg5::Cint)
    ccall((:GDALChunkAndWarpImage, libgdal), CPLErr, (GDALWarpOperationH, Cint, Cint, Cint, Cint), arg1, arg2, arg3, arg4, arg5)
end

"""
    GDALChunkAndWarpMulti(GDALWarpOperationH,
                          int,
                          int,
                          int,
                          int) -> CPLErr
"""
function GDALChunkAndWarpMulti(arg1::GDALWarpOperationH, arg2::Cint, arg3::Cint, arg4::Cint, arg5::Cint)
    ccall((:GDALChunkAndWarpMulti, libgdal), CPLErr, (GDALWarpOperationH, Cint, Cint, Cint, Cint), arg1, arg2, arg3, arg4, arg5)
end

"""
    GDALWarpRegion(GDALWarpOperationH,
                   int,
                   int,
                   int,
                   int,
                   int,
                   int,
                   int,
                   int) -> CPLErr
"""
function GDALWarpRegion(arg1::GDALWarpOperationH, arg2::Cint, arg3::Cint, arg4::Cint, arg5::Cint, arg6::Cint, arg7::Cint, arg8::Cint, arg9::Cint)
    ccall((:GDALWarpRegion, libgdal), CPLErr, (GDALWarpOperationH, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

"""
    GDALWarpRegionToBuffer(GDALWarpOperationH,
                           int,
                           int,
                           int,
                           int,
                           void *,
                           GDALDataType,
                           int,
                           int,
                           int,
                           int) -> CPLErr
"""
function GDALWarpRegionToBuffer(arg1::GDALWarpOperationH, arg2::Cint, arg3::Cint, arg4::Cint, arg5::Cint, arg6, arg7::GDALDataType, arg8::Cint, arg9::Cint, arg10::Cint, arg11::Cint)
    ccall((:GDALWarpRegionToBuffer, libgdal), CPLErr, (GDALWarpOperationH, Cint, Cint, Cint, Cint, Ptr{Void}, GDALDataType, Cint, Cint, Cint, Cint), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11)
end

"""
    GWKGetFilterRadius(GDALResampleAlg eResampleAlg) -> int
"""
function GWKGetFilterRadius(eResampleAlg::GDALResampleAlg)
    ccall((:GWKGetFilterRadius, libgdal), Cint, (GDALResampleAlg,), eResampleAlg)
end

"""
    GWKGetFilterFunc(GDALResampleAlg eResampleAlg) -> FilterFuncType
"""
function GWKGetFilterFunc(eResampleAlg::GDALResampleAlg)
    ccall((:GWKGetFilterFunc, libgdal), FilterFuncType, (GDALResampleAlg,), eResampleAlg)
end

"""
    GWKGetFilterFunc4Values(GDALResampleAlg eResampleAlg) -> FilterFunc4ValuesType
"""
function GWKGetFilterFunc4Values(eResampleAlg::GDALResampleAlg)
    ccall((:GWKGetFilterFunc4Values, libgdal), FilterFunc4ValuesType, (GDALResampleAlg,), eResampleAlg)
end


text:
module test_25_SeveralFigures

using  ModiaPlot

include("test_base/test_25_SeveralFigures.jl")

end


text:
"""
    Categorical(p)
A *Categorical distribution* is parameterized by a probability vector `p` (of length `K`).

```math
P(X = k) = p[k]  \\quad \\text{for } k = 1, 2, \\ldots, K.
```

```julia
Categorical(p)   # Categorical distribution with probability vector p
params(d)        # Get the parameters, i.e. (p,)
probs(d)         # Get the probability vector, i.e. p
ncategories(d)   # Get the number of categories, i.e. K
```
Here, `p` must be a real vector, of which all components are nonnegative and sum to one.

**Note:** The input vector `p` is directly used as a field of the constructed distribution, without being copied.

`Categorical` is simply a type alias describing a special case of a
`DiscreteNonParametric` distribution, so non-specialized methods defined for
`DiscreteNonParametric` apply to `Categorical` as well.

External links:

* [Categorical distribution on Wikipedia](http://en.wikipedia.org/wiki/Categorical_distribution)
"""
const Categorical{P<:Real,Ps<:AbstractVector{P}} = DiscreteNonParametric{Int,P,Base.OneTo{Int},Ps}

function Categorical{P,Ps}(p::Ps; check_args=true) where {P<:Real, Ps<:AbstractVector{P}}
    check_args && @check_args(Categorical, isprobvec(p))
    return Categorical{P,Ps}(Base.OneTo(length(p)), p, check_args=check_args)
end

Categorical(p::Ps; check_args=true) where {P<:Real, Ps<:AbstractVector{P}} =
    Categorical{P,Ps}(p, check_args=check_args)

function Categorical(k::Integer; check_args=true)
    check_args && @check_args(Categorical, k >= 1)
    return Categorical{Float64,Vector{Float64}}(Base.OneTo(k), fill(1/k, k), check_args=check_args)
end

Categorical(probabilities::Real...; check_args=true) = Categorical([probabilities...]; check_args=check_args)

### Conversions

convert(::Type{Categorical{P,Ps}}, x::AbstractVector{<:Real}) where {
    P<:Real,Ps<:AbstractVector{P}} = Categorical{P,Ps}(Ps(x))

### Parameters

ncategories(d::Categorical) = support(d).stop
params(d::Categorical{P,Ps}) where {P<:Real, Ps<:AbstractVector{P}} = (probs(d),)
partype(::Categorical{T}) where {T<:Real} = T

### Statistics

function median(d::Categorical{T}) where {T<:Real}
    k = ncategories(d)
    p = probs(d)
    cp = zero(T)
    i = 0
    while cp < 1/2 && i <= k
        i += 1
        @inbounds cp += p[i]
    end
    i
end

### Evaluation

function cdf(d::Categorical{T}, x::Int) where T<:Real
    k = ncategories(d)
    p = probs(d)
    x < 1 && return zero(T)
    x >= k && return one(T)
    c = p[1]
    for i = 2:x
        @inbounds c += p[i]
    end
    return c
end

function pdf(d::Categorical, x::Real)
    ps = probs(d)
    return insupport(d, x) ? ps[round(Int, x)] : zero(eltype(ps))
end

function _pdf!(r::AbstractArray, d::Categorical{T}, rgn::UnitRange) where {T<:Real}
    vfirst = round(Int, first(rgn))
    vlast = round(Int, last(rgn))
    vl = max(vfirst, 1)
    vr = min(vlast, ncategories(d))
    p = probs(d)
    if vl > vfirst
        for i = 1:(vl - vfirst)
            r[i] = zero(T)
        end
    end
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = p[v]
    end
    if vr < vlast
        for i = (vr - vfirst + 2):length(rgn)
            r[i] = zero(T)
        end
    end
    return r
end

# sampling

sampler(d::Categorical{P,Ps}) where {P<:Real,Ps<:AbstractVector{P}} =
   AliasTable(probs(d))

### sufficient statistics

struct CategoricalStats <: SufficientStats
    h::Vector{Float64}
end

function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}) where T<:Integer
    for i = 1 : length(x)
        @inbounds xi = x[i]
        h[xi] += 1.   # cannot use @inbounds, as no guarantee that x[i] is in bound
    end
    h
end

function add_categorical_counts!(h::Vector{Float64}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    n = length(x)
    if n != length(w)
        throw(DimensionMismatch("Inconsistent array lengths."))
    end
    for i = 1 : n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        h[xi] += wi   # cannot use @inbounds, as no guarantee that x[i] is in bound
    end
    h
end

function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}) where T<:Integer
    CategoricalStats(add_categorical_counts!(zeros(k), x))
end

function suffstats(::Type{<:Categorical}, k::Int, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    CategoricalStats(add_categorical_counts!(zeros(k), x, w))
end

const CategoricalData = Tuple{Int, AbstractArray}

suffstats(::Type{<:Categorical}, data::CategoricalData) = suffstats(Categorical, data...)
suffstats(::Type{<:Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = suffstats(Categorical, data..., w)

# Model fitting

function fit_mle(::Type{<:Categorical}, ss::CategoricalStats)
    Categorical(pnormalize!(ss.h))
end

function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}) where T<:Integer
    Categorical(pnormalize!(add_categorical_counts!(zeros(k), x)), check_args=false)
end

function fit_mle(::Type{<:Categorical}, k::Integer, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    Categorical(pnormalize!(add_categorical_counts!(zeros(k), x, w)), check_args=false)
end

fit_mle(::Type{<:Categorical}, data::CategoricalData) = fit_mle(Categorical, data...)
fit_mle(::Type{<:Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = fit_mle(Categorical, data..., w)

fit_mle(::Type{<:Categorical}, x::AbstractArray{T}) where {T<:Integer} = fit_mle(Categorical, maximum(x), x)
fit_mle(::Type{<:Categorical}, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Integer} = fit_mle(Categorical, maximum(x), x, w)

fit(::Type{<:Categorical}, data::CategoricalData) = fit_mle(Categorical, data)
fit(::Type{<:Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = fit_mle(Categorical, data, w)


text:
#
# Advent of Code
# Day 2: Bathroom Security
#

println("Advent of Code")
println("Day 2: Bathroom Security")
println("")

DEBUG = false
KEYPAD_PART1 = [['0','0','0','0','0'],['0', '1', '2', '3', '0'],['0','4','5','6','0'],['0','7','8','9','0'],['0','0','0','0','0']]
KEYPAD_PART2 = [['0','0','0','0','0','0','0'],['0', '0', '0', '1', '0', '0', '0'],['0','0','2','3','4','0','0'],['0','5','6','7','8','9','0'],['0','0','A','B','C','0','0'],['0','0','0','D','0','0','0'],['0','0','0','0','0','0','0']]

function getNumbers(inputDataBits, keypad)

  numbers = Vector{Char}()
  currentLocation = [3,3] # start at centre

  for controlString in inputDataBits
    if DEBUG println("controlString: $controlString") end

    for direction in controlString
      direction = uppercase(direction)
      if DEBUG println("Direction: $direction") end

      if direction == 'U'
        if DEBUG println("Up from ", currentLocation[1]) end
        if keypad[currentLocation[1]-1][currentLocation[2]] != '0'
          currentLocation[1] -= 1
        end
      elseif direction == 'D'
        if DEBUG println("Down from ", currentLocation[1]) end
        if keypad[currentLocation[1]+1][currentLocation[2]] != '0'
          currentLocation[1] += 1
        end
      elseif direction == 'R'
        if DEBUG println("Right from ", currentLocation[2]) end
        if keypad[currentLocation[1]][currentLocation[2]+1] != '0'
          currentLocation[2] += 1
        end
      elseif direction == 'L'
        if DEBUG println("Left from ", currentLocation[2]) end
        if keypad[currentLocation[1]][currentLocation[2]-1] != '0'
          currentLocation[2] -= 1
        end
      else
        println("Unknown direction: $direction")
        quit()
      end
    end
    if DEBUG println("Current location: ", currentLocation) end
    if DEBUG println("Keypad number: ", KEYPAD[currentLocation[1]][currentLocation[2]]) end

    push!(numbers,keypad[currentLocation[1]][currentLocation[2]])
  end

  return numbers
end

# part 1 data file
# Input data has been saved to this file
inputDataFileName = "day02-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

# part 1 example
#inputData = "ULL\nRRDDD\nLURDL\nUUUUD\n"

# like python split, break the string into bits on commas
inputDataBits = split(inputData)

# in-place map function to strip the white spaces off of each inputDataBit
map!(strip,inputDataBits)

if DEBUG println("inputDataBits\n", inputDataBits) end

part1Code = getNumbers(inputDataBits, KEYPAD_PART1)
println("The part1 code is ", part1Code)

part2Code = getNumbers(inputDataBits, KEYPAD_PART2)
println("The part2 code is ", part2Code)


text:
# Carrier vector for fundamental matrix estimation.
@inline function ∂ₓu(entity::FundamentalMatrix, 𝒟)
    m, mʹ = collect(𝒟)
    𝐦  = 𝑛(m)
    𝐦ʹ = 𝑛(mʹ)
    ∂ₓu(entity, 𝐦 , 𝐦ʹ)
end

@inline function ∂ₓu(entity::FundamentalMatrix, 𝐦::AbstractVector, 𝐦ʹ::AbstractVector)
    𝐞₁ = [1.0 0.0 0.0]'
    𝐞₂ = [0.0 1.0 0.0]'
    [(𝐞₁ ⊗ 𝐦ʹ) (𝐞₂ ⊗ 𝐦ʹ) (𝐦 ⊗ 𝐞₁) (𝐦 ⊗ 𝐞₂)]
end

@inline function uₓ(entity::FundamentalMatrix, 𝒟)
    m, mʹ = collect(𝒟)
    𝐦  = 𝑛(m)
    𝐦ʹ = 𝑛(mʹ)
    uₓ(entity, 𝐦 , 𝐦ʹ)
end

@inline function uₓ(entity::FundamentalMatrix, 𝐦::AbstractVector, 𝐦ʹ::AbstractVector)
    𝐦 ⊗ 𝐦ʹ
end


text:
## for development

using Compat
using DSP
using HDF5

if !isdefined(:WarpedArray)
    include("types.jl")
end

include("rasta.jl")
include("mfccs.jl")
include("warpedarray.jl")
include("feacalc.jl")
include("io.jl")


text:
using DrWatson
@quickactivate "PolarizationFramework"

println("""
        Currently active project is: $(projectname())

        Path of active project: $(projectdir())

        Have fun with your new project!

        You can help us improve DrWatson by opening
        issues on GitHub, submitting feature requests,
        or even opening your own Pull Requests!
        """)

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
using XGPaint
using Healpix

nside = 4096
m = Map{Float64,RingOrder}(nside)
XGPaint.catalog2map!(m, [1.0], [0.2], [0.3])


text:
export
	LatLonAccelDirection,
	AccelSteeringDirection,
    NextState,
    convert2vehicle

include("actions.jl")

function convert2vehicle(veh::Records.Entity{AutomotiveDrivingModels.VehicleState,AutomotiveDrivingModels.BicycleModel,Int64})
    return Entity(veh.state,veh.def.def,veh.id)
end

function convert2vehicle(veh::Records.Entity{AutomotiveDrivingModels.VehicleState,AutomotiveDrivingModels.VehicleDef,Int64})
    return veh
end


text:
using Test
using SmartShockFinder

EXAMPLES_DIR = joinpath(pathof(SmartShockFinder) |> dirname |> dirname, "examples")

@test_nowarn rm(joinpath("datasets"), recursive=true, force=true)

@time @testset "SmartShockFinder" begin
  @time @testset "Mesh" begin
  include("test_meshes.jl")
  end

  @test_nowarn mkpath(joinpath("datasets", "1d"))

  @time @testset "1d" begin
    include("networks1d.jl")
  end

  @time @testset "2d" begin
    include("networks2d.jl")
  end
end

@test_nowarn rm(joinpath("datasets"), recursive=true, force=true)
@test_nowarn rm(joinpath("networks"), recursive=true, force=true)


text:
import MechanicalSketch
import MechanicalSketch: empty_figure, PALETTE, O, HE, finish
import MechanicalSketch: color_with_lumin, color_from_palette, set_scale_sketch
import MechanicalSketch: sethue, settext, WI, EM, circle, arrow
import MechanicalSketch: @import_expand, Quantity

let
BACKCOLOR = color_with_lumin(PALETTE[8], 90);

if !@isdefined N
    @import_expand( ~m, ~N, s)
end

include("test_functions_28.jl")

fi = "Reaction forces.txt"
lfina = joinpath(@__DIR__, "../resource/", fi)
matchlines = matchinglines(lfina, [3s], ["Vert"]);
thisline = matchlines[3]
thisname = split(thisline, "\t")[1]

Fx, Fy, Fz, _, _, _, px, py, pz = parse.(Quantity{Float64}, split(thisline, '\t')[3:end-1])

empty_figure(filename = joinpath(@__DIR__, "test_28.png"),
    backgroundcolor = BACKCOLOR,
    hue = color_from_palette("blue"));

set_scale_sketch(20m, HE)
set_scale_sketch(1500kN, HE)
sethue(color_from_palette("black"))
circle(O, 6269.0mm, :stroke)
circle(O, 5623.0mm, :stroke)

sethue(color_from_palette("green"))
drawforces(O, lfina, "xy", [3s], ["Stopper"]; components = false, labellength = true)
sethue(color_from_palette("red"))
drawforces(O, lfina, "xy", [3s], ["VertWeld"]; components = false, labellength = true)
sethue(color_from_palette("blue"))
drawforces(O, lfina, "xy", [3s], ["Deck stopper"]; components = false, labellength = true)

str = "View from above"
settext(str, O + (-WI/2 + EM, -HE / 2  + 2EM) , markup = true)

finish()
set_scale_sketch()
end


text:
export read_openeew, openeew2seisdata, OpenEEWRecord, splitjson

# structure for holding an individual OpenEEW structure
struct OpenEEWRecord
    country_code::String
    device_id::String
    x::Array{Float32,1}
    y::Array{Float32,1}
    z::Array{Float32,1}
    device_t::Float64
    cloud_t::Float64
    sr::Float64
end

function OpenEEWRecord(d::Dict)
    return OpenEEWRecord(
        d["country_code"],
        d["device_id"],
        d["x"],
        d["y"],
        d["z"],
        d["device_t"],
        d["cloud_t"],
        d["sr"],
)
end

"""
  read_openeew(filepath)

Read OpenEEW data using SeisIO.

# Arguments
- `filepath::String`: Path to OpenEEW .jsonl file.

# Returns
- `O::Array{OpenEEW}`: 3-channel SeisData structure.
"""
function read_openeew(filepath::String)
    # read JSON lines file
    records  = splitjson(read(filepath))

    # convert to Julia struct
    O = OpenEEWRecord.(JSON.parse.(records))
    return openeew2seisdata(O)
end

"""
  openeew2seisdata(O)

Convert `Array` `O` of `OpenEEWRecords` to SeisData structure.
"""
function openeew2seisdata(O::Array{OpenEEWRecord};
    vertical_axis::Union{String,Char}="x",
)
    N = length(O)
    T = eltype(vertical_axis)
    ax = T == Char ? first.(["x","y","z"]) : ["x","y","z"]
    horizontal_axes = setdiff(ax,[vertical_axis])
    vert = Symbol(vertical_axis)
    horiz1 = Symbol(horizontal_axes[1])
    horiz2 = Symbol(horizontal_axes[2])

    # allocate arrays
    x = vcat([getproperty(O[ii],horiz1) for ii = 1:N]...)
    y = vcat([getproperty(O[ii],horiz2) for ii = 1:N]...)
    z = vcat([getproperty(O[ii],vert) for ii = 1:N]...)

    # enforce constraint that channels must be same length
    Nx = length(x)
    Ny = length(y)
    Nz = length(z)
    if Nx != Ny != Nz
        throw(ArgumentError("x, y, and z channels must be same length."))
    end

    # check sample rate same for all records
    @assert length(Set(o.sr for o in O)) == 1
    sr = O[1].sr # sample rate

    # check country code same for all records
    @assert length(Set(o.country_code for o in O)) == 1
    country_code = O[1].country_code

    # check country code same for all records
    @assert length(Set(o.device_id for o in O)) == 1
    device_id = O[1].device_id

    # allocate time variable and re-sort by time
    t = vcat([sample(o.device_t,o.sr,length(o.x)) for o in O]...)
    ind = sortperm(t)
    t .= t[ind]
    x .= x[ind]
    y .= y[ind]
    z .= z[ind]

    # find gaps and interpolate data onto an even grid-spacing
    # due to time-delays, the last few points will be lost
    newt = resample_time(t,sr)
    bad = ceil(Int,(newt[end] - t[end]) * sr)
    bad = bad < 0 ? 0 : bad
    itpx = interpolate((t,), x, Gridded(Linear()))
    itpy = interpolate((t,), y, Gridded(Linear()))
    itpz = interpolate((t,), z, Gridded(Linear()))
    newx = itpx(newt[1:end-bad])
    newy = itpy(newt[1:end-bad])
    newz = itpz(newt[1:end-bad])
    newt = newt[1:end-bad]
    return openeew2seisdata(newx,newy,newz,newt,sr,country_code,device_id)
end

function sample(t::Real,sr::Real,N::Int)
    return round.(t .- (1/sr) .* collect(N-1:-1:0),digits=3)
end

function resample_time(t::AbstractArray,sr::Real)
    newt = Array{Float64}(undef,0)
    gaps = get_gaps(t,sr)
    for ii = 1:length(gaps)-1
            gapt = t[gaps[ii]+1:gaps[ii+1]]
            Ngap = length(gapt)
            resampt = round.(collect((0:Ngap-1) ./ sr .+ gapt[1]),digits=3)
            append!(newt,resampt)
    end
    return newt
end

function get_gaps(t::AbstractArray,sr::Real)
    tdiff = abs.(diff(t))
    gaps = vcat([0; findall(tdiff .>= 2 / sr )])
    append!(gaps,length(t))
    return gaps
end

"""
  openeew2seisdata(x,y,z,t,sr,country_code,device_id)

Convert `OpenEEW` data to SeisData structure.
"""
function openeew2seisdata(
    x::AbstractArray,
    y::AbstractArray,
    z::AbstractArray,
    t::AbstractArray,
    sr::Real,
    country_code::String,
    device_id::String,
)
    @assert length(x) == length(y) == length(z) == length(t)
    S = SeisData(3)

    # use SEED naming convetion
    S.id[1] = "$country_code.$device_id..ENX"
    S.id[2] = "$country_code.$device_id..ENY"
    S.id[3] = "$country_code.$device_id..ENZ"
    S.name[1] = "$country_code.$device_id..ENX"
    S.name[2] = "$country_code.$device_id..ENY"
    S.name[3] = "$country_code.$device_id..ENZ"

    # add sampling rate
    S.fs[1] = sr
    S.fs[2] = sr
    S.fs[3] = sr

    # check for gaps with t
    tmat = SeisIO.t_collapse(convert.(Int,t .* 1e6),sr)
    S.t[1] = tmat
    S.t[2] = tmat
    S.t[3] = tmat
    S.x[1] = x
    S.x[2] = y
    S.x[3] = z
    return S
end

function splitjson(A::AbstractArray)
    splits = split(String(A),"\n")

    # remove empty last record
    if splits[end] == ""
        deleteat!(splits,length(splits))
    end

    return splits
end


text:
include("../src/classification.jl")
using Main.Classification
using Random
using PyCall
using Statistics

metrics = pyimport("sklearn.metrics");
Random.seed!(5);

sklsets = pyimport("sklearn.datasets")
observables, target = sklsets.make_classification(n_samples=5000, n_features=20, n_informative=15,
n_redundant=1, n_repeated=0, n_classes=2, random_state=1, flip_y=0.05, class_sep=1.0, shuffle=true)
(Xtrain, Ytrain), (Xtest, Ytest) = traintestsplit(observables, target, 0.8, 5000);

gd = GDContext(learningrate=0.1, tolerance=1e-7, maxiterations=1000)
clf = LogisticClassifier(GradientDescent(gd))
add_designmatrix!(clf, Xtrain)
@time fit!(clf, Ytrain)
@show clf.optimizer.converged
@show mean(Ytest .== predict(clf, Xtest))


text:
# Spherical harmonic I/O, storage, and conversions

@testset "SHCilmToCindex and SHCindexToCilm" begin
    Random.seed!(0)
    # Choose random values
    lmax = 2
    cilm = randn(2, lmax + 1, lmax + 1)

    cindex = SHCilmToCindex(cilm)
    @test size(cindex) == (2, (lmax + 1) * (lmax + 2) ÷ 2)
    # These loops test all values in cindex
    for l in 0:lmax, m in 0:l, i in 1:(1 + (m > 0))
        @test cindex[i, l * (l + 1) ÷ 2 + m + 1] == cilm[i, l + 1, m + 1]
    end

    cilm′ = SHCindexToCilm(cindex)
    for l in 0:lmax, m in 0:l, i in 1:(1 + (m > 0))
        @test isapprox(cilm′[i, l + 1, m + 1], cilm[i, l + 1, m + 1])
    end

    cindex′ = SHCilmToCindex(cilm′)
    @test isapprox(cindex′, cindex)
end

@testset "SHCilmToVector SHVectorToCilm, and YlmIndexVector" begin
    Random.seed!(0)
    # Choose random values
    lmax = 2
    cilm = randn(2, lmax + 1, lmax + 1)

    vector = SHCilmToVector(cilm, lmax)
    @test length(vector) == (lmax + 1)^2
    coverage = falses(length(vector))
    for l in 0:lmax, m in 0:l, i in 1:(1 + (m > 0))
        index = YlmIndexVector(i, l, m)
        @test !coverage[index]
        coverage[index] = true
        @test vector[index] == cilm[i, l + 1, m + 1]
    end
    @test all(coverage)

    cilm′ = SHVectorToCilm(vector, lmax)
    for l in 0:lmax, m in 0:l, i in 1:(1 + (m > 0))
        @test isapprox(cilm′[i, l + 1, m + 1], cilm[i, l + 1, m + 1])
    end

    vector′ = SHCilmToVector(cilm′, lmax)
    @test isapprox(vector′, vector)
end


text:
"""
Returns path from source to sink.

It assumes that the subproblem (being passed in) has been solved to a solution
with a positive objective value. Returns an empty vector otherwise.

### Arguments
* `sp`: SubProblem whose optimal solution we're retrieving as a path.

### Returns
A vector of ints corresponding to stations along the path.
"""
function getpath(sp::SubProblem)
    if round(JuMP.getobjectivevalue(sp.model)) > 0
        visited = falses(sp.np.nstations)
        source_val = findfirst(round.(JuMP.getvalue(sp.src)) .> 0.1)
        sink_val = findfirst(round.(JuMP.getvalue(sp.snk)) .> 0.1)
        return getpath(
            source_val, source_val, sink_val, sp.edg, visited, sp.outneighbors
        )
    end

    Int[]
end

"""
Helper function for getpath(), where `cur` can be any starting node (in case of
subtours).
"""
function getpath(
        cur::Int,
        source_val::Int,
        sink_val::Int,
        edg::JuMP.JuMPDict{JuMP.Variable},
        visited::BitArray,
        outneighbors::Vector{Vector{Int}};
        maxiterations = 1000
    )
    pathnodes = Int[]
    push!(pathnodes, cur)
    visited[cur] = true
    for i in 1:maxiterations
        for nxt in outneighbors[cur]
            if round(JuMP.getvalue(edg[cur,nxt])) > 0.1
                if in(nxt, pathnodes)
                    return pathnodes
                elseif nxt == sink_val
                    push!(pathnodes, nxt)
                    visited[nxt] = true
                    return pathnodes
                end
                push!(pathnodes, nxt)
                visited[nxt] = true
                cur = nxt
                break
            end
        end
    end

    pathnodes
end


text:
@testset "Evolution generator construction" begin
  H = [1. 1. + im 3.; 1. - im 1. im; 3. -im 1.]
  L1 = sparse([1. + 0im 2. 3.; 4. 5. 6.; 6. 7. -6.])
  L2 = sparse([0. + 0im 0. 1.; 0. 0. 0.; 0. 0. 0.])
  resultnoomega = sparse([-52.0+0.0im  -29.0+1.0im    7.5+3.0im  -29.0-1.0im    4.0+0.0im     6.0+0.0im    7.5-3.0im     6.0+0.0im   10.0+0.0im;
                          -29.0+1.0im  -60.5+0.0im   10.0+0.0im    8.0+0.0im  -21.0-1.0im    12.0+0.0im   12.0+0.0im    19.5-3.0im   18.0+0.0im;
                           10.5+3.0im    9.0+0.0im  -73.5+0.0im   12.0+0.0im   14.0+0.0im   -43.0-1.0im   18.0+0.0im    21.0+0.0im  -13.5-3.0im;
                          -29.0-1.0im    8.0+0.0im   12.0+0.0im  -60.5+0.0im  -21.0+1.0im    19.5+3.0im   10.0+0.0im    12.0+0.0im   18.0+0.0im;
                           16.0+0.0im  -13.0-1.0im   24.0+0.0im  -13.0+1.0im  -53.0+0.0im    34.0+0.0im   24.0+0.0im    34.0+0.0im   36.0+0.0im;
                           24.0+0.0im   28.0+0.0im  -57.0-1.0im   34.5+3.0im   37.0+0.0im  -110.0+0.0im   36.0+0.0im    42.0+0.0im  -32.0+0.0im;
                           10.5-3.0im   12.0+0.0im   18.0+0.0im    9.0+0.0im   14.0+0.0im    21.0+0.0im  -73.5+0.0im   -43.0+1.0im  -13.5+3.0im;
                           24.0+0.0im   34.5-3.0im   36.0+0.0im   28.0+0.0im   37.0+0.0im    42.0+0.0im  -57.0+1.0im  -110.0+0.0im  -32.0+0.0im;
                           36.0+0.0im   42.0+0.0im  -31.5-3.0im   42.0+0.0im   49.0+0.0im   -40.0+0.0im  -31.5+3.0im   -40.0+0.0im  -46.0+0.0im])
  @testset "Standard usage" begin
    #no locH case

    @test evolve_generator(H, [L1, L2]) == resultnoomega
    @test evolve_generator(H, [L1, L2], 1/2) == resultnoomega/2
  end
  @testset "Nonmoralized usage" begin
    globalH = nm_glob_ham(H)
    Lnonmoral1 = nm_lind(L1)
    locH1 = nm_loc_ham(Lnonmoral1[2])
    @test evolve_generator(globalH, SparseMatrixCSC{ComplexF64}[], locH1, 1/3) ≈ evolve_generator((1-1/3)*globalH+1/3*locH1,SparseMatrixCSC{ComplexF64}[])
    @test evolve_generator(globalH, [Lnonmoral1[1]], locH1, 1/3) ≈ evolve_generator(globalH+1/2*(locH1),[Lnonmoral1[1]],1/3)

  end

  @testset "Error tests" begin
    @test_throws MethodError evolve_generator(H, [L1, L2], 1im)
    @test_throws ArgumentError evolve_generator(H, [L1, L2], -1)
    @test_throws ArgumentError evolve_generator(H, [L1, L2], 3)
  end
end

@testset "Classical lindblad generators" begin
  H = [1. 1.0+im ; 1.0-im im]
 result = [sparse([1.0+0im 0 ; 0 0 ]),
            sparse([0.0+0im 1.0+im ; 0 0 ]),
            sparse([0.0+0im 0 ; 1.0-im 0 ]),
            sparse([0.0+0im 0 ; 0 im ])]
  resultwithepsilon = [sparse([0.0+0im 1.0+im ; 0 0 ]),
           sparse([0.0+0im 0 ; 1.0-im 0 ])]
  @testset "Standard usage" begin
    @test local_lind(H) == result
    @test local_lind(sparse(H)) == result
    @test local_lind(H, epsilon = 1.1) == resultwithepsilon
    @test local_lind(sparse(H), epsilon = 1.1) == resultwithepsilon
  end

  @testset "Type tests" begin
    @test typeof(local_lind(H)) == Vector{SparseMatrixCSC{ComplexF64}}
    @test typeof(local_lind(sparse(H))) == Vector{SparseMatrixCSC{ComplexF64}}
  end

  @testset "Error tests" begin
    @test_throws ArgumentError local_lind(H, epsilon = -1.)
    @test_throws TypeError local_lind(H, epsilon = -1im)
  end
end