text:
using LinearAlgebra, Plots

#########################
##  Double Integrator  ##
#########################

model,obj = Dynamics.double_integrator

solver = Solver(model,obj,N=51)
n,m,N = get_sizes(solver)
U0 = ones(m,N-1)

solver.opts.verbose = true
res,stats = solve(solver,U0)
plot(to_array(res.U)')

opts = SolverOptions()
tol = 1e-8
opts.cost_tolerance = tol
opts.constraint_tolerance = tol
opts.cost_tolerance_intermediate = tol

# Constrained
obj_con = ConstrainedObjective(obj,u_min=-0.5,u_max=0.5)
solver = Solver(model,obj_con,N=51, opts=opts)
solver.opts.verbose = true
solver.opts.outer_loop_update_type = :default
solver.opts.penalty_update_frequency = 1
solver.opts.use_nesterov = true
res,stats = solve(solver,U0)
stats["iterations"]

solver = Solver(model,obj_con,N=51, opts=opts)
solver.state.penalty_only
solver.opts.verbose = true
solver.opts.outer_loop_update_type = :feedback
solver.opts.penalty_update_frequency = 1
solver.opts.constraint_decrease_ratio = 0.25
solver.opts.constraint_tolerance_intermediate = 1e-6
solver.opts.use_nesterov
res_m,stats_m = solve(solver,U0)
solver.state.penalty_only
stats_m["iterations"]

stats_m["cost"][end] - stats["cost"][end]

#########################
##     Pendulum        ##
#########################

model,obj = Dynamics.pendulum
costfun = copy(obj.cost)
costfun.Qf .*= 0
xf = [pi; 0]
obj_con = ConstrainedObjective(costfun,obj.tf,obj.x0,xf)
obj_con2 = Dynamics.pendulum_constrained[2]

opts = SolverOptions()
tol = 1e-8
opts.cost_tolerance = tol
opts.constraint_tolerance = tol
opts.cost_tolerance_intermediate = 1e-4
opts.constraint_tolerance_intermediate = 1e-2
opts.iterations_outerloop = 200

solver = Solver(model,obj_con,N=101,opts=opts)
solver.opts.verbose = true
n,m,N = get_sizes(solver)
U0 = ones(m,N-1)

solver.opts.outer_loop_update_type = :default
solver.opts.penalty_update_frequency = 1
solver.opts.use_nesterov = true
res,stats = solve(solver,U0)
stats["iterations"]
# no nesterov: 187
#    nesterov: 167

solver = Solver(model,obj_con,N=101,opts=opts)
TrajectoryOptimization.get_num_terminal_constraints(solver)
solver.opts.outer_loop_update_type = :feedback
solver.opts.penalty_update_frequency = 1
solver.opts.verbose = false
solver.opts.constraint_decrease_ratio = 0.25
solver.opts.use_nesterov = false
res_m,stats_m = solve(solver,U0)
stats_m["iterations"]
res_m.X[N]
# no nesterov: 176
#    nesterov: 176

stats_m["cost"][end] - stats["cost"][end]

plot(to_array(res.X)')
plot!(to_array(res_m.X)',linewidth=2)

########################
##     Dubins Car     ##
########################

model,obj = Dynamics.dubinscar
n,m = model.n, model.m
dt = 0.1

opts = SolverOptions()
opts.cost_tolerance = 1e-8
opts.constraint_tolerance = 1e-8
opts.cost_tolerance_intermediate = 1e-3

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]
obj_con = ConstrainedObjective(obj,x_min=x_min,x_max=x_max)

solver = Solver(model,obj_con,N=51,opts=opts)
n,m,N = get_sizes(solver)
U0 = ones(m,N-1)

solver.opts.outer_loop_update_type = :default
solver.opts.verbose = true
solver.opts.use_nesterov = false
res,stats = solve(solver,U0)
stats["iterations"]

solver = Solver(model,obj_con,N=51,opts=opts)
solver.state.penalty_only
solver.opts.outer_loop_update_type = :feedback
solver.opts.constraint_decrease_ratio = 0.9
solver.opts.verbose = true
solver.opts.constraint_tolerance_intermediate = 1e-4
solver.opts.use_nesterov = true
res_a,stats_a = solve(solver,U0)
stats_a["iterations"]

(stats["cost"][end] - stats_a["cost"][end])

plot(to_array(res.X)')
plot!(to_array(res_a.X)',linewidth=2)

########################
## Obstacle Avoidance ##
########################
using TrajectoryOptimization: circle_constraint, plot_obstacles, plot_trajectory!

x0 = [2.5;2.5;0.]
xf = [7.5;2.5;0.]
tf =  3.0
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

r = 0.5
s1 = 30; s2 = 50; s3 = 15
circles = []
for i in range(0,stop=5,length=s1)
    push!(circles,(0.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(5.,i,r))
end
for i in range(0,stop=5,length=s1)
    push!(circles,(10.,i,r))
end
for i in range(0,stop=10,length=s2)
    push!(circles,(i,0.,r))
end
for i in range(0,stop=3,length=s3)
    push!(circles,(i,5.,r))
end
for i in range(5,stop=8,length=s3)
    push!(circles,(i,5.,r))
end

n_circles = 3*s1 + s2 + 2*s3

# plt = plot(title="Escape",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
# plot_obstacles(circles,:red)
# display(plt)

function cI(c,x,u)
    for i = 1:n_circles
        c[i] = circle_constraint(x,circles[i][1],circles[i][2],circles[i][3])
    end
    c
end

obj_con_obstacles = TrajectoryOptimization.ConstrainedObjective(obj,cI=cI)

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-6
opts.R_infeasible = 1.0

solver = Solver(model, obj_con_obstacles, dt=dt, opts=opts)
n,m,N = get_sizes(solver)

# -Initial state and control trajectories
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = TrajectoryOptimization.interp_rows(N,tf,Array(X_guess'))
# X0 = line_trajectory(solver_uncon)
U0 = rand(m,N)

solver = Solver(model, obj_con_obstacles, dt=dt, opts=opts)
solver.opts.use_nesterov = true
solver.opts.outer_loop_update_type = :default
res, stats = TrajectoryOptimization.solve(solver,X0,U0)
stats["iterations"]
# w/o nesterov: 113
# w/  nesterov:  87

solver = Solver(model, obj_con_obstacles, dt=dt, opts=opts)
solver.opts.outer_loop_update_type = :feedback
solver.opts.use_nesterov = false
solver.opts.constraint_tolerance_intermediate = 1e-4
res_a, stats_a = solve(solver,X0,U0)
stats_a["iterations"]
# w/o nesterov: 128
# w/  nesterov: 107

# w/o nesterov: 146
# w/  nesterov:  82

stats["cost"][end] - stats_a["cost"][end]

# Escape
plt = plot(title="Escape",aspect_ratio=:equal)
plot_obstacles(circles)
plot_trajectory!(to_array(res.X),width=2,color=:purple,label="Default",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot_trajectory!(to_array(res_a.X),width=2,color=:black,linestyle=:dash,label="Accelerated",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
plot!(X0[1,:],X0[2,:],label="Infeasible Initialization",width=1,color=:black,linestyle=:dash)
display(plt)


text:
#--------------------------------------------------------------------#
#     Input data used for heat and moisture transfer calculation     #
#--------------------------------------------------------------------#

#DODATI RUBNE UVJETE

#########
#Timestep if necessary#
##########

#dt=1        #Timestep in h
#T=8760      #Time domain

###########
#Mesh size#
###########

dl=0.005 #defined in m

########################
#Geometry and materials#
########################

#Elements -- matrix of vectors E_1=[material_key,x_0,y_0,x_1,y_1] to E_n

ne=5 #number of elements

E=zeros(Float64,ne,9) #Geometry and material matrix + boundary conditions

#Element 1 (Concrete)
E[1,1]=12
E[1,2]=0
E[1,3]=0
E[1,4]=0.4
E[1,5]=0.8
E[1,6]=1
E[1,7]=1
E[1,8]=1
E[1,9]=1

##############################
#Boundary condition materials#
##############################

#Element 17 (dolje)
E[2,1]=100
E[2,2]=0
E[2,3]=-dl
E[2,4]=0.4
E[2,5]=0
E[2,6]=17
E[2,7]=(10000)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[2,8]=0
E[2,9]=1

#Element 18 (lijevo)
E[3,1]=101
E[3,2]=-dl
E[3,3]=0
E[3,4]=0
E[3,5]=0.8
E[3,6]=18
E[3,7]=(10000)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[3,8]=0
E[3,9]=1

#Element 19 (gore)
E[4,1]=102
E[4,2]=0
E[4,3]=0.8
E[4,4]=0.4
E[4,5]=0.8+dl
E[4,6]=19
E[4,7]=(10000)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[4,8]=20
E[4,9]=1

#Element 20 (desno)
E[5,1]=103
E[5,2]=0.4
E[5,3]=0
E[5,4]=0.4+dl
E[5,5]=0.8
E[5,6]=20
E[5,7]=(10000)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[5,8]=20
E[5,9]=-1

ext: jl


text:
#
# Denoise an SPD Example with Douglas Rachford
#
# L2-TV functional with anisotropic TV
#
#
#
using Manopt
using Images, CSV, DataFrames, LinearAlgebra, JLD2, Dates
#
# Settings
ExportOrig = false
ExportResult = true
ExportTable = true
resultsFolder = "src/examples/Total_Variation/SPD_TV/"
comparisonData = "ImageCPPA-CostValue.jld2"
experimantName = "ImageDR"
if !isdir(resultsFolder)
    mkdir(resultsFolder)
end
#
# Manifold & Data
f = artificialSPDImage2(32)
if ExportOrig
  renderAsymptote(resultsFolder*experimantName*"orig.asy", asyExportSPDData; data=f, scaleAxes=(7.5,7.5,7.5))
end
#
# Parameters
η  = 0.58
λ = 0.93
α = 6.
#
# Build Problem for L2-TV
pixelM = SymmetricPositiveDefinite(3);
M = Power(pixelM,size(f))
d = length(size(f))
rep(d) = (d>1) ? [ones(Int,d)...,d] : d
fidelity(x) = 1/2*distance(M,x,f)^2
Λ(x) = forwardLogs(M,x) # on T_xN
prior(x) = norm(norm.(Ref(pixelM),getValue(repeat(x,rep(d)...)),getValue(Λ(x))),1)
#
# Setup & Optimize
print("--- Douglas–Rachford with η: ",η," and λ: ",λ," ---\n")
cost(x) = fidelity(x) + α*prior(x)
prox1 = (η,x) -> PowPoint(cat( proxDistance(M,η,f,x[1]), proxParallelTV(M,α*η,getValue(x[2:5])), dims=1))
prox2 = (η,x) -> PowPoint(fill(mean(M,getValue(x);stoppingCriterion=stopAfterIteration(20)),5))
sC = stopAfterIteration(400)
try
  costFunctionThreshold = load(resultsFolder*comparisonData)["compareCostFunctionValue"]
  global sC = stopWhenCostLess(costFunctionThreshold)
  @info "Comparison to CPPA (`SPDImage_CPPA.jl`) and its cost of $costFunctionThreshold."
catch y
  if isa(y, SystemError)
    @info "Comparison to CPPA only possible after runninng `SPDImage_CPPA.jl` its cost was stored."
  end
end
x0 = f
@time y, r = DouglasRachford(M, cost, [prox1,prox2], f;
  λ = i -> η, α = i -> λ, # map from Paper notation of BPS16 to toolbox notation
  debug = [:Iteration," | ", :Change, " | ", :Cost,"\n",10,:Stop],
  record = [:Iteration, :Cost ],
  stoppingCriterion = sC,
  parallel=5
)
#
# Result
numIter = length(r)
if ExportResult
  renderAsymptote(resultsFolder*experimantName*"img-result-$(numIter)-α$(replace(string(α), "." => "-")).asy", asyExportSPDData; data=y, render=4, scaleAxes=(7.5,7.5,7.5) )
end
if ExportTable
  A = cat( [ri[1] for ri in r], [ri[2] for ri in r]; dims=2 )
  CSV.write(resultsFolder*experimantName*"-Cost.csv",  DataFrame(A), writeheader=false);
end


text:
using PyCall
using Plots
using DataFrames
gr()

pth = "/home/castaned/Documents/ROTORCmodels/viscalc_suite/"
unshift!(PyVector(pyimport("sys")["path"]), pth)
@pyimport main_modules as mmod
@pyimport legendre_interp as lint

vels = 1:4
zps = Array{Float64}[]
mselect = "1p875"
ells = [6]

for vel in vels
  for m in ells
    amode = mmod.emode(mselect,Int(vel),"6 g4");
    xi_r,xi_t,dt_t,zg,r,zp,sig,cs = mmod.load_MODE_file([amode[:mass]],amode[:vel_idx],amode[:index],amode[:parity],amode[:frequency],true,false);
    if amode[:parity] == "EVEN"
          zp_fine = lint.leg_interp(zp[end-10:end,:],8,"EVEN")
      else
          zp_fine = lint.leg_interp(zp[end-10:end,:],8,"OE")
    end
    push!(zps,zp_fine[1,:])
  end
end

df = DataFrame()
for i in 1:length(zps)
  df[string("v",i)] = zps[i]
end
#using Plots
#gr()
# plotlyjs()
# pl1 = plot(linspace(0,90,length(zps[1][:])),zps[1][:],xlabel="R[zone]",ylabel="zp",legend=false,linewidth=1.5)
# plot!(linspace(0,90,length(zps[5][:])),zps[5][:],xlabel="R[zone]",ylabel="zp",legend=false,linewidth=1.5)
# plot!(linspace(0,90,length(zps[3][:])),zps[3][:],xlabel="R[zone]",ylabel="zp",legend=false,linewidth=1.5)
# PlotlyJS.relayout!(pl1.o, margin = Dict("l"=>100, "b"=>80, "r"=>8, "t"=>20))

# using PGFPlots
# xvec = 1:length(xi_r[:,1])
# plts = [Plots.Linear(1:length(xi_r[:,i]),xi_r[:,i],mark="none")::PGFPlots.Plots.Linear for i in [1,4,8]]
# p =Axis(plts,title = L"M=1.875M$\odot$ V=0 km s$^{-1}$",ylabel=L"$\xi_{r}$",xlabel="R[zone]")
# save("myfile.pdf", p)

# using Gadfly
# pl1=Gadfly.plot([layer(y=xi_r[:,i],x=1:length(xi_r[:,i]), Geom.line,Theme(line_width=0.6mm,default_color=Gadfly.distinguishable_colors(8)[i])) for i in [1,4,8]]...,Guide.xlabel("R[zone]"),Guide.ylabel("ξ_r"))
# draw(PDF("gadflyplot.pdf", 10inch, 6inch), pl1)

max_stars_repo_name: dfcastap/visuite

max_issues_repo_name: dfcastap/visuite

max_forks_repo_name: dfcastap/visuite


text:
using Multisets, IterTools

function print_map_est(output::SisPosteriorMcmcOutput; top_num::Int=5)
    d = proportionmap(output.S_sample)
    # for key in keys(d)
    #     d[key] /= length(output.S_sample)
    # end
    props = sort(collect(d), by=x->x[2], rev=true)
    title = "\nPosterior probability of mode"
    println(title)
    println("-"^length(title), "\n")
    for i in 1:min(top_num, length(d))
        println(props[i][2],"  ", props[i][1])
    end
    println("\n...showing top $(min(top_num, length(d))) interaction sequences.")
end

function print_map_est(sample::InteractionSequenceSample; top_num::Int=5)
    d = proportionmap(sample)
    # for key in keys(d)
    #     d[key] /= length(output.S_sample)
    # end
    props = sort(collect(d), by=x->x[2], rev=true)
    title = "\nPosterior probability of mode"
    println(title)
    println("-"^length(title), "\n")
    for i in 1:min(top_num, length(d))
        println(props[i][2],"  ", props[i][1])
    end
    println("\n...showing top $(min(top_num, length(d))) interaction sequences.")
end

function mean_dists_summary(
    samples::Vector{InteractionSequenceSample{T}},
    d::Union{Metric,Metric}
    ) where {T<:Union{Int,String}}

    tot_var = sample_frechet_var(
            vcat(samples...),
            d,
            with_memory=true, show_progress=true
        )
    mean_of_var = mean(
        x->sample_frechet_var(
            x,
            d,
            with_memory=true, show_progress=true), samples
    )
    means = map(x->sample_frechet_mean(x, d, with_memory=true)[1], samples)
    var_of_mean = sample_frechet_var(means, d)

    out = Dict{String,Float64}(
        "tot_var" => tot_var,
        "mean_of_vars" => mean_of_var,
        "var_of_means" => var_of_mean
    )
    return out
end

function mean_dists_summary(
    chains::Vector{SisPosteriorMcmcOutput}
    )
    samples = [x.S_sample for x in chains]
    d = chains[1].posterior.dist
    return mean_dists_summary(samples, d)
end

function print_map_est(output::Union{SimPosteriorMcmcOutput,SimPosteriorModeConditionalMcmcOutput}; top_num::Int=5)
    multiset_sample = Multiset.(output.S_sample)
    d = proportionmap(multiset_sample)
    # for key in keys(d)
    #     d[key] /= length(output.S_sample)
    # end
    props = sort(collect(d), by=x->x[2], rev=true)
    title = "\nPosterior probability of mode"
    println(title)
    println("-"^length(title), "\n")
    for i in 1:min(top_num, length(d))
        println(props[i][2],"  ", props[i][1])
    end
    println("\n...showing top $(min(top_num, length(d))) interaction multisets.")
end


text:
###############################################################################
#
# Leaf fluorescence-related parameter set
#
###############################################################################
#= AbstractFluoModelParaSet type tree
---> FluoParaSet
=#
"""
    abstract type AbstractFluoModelParaSet{FT}

Hierarchy of the `AbstractFluoModelParaSet`:
- [`FluoParaSet`](@ref)
"""
abstract type AbstractFluoModelParaSet{FT} end

"""
    mutable struct CytoFluoParaSet{FT}

A `AbstractFluoModelParaSet` type paramter set.

# Fields
$(TYPEDFIELDS)
"""
struct CytoFluoParaSet{FT<:AbstractFloat} <: AbstractFluoModelParaSet{FT}
    "Fluorescence model coefficient"
    Kr1::FT
    "Fluorescence model coefficient"
    Kr2::FT
    "Fluorescence model coefficient"
    Kr3::FT
end

"""
    mutable struct FluoParaSet{FT}

A `AbstractFluoModelParaSet` type paramter set.

# Fields
$(TYPEDFIELDS)
"""
struct FluoParaSet{FT<:AbstractFloat} <: AbstractFluoModelParaSet{FT}
    "Fluorescence model coefficient"
    Kr1::FT
    "Fluorescence model coefficient"
    Kr2::FT
    "Fluorescence model coefficient"
    Kr3::FT
end

max_stars_repo_name: climate-machine/Land


text:
using LikelihoodfreeInference
using Test

include("kernels.jl")


text:
convert(::Type{Chains}, x::T) where T <: Fit = get_mcmc_chains(x; include_subject_specific=false)

## convert Fit object to MCMCChains.Chains object
function get_mcmc_chains(x::T; include_subject_specific=false) where T <: Fit
    get_mcmc_chains(x.chains; include_subject_specific=include_subject_specific)
end

function get_mcmc_chains(x::T; include_subject_specific=false) where T <: Chain
    get_mcmc_chains([x]; include_subject_specific=include_subject_specific)
end

## convert vector of Chain objects to MCMCChains.Chains object
function get_mcmc_chains(x::AbstractVector{T}; include_subject_specific=false) where T <: Chain
    K = length(x)
    M = length(first(x))

    x0 = first(x).init
    N = length(vectorize(x0; include_subject_specific=include_subject_specific))
    ch = Array{Float64,3}(undef, M, N, K)

    for k in 1:K
        s = get_samples(x[k])
        for m in 1:M
            ch[m,:,k] = vectorize(s[m]; include_subject_specific=include_subject_specific)
        end
    end
    Chains(ch, variables(x0; include_subject_specific=include_subject_specific))
end


text:
using LazySets.Arrays: _rationalize

for N in [Float32, Float64]
    x = [N(1), N(2), N(3)]
    out = [1//1, 2//1, 3//1]
    @test _rationalize(x) == out
    v = _rationalize(BigInt, x)
    @test v[1] isa Rational{BigInt}
    @test _rationalize(x, tol=2*eps(N)) == out
end


text:
using Polyhedra

using SparseArrays
using LinearAlgebra
using Test

using StaticArrays

include("utils.jl")

include("solvers.jl")

include("elements.jl")
include("comp.jl")
include("representation.jl")
include("polyhedron.jl")

include("redundancy.jl")
include("doubledescription.jl")

include("interval.jl")

include("polyhedra_to_lpqp.jl")
include("default.jl")

include("show.jl")

include("polyhedra.jl")
for (arith, T) in (("floating point", Float64), ("exact", Rational{BigInt}))
    @testset "Polyhedra tests in $arith arithmetic" begin
        polyhedratest(Polyhedra.DefaultLibrary{T}(lp_solver), ["board"])
    end
end


text:
using DimensionalData, Aqua, SafeTestsets

if VERSION >= v"1.5.0"
    # This is catching some unambiguous constructors for T<:Metadata.
    # Aqua.test_ambiguities([DimensionalData, Base, Core])
    Aqua.test_unbound_args(DimensionalData)
    Aqua.test_undefined_exports(DimensionalData)
    Aqua.test_project_extras(DimensionalData)
    Aqua.test_stale_deps(DimensionalData)
    Aqua.test_deps_compat(DimensionalData)
    Aqua.test_project_toml_formatting(DimensionalData)
    Aqua.test_project_extras(DimensionalData)
    Aqua.test_stale_deps(DimensionalData)
end

# For coverage stats only
DimensionalData._precompile()

@time @safetestset "interface" begin include("interface.jl") end
@time @safetestset "metadata" begin include("metadata.jl") end
@time @safetestset "name" begin include("name.jl") end

@time @safetestset "lookup" begin include("lookup.jl") end
@time @safetestset "selector" begin include("selector.jl") end

@time @safetestset "coord" begin include("coord.jl") end
@time @safetestset "dimension" begin include("dimension.jl") end
@time @safetestset "primitives" begin include("primitives.jl") end
@time @safetestset "format" begin include("format.jl") end

@time @safetestset "array" begin include("array.jl") end
@time @safetestset "stack" begin include("stack.jl") end
@time @safetestset "indexing" begin include("indexing.jl") end
@time @safetestset "methods" begin include("methods.jl") end
@time @safetestset "broadcast" begin include("broadcast.jl") end
@time @safetestset "matmul" begin include("matmul.jl") end

@time @safetestset "dimindices" begin include("dimindices.jl") end
@time @safetestset "adapt" begin include("adapt.jl") end
@time @safetestset "set" begin include("set.jl") end
@time @safetestset "show" begin include("show.jl") end
@time @safetestset "tables" begin include("tables.jl") end
@time @safetestset "utils" begin include("utils.jl") end

if Sys.islinux()
    # Unfortunately this can hang on other platforms.
    # Maybe ram use of all the plots on the small CI machine? idk
    @time @safetestset "plotrecipes" begin include("plotrecipes.jl") end
end


text:
# in module DTW...

"""
A Component contains two `Int64` and a weight `Float64`. These represent the
quantity subtracted from the position (i, j), e.g.:

    `weight * d[i - row, j - col]`
"""
struct Component
    row::Int64 # row index difference with the current position
    col::Int64 # col index difference with the current position
    weight::Float64
end

function _consistency_check(components::Vector{Component})::Vector{Component}
    # last component must be [0, 0]
    if components[end].row != 0 || components[end].row != 0
        error("Last components of rules must be [0, 0]")
    end

    # components in decreasing rows and cols
    prev_row = components[1].row
    prev_col = components[1].col
    @views for comp in components[2:end]
        if comp.row > prev_row
            error("Rules must have rows in decreasing order")
        end
        if comp.col > prev_col
            error("Rules must have cols in decreasing order")
        end
    end

    return components
end

"""
A Rule contains a `Vector{Component}` so that:

1. the last component should always be 0, 0
2. both rows and cols should be in decreasing order

The components are summed so that a rule is created
"""
struct Rule
    # the steps starting from the end
    components::Vector{Component}
    weight::Float64
    bias::Float64
    # a consistency check for rules
    Rule(x) = new(_consistency_check(x), 1.0, 0)
    Rule(x; weight=1.0, bias=0.0) = new(_consistency_check(x), weight, bias)
    Rule(x, weight, bias) = new(_consistency_check(x), weight, bias)
end

"""
A StepPattern contains a `Vector{Rule}` sorted from the most important to the least important

You can print it by using `recursion_rule`
"""
struct StepPattern
    rules::Vector{Rule}
    norm_hint::String
    weight::Float64

    # a constructor with no normalization hint and no fede weight
    StepPattern(x) = new(x, "NA", 1.0)

    # a constructor for dtw-python objects with normalization hint
    StepPattern(x, y::String) = new(x, y, 1.0)

    # a constructor for fede objects with weight
    StepPattern(x, y::Float64) = new(x, "NA", y)

    # a constructor for fede objects with weight and normalization hint
    StepPattern(x, y, z) = new(x, y, z)
end

"""
Returns a new same step-pattern with row and cols inverted (transposed)
"""
function tr(x::StepPattern)::StepPattern
    rules = Vector{Rule}(undef, length(x.rules))
    for (r, rule) in enumerate(x.rules)
        components = Vector{Component}(undef, length(rule.components))
        for (c, comp) in enumerate(rule.components)
            components[c] = Component(comp.col, comp.row, comp.weight)
        end
        rules[r] = Rule(components)
    end
    return StepPattern(rules)
end

"""
A function to convert dtw-python step-patterns
"""
function _c(args...)
    num_rule = 0
    row = col = weight = 0
    rules = Rule[]
    components = Component[]
    for (i, arg) in enumerate(args)
        el = i % 4
        if el == 1
            # a new component
            # push back the previous component if any
            if num_rule > 0
                push!(components, Component(row, col, weight))
            end

            if arg != num_rule
                # a new rule
                # push back the previous rule if any
                if num_rule > 0
                    push!(rules, Rule(components))
                end

                # prepare a new vector of components for this new rule
                components = Component[]
                num_rule = arg
            end
        elseif el == 2
            # the row of this component
            row = arg
        elseif el == 3
            # the col of this component
            col = arg
        elseif el == 0
            # the weight of this component
            weight = abs(arg)
        end
    end
    # push back last component
    push!(components, Component(row, col, weight))
    # push back last rule
    push!(rules, Rule(components))
    return rules
end

#! format: off
"""
an asymmetric pattern which favours the vertical paths (changing row is
easier than changing column)
"""
const fede_asymmetric = StepPattern(
                             [
                              Rule([
                                    Component(1, 1, 1),
                                    Component(0, 0, 3),
                                   ]),

                              Rule([
                                    Component(1, 0, 1),
                                    Component(0, 0, 1),
                                   ]),

                              Rule([
                                    Component(0, 1, 1),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(1, 0, 2),
                                    Component(0, 0, 1),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(1, 1, 2),
                                    Component(0, 0, 1),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(0, 0, 3),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(0, 1, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(1, 1, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(0, 0, 4),
                                   ]),
                             ]
                            )

"""
a symmetric pattern for DTW
"""
const fede_symmetric = StepPattern(
                        [
                              Rule([
                                    Component(1, 1, 1),
                                    Component(0, 0, 3),
                                   ]),

                              Rule([
                                    Component(1, 0, 1),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(0, 1, 1),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(1, 0, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(1, 1, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(2, 1, 1),
                                    Component(0, 0, 4),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(0, 1, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(1, 1, 2),
                                    Component(0, 0, 2),
                                   ]),

                              Rule([
                                    Component(1, 2, 1),
                                    Component(0, 0, 4),
                                   ]),
                        ]
                       )

"""
Given a vector of weights creates a step-pattern using standard rules
"""
function weight_step_pattern(w::AbstractVector{Float64})::StepPattern
    return StepPattern([
        Rule([
                Component(1, 1, w[1]),
                Component(0, 0, w[2])
            ],
            w[3], w[4] # weight and bias
        ),

        Rule([
                Component(1, 0, w[5]),
                Component(0, 0, w[6])
            ],
            w[7], w[8] # weight and bias
        ),

        Rule([
                Component(0, 1, w[9]),
                Component(0, 0, w[10])
            ],
            w[11], w[12] # weight and bias
        )],
        w[13]
      )

    # return StepPattern([
    #     Rule([
    #           Component(1, 1, w[1]),
    #           Component(0, 0, w[2]),
    #           ],
    #          1, 0),
    #          # w[3], w[4]),

    #     Rule([
    #           Component(1, 0, w[3]),
    #           Component(0, 0, w[4]),
    #           ],
    #          1, 0),
    #          # w[7], w[8],

    #     Rule([
    #           Component(0, 1, w[5]),
    #           Component(0, 0, w[6]),
    #           ],
    #          1, 0),
    #          # w[11], w[12]),

    #     Rule([
    #           Component(2, 1, w[7]),
    #           Component(1, 0, w[8]),
    #           Component(0, 0, w[9]),
    #           ],
    #          1, 0),
    #          # w[16], w[17]),

    #     Rule([
    #           Component(2, 1, w[10]),
    #           Component(1, 1, w[11]),
    #           Component(0, 0, w[12]),
    #           ],
    #          1, 0),
    #          # w[21], w[22]),

    #     Rule([
    #           Component(2, 1, w[13]),
    #           Component(0, 0, w[14]),
    #           ],
    #          1, 0),
    #          # w[25], w[26]),

    #     Rule([
    #           Component(1, 2, w[15]),
    #           Component(0, 1, w[16]),
    #           Component(0, 0, w[17]),
    #           ],
    #          1, 0),
    #          # w[30], w[31]),

    #     Rule([
    #           Component(1, 2, w[18]),
    #           Component(1, 1, w[19]),
    #           Component(0, 0, w[20]),
    #           ],
    #          1, 0),
    #          # w[35], w[36]),

    #     Rule([
    #           Component(1, 2, w[21]),
    #           Component(0, 0, w[22]),
    #          ],
    #          1, 0),
    #          # w[39], w[40]),
    #    ],
    #    w[23])
end

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

## Everything here is semi auto-generated from the python source, which is
## auto-generated from the R source. Don't edit!

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

##
## Various step patterns, defined as internal variables
##
## First column: enumerates step patterns.
## Second   	 step in query index
## Third	 step in reference index
## Fourth	 weight if positive, or -1 if starting point
##
## For cite{} see dtw.bib in the package
##

## Widely-known variants

"""
White-Neely symmetric (default)
aka Quasi-symmetric cite{White1976}
normalization: no (N+M?)
"""
const symmetric1 = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 1,
    2, 0, 1, -1,
    2, 0, 0, 1,
    3, 1, 0, -1,
    3, 0, 0, 1
))

"""
Normal symmetric
normalization: N+M
"""
const symmetric2 = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 2,
    2, 0, 1, -1,
    2, 0, 0, 1,
    3, 1, 0, -1,
    3, 0, 0, 1
), "N+M")

"""
classic asymmetric pattern: max slope 2, min slope 0
normalization: N
"""
const asymmetric = StepPattern(_c(
    1, 1, 0, -1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 0, 1
), "N")

# % item{code{symmetricVelichkoZagoruyko}}{symmetric, reproduced from %
# [Sakoe1978]. Use distance matrix code{1-d}}
#

"""
normalization: max[N,M]
note: local distance matrix is 1-d
cite{Velichko}
"""
const _symmetricVelichkoZagoruyko = StepPattern(_c(
    1, 0, 1, -1,
    1, 0, 0, 0,
    2, 1, 1, -1,
    2, 0, 0, -1.001,
    3, 1, 0, -1,
    3, 0, 0, 0,
    ))

# % item{code{asymmetricItakura}}{asymmetric, slope contrained 0.5 -- 2
# from reference [Itakura1975]. This is the recursive definition % that
# generates the Itakura parallelogram }
#

"""
Itakura slope-limited asymmetric cite{Itakura1975}
Max slope: 2 min slope: 1/2
normalization: N
"""
const _asymmetricItakura = StepPattern(_c(
    1, 1, 2, -1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 2, 1, -1,
    3, 1, 0, 1,
    3, 0, 0, 1,
    4, 2, 2, -1,
    4, 1, 0, 1,
    4, 0, 0, 1
))

#############################
## Slope-limited versions
##
## Taken from Table I, page 47 of "Dynamic programming algorithm
## optimization for spoken word recognition," Acoustics, Speech, and
## Signal Processing, vol.26, no.1, pp. 43-49, Feb 1978 URL:
## http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1163055
##
## Mostly unchecked

"""
Row P=0
"""
symmetricP0 = symmetric2

"""
normalization: N ?
"""
const asymmetricP0 = StepPattern(_c(
    1, 0, 1, -1,
    1, 0, 0, 0,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 0, -1,
    3, 0, 0, 1
), "N")

"""
alternative implementation
"""
const _asymmetricP0b = StepPattern(_c(
    1, 0, 1, -1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 0, -1,
    3, 0, 0, 1
), "N")

"""
Row P=1/2
"""
const symmetricP05 = StepPattern(_c(
    1, 1, 3, -1,
    1, 0, 2, 2,
    1, 0, 1, 1,
    1, 0, 0, 1,
    2, 1, 2, -1,
    2, 0, 1, 2,
    2, 0, 0, 1,
    3, 1, 1, -1,
    3, 0, 0, 2,
    4, 2, 1, -1,
    4, 1, 0, 2,
    4, 0, 0, 1,
    5, 3, 1, -1,
    5, 2, 0, 2,
    5, 1, 0, 1,
    5, 0, 0, 1
), "N+M")

const asymmetricP05 = StepPattern(_c(
    1, 1, 3, -1,
    1, 0, 2, 1 / 3,
    1, 0, 1, 1 / 3,
    1, 0, 0, 1 / 3,
    2, 1, 2, -1,
    2, 0, 1, .5,
    2, 0, 0, .5,
    3, 1, 1, -1,
    3, 0, 0, 1,
    4, 2, 1, -1,
    4, 1, 0, 1,
    4, 0, 0, 1,
    5, 3, 1, -1,
    5, 2, 0, 1,
    5, 1, 0, 1,
    5, 0, 0, 1
), "N")

"""
Row P=1
Implementation of Sakoe's P=1, Symmetric algorithm
"""
const symmetricP1 = StepPattern(_c(
    1, 1, 2, -1,  # First branch: g(i-1,j-2)+
    1, 0, 1, 2,  # + 2d(i  ,j-1)
    1, 0, 0, 1,  # +  d(i  ,j)
    2, 1, 1, -1,  # Second branch: g(i-1,j-1)+
    2, 0, 0, 2,  # +2d(i,  j)
    3, 2, 1, -1,  # Third branch: g(i-2,j-1)+
    3, 1, 0, 2,  # + 2d(i-1,j)
    3, 0, 0, 1  # +  d(  i,j)
), "N+M")

const asymmetricP1 = StepPattern(_c(
    1, 1, 2, -1,
    1, 0, 1, .5,
    1, 0, 0, .5,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 2, 1, -1,
    3, 1, 0, 1,
    3, 0, 0, 1
), "N")

"""
Row P=2
"""
const symmetricP2 = StepPattern(_c(
    1, 2, 3, -1,
    1, 1, 2, 2,
    1, 0, 1, 2,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 2,
    3, 3, 2, -1,
    3, 2, 1, 2,
    3, 1, 0, 2,
    3, 0, 0, 1
), "N+M")

const asymmetricP2 = StepPattern(_c(
    1, 2, 3, -1,
    1, 1, 2, 2 / 3,
    1, 0, 1, 2 / 3,
    1, 0, 0, 2 / 3,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 3, 2, -1,
    3, 2, 1, 1,
    3, 1, 0, 1,
    3, 0, 0, 1
), "N")

################################
## Taken from Table III, page 49.
## Four varieties of DP-algorithm compared

## 1st row:  asymmetric

## 2nd row:  symmetricVelichkoZagoruyko

## 3rd row:  symmetric1

## 4th row:  asymmetricItakura

#############################
## Classified according to Rabiner
##
## Taken from chapter 2, Myers' thesis [4]. Letter is
## the weighting function:
##
##      rule       norm   unbiased
##   a  min step   ~N     NO
##   b  max step   ~N     NO
##   c  x step     N      YES
##   d  x+y step   N+M    YES
##
## Mostly unchecked

# R-Myers     R-Juang
# type I      type II
# type II     type III
# type III    type IV
# type IV     type VII

const typeIa = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1,
    1, 0, 0, 0,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, 1,
    3, 0, 0, 0
))

const typeIb = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, 1,
    3, 0, 0, 1
))

const typeIc = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, 1,
    3, 0, 0, 0
), "N")

const typeId = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 2,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 2,
    3, 1, 2, -1,
    3, 0, 1, 2,
    3, 0, 0, 1
), "N+M")

## ----------
## smoothed variants of above

const typeIas = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, .5,
    1, 0, 0, .5,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, .5,
    3, 0, 0, .5
))

const typeIbs = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, 1,
    3, 0, 0, 1
))

const typeIcs = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 1, 2, -1,
    3, 0, 1, .5,
    3, 0, 0, .5
), "N")

const typeIds = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 1.5,
    1, 0, 0, 1.5,
    2, 1, 1, -1,
    2, 0, 0, 2,
    3, 1, 2, -1,
    3, 0, 1, 1.5,
    3, 0, 0, 1.5
), "N+M")

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

const typeIIa = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 1,
    2, 1, 2, -1,
    2, 0, 0, 1,
    3, 2, 1, -1,
    3, 0, 0, 1
))

const typeIIb = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 1,
    2, 1, 2, -1,
    2, 0, 0, 2,
    3, 2, 1, -1,
    3, 0, 0, 2
))

const typeIIc = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 1,
    2, 1, 2, -1,
    2, 0, 0, 1,
    3, 2, 1, -1,
    3, 0, 0, 2
), "N")

const typeIId = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 2,
    2, 1, 2, -1,
    2, 0, 0, 3,
    3, 2, 1, -1,
    3, 0, 0, 3
), "N+M")

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

## Rabiner [3] discusses why this is not equivalent to Itakura's

const typeIIIc = StepPattern(_c(
    1, 1, 2, -1,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 1,
    3, 2, 1, -1,
    3, 1, 0, 1,
    3, 0, 0, 1,
    4, 2, 2, -1,
    4, 1, 0, 1,
    4, 0, 0, 1
), "N")

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

## numbers follow as production rules in fig 2.16

const typeIVc = StepPattern(_c(
    1, 1, 1, -1,
    1, 0, 0, 1,
    2, 1, 2, -1,
    2, 0, 0, 1,
    3, 1, 3, -1,
    3, 0, 0, 1,
    4, 2, 1, -1,
    4, 1, 0, 1,
    4, 0, 0, 1,
    5, 2, 2, -1,
    5, 1, 0, 1,
    5, 0, 0, 1,
    6, 2, 3, -1,
    6, 1, 0, 1,
    6, 0, 0, 1,
    7, 3, 1, -1,
    7, 2, 0, 1,
    7, 1, 0, 1,
    7, 0, 0, 1,
    8, 3, 2, -1,
    8, 2, 0, 1,
    8, 1, 0, 1,
    8, 0, 0, 1,
    9, 3, 3, -1,
    9, 2, 0, 1,
    9, 1, 0, 1,
    9, 0, 0, 1
), "N")

#############################
##
"""
Mori's asymmetric step-constrained pattern. Normalized in the
reference length.

Mori, A. Uchida, S. Kurazume, R. Taniguchi, R. Hasegawa, T. &
Sakoe, H. Early Recognition and Prediction of Gestures Proc. 18th
International Conference on Pattern Recognition ICPR 2006, 2006, 3,
560-563
"""
const mori2006 = StepPattern(_c(
    1, 2, 1, -1,
    1, 1, 0, 2,
    1, 0, 0, 1,
    2, 1, 1, -1,
    2, 0, 0, 3,
    3, 1, 2, -1,
    3, 0, 1, 3,
    3, 0, 0, 3
), "M")

"""
Completely unflexible: fixed slope 1. Only makes sense with
open.begin and open.end

N.B. open begin and open end are not implemented yet in Julia
"""
const rigid = StepPattern(_c(1, 1, 1, -1,
                       1, 0, 0, 1), "N")
##########################################################################################
##########################################################################################

function get_default_step_patterns()
    return (
        # fede symmetric pattern
        fede_symmetric,

        # fede asymmetric pattern
        fede_asymmetric,
        tr(fede_asymmetric),

        # dtw-python symmetric pattern
        symmetric1,
        symmetric2,
        symmetricP05,
        symmetricP1,
        symmetricP2,
        _symmetricVelichkoZagoruyko,

        # dtw-python asymmetric pattern
        asymmetric,
        tr(asymmetric),
        asymmetricP0,
        tr(asymmetricP0),
        asymmetricP05,
        tr(asymmetricP05),
        asymmetricP1,
        tr(asymmetricP1),
        asymmetricP2,
        tr(asymmetricP2),
        mori2006,
        tr(mori2006),
        _asymmetricItakura,
        tr(_asymmetricItakura),
    )
end
#! format: on

"""
Just prints the recursion rule from a pattern
"""
function recursion_rule(step_pattern::StepPattern)
    println("d[i, j] = min(")
    for rule in step_pattern.rules
        print("\t$(rule.weight) * [")
        for (i, comp) in enumerate(rule.components)
            print("$(comp.weight) * d[i-$(comp.row), j-$(comp.col)]")
            if i < length(rule.components)
                print(" + ")
            else
                print("] + $(rule.bias),\n")
            end
        end
    end
    return println(")")
end


text:
using MCMCChains, Test

ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "EightSchools")
cd(ProjDir) do

  isdir("tmp") &&
    rm("tmp", recursive=true);

  include(joinpath(ProjDir, "schools8.jl"))

  if rc == 0
    s = summarize(chns)
    @test s[:mu, :mean][1] ≈ 7.7 atol=1.0
  end

  isdir("tmp") &&
    rm("tmp", recursive=true);

end # cd


text:
function tracker(filename::String, lnn::LineNumberNode)
    return lock(Main.MicroCoverage_lock) do
        new_tracker_length = length(Main.MicroCoverage_tracker) + 1
        @debug("Tracker length: $(new_tracker_length)")
        push!(Main.MicroCoverage_tracker,
              false)
        push!(Main.MicroCoverage_tracker_linenumbernodes,
              LineNumberNode(lnn.line, filename))
        return :(Main.MicroCoverage_tracker[$(new_tracker_length)] = true)
    end
end

function insert_prepended_tracker(filename::String,
                                  expr::Expr,
                                  lnn::LineNumberNode;
                                  prepend_tracker::Bool = true)
    if prepend_tracker
        return Expr(:block, tracker(filename, lnn), expr)
    else
        return expr
    end
end


text:
function get_b1_bot(obj::FieldNeQuickParameters)
    return jcall(obj, "getB1Bot", RealFieldElement, ())
end

function get_b1_top(obj::FieldNeQuickParameters)
    return jcall(obj, "getB1Top", RealFieldElement, ())
end

function get_b2_bot(obj::FieldNeQuickParameters)
    return jcall(obj, "getB2Bot", RealFieldElement, ())
end

function get_be_bot(obj::FieldNeQuickParameters)
    return jcall(obj, "getBEBot", RealFieldElement, ())
end

function get_be_top(obj::FieldNeQuickParameters)
    return jcall(obj, "getBETop", RealFieldElement, ())
end

function get_h0(obj::FieldNeQuickParameters)
    return jcall(obj, "getH0", RealFieldElement, ())
end

function get_hm_e(obj::FieldNeQuickParameters)
    return jcall(obj, "getHmE", RealFieldElement, ())
end

function get_hm_f1(obj::FieldNeQuickParameters)
    return jcall(obj, "getHmF1", RealFieldElement, ())
end

function get_hm_f2(obj::FieldNeQuickParameters)
    return jcall(obj, "getHmF2", RealFieldElement, ())
end

function get_layer_amplitudes(obj::FieldNeQuickParameters)
    return jcall(obj, "getLayerAmplitudes", Vector{RealFieldElement}, ())
end

function get_nm_f2(obj::FieldNeQuickParameters)
    return jcall(obj, "getNmF2", RealFieldElement, ())
end


text:
include("iron_line_profile.jl")
include("arbitrary_f/iron_line_profile_approx.jl")

# x_vals_exact, bins_exact = iron_line_profile(
#                                             tolerance=1e-12,
#                                             # size_multiplier=3,
#                                             # fov=12,
#                                             # dtmax=5,
#                                             obs_angle=30,
#                                             output = "data"
#                                             )

vars = [-2,0,2]

labels = []
x_vals_vals = []
bins_vals = []

for var in vars
    x_vals_approx, bins_approx = iron_line_profile_approx(
                                                tolerance=1e-12,
                                                size_multiplier=3,
                                                fov=12,
                                                dtmax=2,
                                                obs_angle=30,
                                                output = "data",
                                                ϵ_3 = var ,
                                                spin = 0.8,
                                                normalised = false
                                                )
    push!(bins_vals, bins_approx)
    push!(x_vals_vals, x_vals_approx)
    push!(labels, "(ϵ_3 = $var)")
end
fixed_labels=permutedims(labels)
plt1 = plot(x_vals_vals, bins_vals, label=fixed_labels)
display(plt1)
# plot!(x_vals_exact, bins_exact, label="Page Thorne f",legend=:topleft)

labels = []
x_vals_vals = []
bins_vals = []

for var in vars
    x_vals_approx, bins_approx = iron_line_profile_approx(
                                                tolerance=1e-12,
                                                size_multiplier=3,
                                                fov=12,
                                                dtmax=2,
                                                obs_angle=30,
                                                output = "data",
                                                α_22 = var ,
                                                spin = 0.8 ,
                                                normalised = false
                                                )
    push!(bins_vals, bins_approx)
    push!(x_vals_vals, x_vals_approx)
    push!(labels, "(α_22 = $var)")
end
fixed_labels=permutedims(labels)
plt2 = plot(x_vals_vals, bins_vals, label=fixed_labels)
display(plt2)


text:
using FLoopsBase
using Test

@testset "FLoopsBase.jl" begin
    # Write your tests here.
end


text:
using Knet
using ArgParse
using JLD,MAT,JSON,HDF5

include("util.jl")
include("model.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "CNN-HMM hybrid ASR system in Knet by Ilker Kesen."

    @add_arg_table s begin
        ("--features"; required=true; help="MFSC features file")
        ("--jsonfile"; required=true)
        ("--model"; default=nothing)
        ("--probs"; default=nothing)
        ("--savefile"; default=nothing)
        ("--batchsize"; default=400; arg_type=Int64)
        ("--atype"; default=(gpu()>=0?"KnetArray{Float32}":"Array{Float32}"))
        ("--period"; default=100; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:atype] = eval(parse(o[:atype]))
    atype = o[:atype]

    # load data
    jsondata = JSON.parsefile(o[:jsonfile])
    p2i, i2p = build_vocab(jsondata)

    # make data
    trnids, tstids = [], []
    # trn = make_data(o[:features], jsondata; extra=trnids)
    tst = make_data(o[:features], jsondata, "test"; extra=tstids)

    # load model
    posteriors = load(o[:probs], "posteriors")
    priors = load(o[:probs], "priors")
    w = load(o[:model], "w")

    # convert atype
    # posteriors = convert(atype, posteriors)
    priors = convert(atype, priors)
    w = map(wi->convert(atype,wi), w)

    solutions = Dict()
    # @time for (split,data) in zip(("train","test"),(trn,tst))
    #  getid(i) = split=="train"?trnids[i]:tstids[i]
    data = tst
    getid(i) = tstids[i]
    for (i,sample) in enumerate(data)
        x = make_input(sample)
        ygold = make_output(sample, p2i)
        ypred = predict(w,convert(atype,x))
        ysoft = softmax(ypred)
        ynorm = ysoft ./ priors # acoustic model
        ynorm = Array(ynorm)

        # viterbi
        prob, track = viterbi(ynorm,posteriors)
        solutions[getid(i)] = Dict("logprob"=>prob,"states"=>map(t->i2p[t], track))

        println("$i generation so far.")
        if i % o[:period] == 0
            println("$i generation so far ", now())
        end
    end
    # end

    ff = open(abspath(o[:savefile]), "w")
    write(ff, JSON.json(solutions))
    close(ff)
    println("done")
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
module FizzBuzzParametricDispatch

export fizzbuzz

struct T{M,N} end

fizzbuzz(::T{0,N}) where N = "Fizz"
fizzbuzz(::T{M,0}) where M = "Buzz"
fizzbuzz(::T{0,0})         = "FizzBuzz"
fizzbuzz(::T)              = nothing

fizzbuzz(x) = something(fizzbuzz(T{x % 3, x % 5}()), string(x))

end


text:
function parse_ranking(s::String)
    d = Dict{String,Int}()
    if match(r"[\s,\(\)~∼≻]",s) === nothing
        Arr = split(s,"")
        for i = 1:length(Arr)
            haskey(d,Arr[i]) && @error("Duplicated identifier!")
            d[Arr[i]] = 1 - i
        end
    else
        parse_tokens!(s, d)
    end
    return d
end

function parse_tokens!(s::String, d::Dict{String,Int})
    in_equi_p = false
    weight_step = 1
    rel_weight = 0
    rk_expr = collect(tokenize(s))
    for tok in rk_expr
        if Tokens.kind(tok) == Tokens.IDENTIFIER
            token_string = Tokens.untokenize(tok)
            haskey(d,token_string) && @error("Duplicated identifier!")
            d[token_string] = rel_weight
            rel_weight -= weight_step
        elseif Tokens.kind(tok) == Tokens.LPAREN
            if in_equi_p
                @error("Unexpected opening parenthesis")
            else
                in_equi_p = true
                weight_step = 0
            end
        elseif Tokens.kind(tok) == Tokens.RPAREN
            if in_equi_p
                in_equi_p = false
                weight_step = 1
                rel_weight -= weight_step
            else
                @error("Unexpected closing parenthesis")
            end
        elseif Tokens.kind(tok) == Tokens.OP && Tokens.untokenize(tok) == "~"  ||
               Tokens.kind(tok) == Tokens.ERROR && Tokens.untokenize(tok) == "∼" # \sim
            rel_weight += weight_step
        else
            continue
        end
    end
    if in_equi_p
        @warn("Opened parenthesis is not closed")
    end
    if length(d) == 0
        @error("No identifier found!")
    end
end


text:
module PortAudioBindings
using Unitful
using Unitful: s

import Libdl

# Load in `deps.jl`, complaining if it does not exist

#const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
#if !isfile(depsjl_path)
#    error("PortAudioBindings was not build properly. Please run Pkg.build(\"PortAudioBindings\").")
#end
#include(depsjl_path)
## Module initialization function
#function __init__()
#    check_deps()
#end

using CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t
const libportaudio = "libportaudio"
include("libportaudio_manual.jl")
include("libportaudio_common.jl")
include("libportaudio_api.jl")

# export everything
foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "Pa") || startswith(string(s), "pa")
        @eval export $s
    end
end

using SampleArrays

export paerr, PortAudio, listdevices, device_id_by_name, device_id_default_in, device_id_default_out
export play, record, recordresponse

function paerr(err)
    if err != 0
        txt = unsafe_string(Pa_GetErrorText(err))
        @error "PortAudio Error ($err): $txt"
    end
    err
end

mutable struct PortAudio
    # TODO add some state variables, e.g., default devices, channels, sample rate?
    function PortAudio()
        err = Pa_Initialize() |> paerr
        if err != 0
            error("failed to initialize PortAudio!")
        end
        pa = new()

        function f(pa_)
            @async println("finalizing PortAudio")
            err = Pa_Terminate() |> paerr
            if err != 0
                error("failed to initialize PortAudio!")
            end
        end

        finalizer(f, pa)
    end
end

function listdevices(::PortAudio)
    ndevices = Pa_GetDeviceCount()
    defin = Pa_GetDefaultInputDevice()
    defout = Pa_GetDefaultOutputDevice()
    devs = NamedTuple[]
    for i in 0:ndevices-1
        devinfo = unsafe_load(Pa_GetDeviceInfo(i))
        push!(devs, (
            name = unsafe_string(devinfo.name),
            nins = devinfo.maxInputChannels,
            nouts = devinfo.maxOutputChannels,

            rate = devinfo.defaultSampleRate * u"Hz",
            inlowlatency = devinfo.defaultLowInputLatency * u"s",
            outlowlatency = devinfo.defaultLowOutputLatency * u"s",
            inhighlatency = devinfo.defaultHighInputLatency * u"s",
            outhighlatency = devinfo.defaultHighOutputLatency * u"s",
            defaultin = i == defin,
            defaultout = i == defout,
        ))
        deviceName = unsafe_string(devinfo.name)
    end
    devs
end

function device_id_by_name(pa::PortAudio, name::AbstractString)
    devs = listdevices(pa)
    id = findfirst(e -> e.name == name, devs)
    if isnothing(id)
        error("PortAudio: unknown device: $name")
    end
    id-1
end

function device_id_default_in(::PortAudio)
    Pa_GetDefaultInputDevice()
end

function device_id_default_out(::PortAudio)
    Pa_GetDefaultOutputDevice()
end

function play(::PortAudio, x::SampleArray; outdev=nothing, bsize=512)
    nchannels_, nframes_, rate_ = nchannels(x), nframes(x), rate(x)

    outdev = isnothing(outdev) ? Pa_GetDefaultOutputDevice() : outdev
    outinfo = unsafe_load(Pa_GetDeviceInfo(outdev));

    @debug "output: $(unsafe_string(outinfo.name)), #channels: $(nchannels_)"

    outParams = PaStreamParameters(outdev, nchannels_, paFloat32,
        outinfo.defaultHighOutputLatency, C_NULL)

    stream = Ref(Ptr{PaStream}(C_NULL))

    Pa_OpenStream(
        stream,
        C_NULL,
        Ref(outParams),
        rate_,
        bsize,
        paClipOff,
        C_NULL,
        C_NULL) |> paerr

    samplebuf = zeros(Float32, nchannels_, bsize)

    Pa_StartStream(stream[]) |> paerr

    samples = data(x)'

    off = 1
    while off <= nframes_
        high = min(off + bsize -1, nframes_)
        highbuf = high-off+1
        samplebuf[:, 1:highbuf] .= samples[:, off:high]
        samplebuf[:, highbuf+1:end] .= 0
        Pa_WriteStream(stream[], pointer(samplebuf), bsize) |> paerr
        off += bsize
    end

    Pa_StopStream(stream[]) |> paerr
    Pa_CloseStream(stream[]) |> paerr
end

function record(::PortAudio; s=1.0, rate=44100.0, indev=nothing, nins=nothing, bsize=512)
    nblocks = ceil(Int, s*rate/bsize)

    indev = isnothing(indev) ? Pa_GetDefaultInputDevice() : indev
    ininfo = unsafe_load(Pa_GetDeviceInfo(indev))
    nins = isnothing(nins) ? ininfo.maxInputChannels : nins

    inparams = PaStreamParameters(indev, nins, paFloat32,
        ininfo.defaultHighInputLatency, C_NULL)

    stream = Ref(Ptr{PaStream}(C_NULL))

    err = Pa_OpenStream(
        stream,
        Ref(inparams),
        C_NULL,
        rate,
        bsize,
        paClipOff,
        C_NULL,
        C_NULL) |> paerr

    if err != 0
        return
    end

    samplebuf = zeros(Float32, nins, bsize)
    samples = zeros(Float32, nins, nblocks * bsize)

    err = Pa_StartStream(stream[]) |> paerr

    off = 1
    for i in 1:nblocks
        Pa_ReadStream(stream[], pointer(samplebuf), bsize) |> paerr
        samples[:, off:off+bsize-1] .= samplebuf
        off += bsize
    end

    Pa_StopStream(stream[]) |> paerr
    Pa_CloseStream(stream[]) |> paerr

    return SampleArray(samples', rate)
end

function recordresponse(::PortAudio, x::SampleArray; indev=nothing, outdev=nothing, nins=nothing, post::Time=1.0s, bsize=512)
    post = tos(post)

    nouts, nframes_, rate_ = nchannels(x), nframes(x), rate(x)

    indev = isnothing(indev) ? Pa_GetDefaultInputDevice() : indev
    ininfo = unsafe_load(Pa_GetDeviceInfo(indev));
    outdev = isnothing(outdev) ? Pa_GetDefaultOutputDevice() : outdev
    outinfo = unsafe_load(Pa_GetDeviceInfo(outdev));
    nins = isnothing(nins) ? ininfo.maxInputChannels : nins

    inparams = PaStreamParameters(indev, nins, paFloat32, ininfo.defaultHighInputLatency, C_NULL)
    @debug "input: $(unsafe_string(ininfo.name)), #channels: $(nins)"

    outparams = PaStreamParameters(outdev, nouts, paFloat32, ininfo.defaultHighOutputLatency, C_NULL)
    @debug "output: $(unsafe_string(outinfo.name)), #channels: $(nouts)"

    stream = Ref(Ptr{PaStream}(C_NULL))

    err = Pa_OpenStream(
        stream,
        Ref(inparams),
        Ref(outparams),
        rate_,
        bsize,
        paClipOff,
        C_NULL,
        C_NULL) |> paerr

    if err != 0
        return
    end

    nframesrec = nframes_ + ceil(Int, post * rate_)
    @debug "input signal frames: $nframes_, frames to record: $nframesrec"
    sampleinbuf = zeros(Float32, nins, bsize)
    sampleoutbuf = zeros(Float32, nouts, bsize)
    y = similar(x, nframesrec, nins)
    samples = data(y)'
    samples .= 0

    Pa_StartStream(stream[]) |> paerr

    insignal = data(x)'

    off = 1
    werr = 0
    woffs = []
    rerr = 0
    while off <= nframes_
        high = min(off + bsize -1, nframes_)
        highbuf = high-off+1
        sampleoutbuf[:, 1:highbuf] .= insignal[:, off:high]
        sampleoutbuf[:, highbuf+1:end] .= 0
        # @show high
        # @show highbuf
        # @show off
        # @show 1:highbuf, off:high

        err = Pa_WriteStream(stream[], pointer(sampleoutbuf), bsize)
        if err != 0
            werr = err
            push!(woffs, off)
        end

        err = Pa_ReadStream(stream[], pointer(sampleinbuf), bsize)
        if err != 0
            rerr = err
        end

        samples[:, off:off+bsize-1] .= sampleinbuf
        # if off > 1024
            # break
        # end
        off += bsize
    end

    while off <= nframesrec
        err = Pa_ReadStream(stream[], pointer(sampleinbuf), bsize)
        if err != 0
            rerr = err
        end
        high = min(off + bsize -1, nframesrec)
        highbuf = high-off+1
        samples[:, off:high] .= sampleinbuf[:, 1:highbuf]

        off += bsize
    end

    paerr(werr)
    @error "offsets: $woffs"
    paerr(rerr)

    Pa_StopStream(stream[]) |> paerr
    i = 1
    while true
        err = Pa_IsStreamActive(stream[])
        if err != 1
            break
        end
        @show i, err
        Pa_Sleep(1000);
        i += 1
    end

    # Pa_AbortStream(stream[]) |> paerr
    Pa_CloseStream(stream[]) |> paerr

    y
end

end # module


text:
module BatchIterators

using Statistics

export BatchIterator
export choose_batchsize
export centered_batch_iterator

"""
	BatchIterator(X; batchsize = nothing, limit=size(X,2))

Wrapper allowing to iterate over batches of `batchsize` columns of `X`. `X` can be of any type supporting `size` and 2d indexing. When `limit` is provided, iteration is restricted to the columns of `X[:, 1:limit]`.
"""
struct BatchIterator{T}
	X::T
	length::Int		# Number of batches
	bsz::Int		# Batch size
	limit::Int
	function BatchIterator(X; batchsize=nothing, limit=size(X,2))
		@assert limit > 0 && limit ≤ size(X,2)
		bsz = (batchsize == nothing) ? choose_batchsize(size(X,1), limit) : batchsize
		nb = ceil(Int, limit/bsz)
		new{typeof(X)}(X, nb, bsz, limit)
	end
end

view_compatible(::Any) = false
view_compatible(::Array) = true
view_compatible(bi::BatchIterator) = view_compatible(bi.X)

#######################################################################
#                             Iteration                               #
#######################################################################

function Base.getindex(it::BatchIterator, i)
	d = i - it.length		# > 0 means overflow, == 0 means last batch
	cbsz = (d == 0) ? mod(it.limit - 1, it.bsz) + 1 : it.bsz		# Size of current batch
	if (i<1 || d > 0)
		@error "Out of bounds."
	else
		# TODO using views here might impact type stability.
		view_compatible(it) ? (@view it.X[:, (i-1)*it.bsz+1:(i-1)*it.bsz+cbsz]) : it.X[:, (i-1)*it.bsz+1:(i-1)*it.bsz+cbsz]
	end
end
Base.length(it::BatchIterator)  = it.length
function Base.iterate(it::BatchIterator{T}, st = 0) where T
	st = st + 1				# new state
	d = st - it.length		# > 0 means overflow, == 0 means last batch
	(d > 0) ?  nothing : (it[st], st)
end

"""
	centered_batch_iterator(X; kwargs...)

Similar to BatchIterator, but performs first one pass over the data to compute the mean, and centers the batches.
"""
function centered_batch_iterator(X; kwargs...)
	bi = BatchIterator(X; kwargs...)
	μ = vec(mean(mean(b, dims=2) for b in BatchIterator(X)))
	(b .- μ for b in bi)
end

#######################################################################
#                              Utilities                              #
#######################################################################

"""
	choose_batchsize(d, n; maxmemGB = 1.0, maxbatchsize = 2^14, sizeoneB = d*sizeof(Float64))

Computes the size (nb. of columns) of a batch, so that each column of the batch can be converted to a vector of size `sizeoneB` (in bytes) with a total memory constrained by `maxmemGB` (gigabytes).
"""
function choose_batchsize(d, n;
				   maxmemGB = 1.0,
				   maxbatchsize = 2^14,
				   sizeoneB = d*sizeof(Float64),
				   forcepow2 = true)

	fullsizeGB = n * sizeoneB/1024^3		# Size of the sketches of all samples
	batchsize = (fullsizeGB > maxmemGB) ? ceil(Int, n/ceil(Int, fullsizeGB/maxmemGB)) : n
	batchsize = min(batchsize, maxbatchsize)
	(forcepow2 && batchsize != n) ? prevpow(2, batchsize) : batchsize
end

end # module


text:
export MixedSystem

"""
    MixedSystem

A data structure which uses for `evaluate!` and `evaluate_and_jacobian!` a [`CompiledSystem`](@ref)
and for `taylor!` an [`InterpretedSystem`](@ref).
"""
struct MixedSystem{T,ID} <: AbstractSystem
    compiled::CompiledSystem{ID}
    interpreted::InterpretedSystem{T}
end
MixedSystem(F::System; kwargs...) =
    MixedSystem(CompiledSystem(F; kwargs...), InterpretedSystem(F; kwargs...))

Base.size(F::MixedSystem) = size(F.compiled)
ModelKit.variables(F::MixedSystem) = variables(F.interpreted)
ModelKit.parameters(F::MixedSystem) = parameters(F.interpreted)
ModelKit.variable_groups(F::MixedSystem) = variable_groups(F.interpreted)
ModelKit.System(F::MixedSystem) = System(F.interpreted)
Base.:(==)(F::MixedSystem, G::MixedSystem) = F.interpreted == G.interpreted

function Base.show(io::IO, F::MixedSystem)
    print(io, "Mixed: ")
    show(io, System(F))
end

(F::MixedSystem)(x, p = nothing) = F.interpreted(x, p)
ModelKit.evaluate!(u, F::MixedSystem, x::AbstractVector, p = nothing) =
    evaluate!(u, F.compiled, x, p)
ModelKit.evaluate_and_jacobian!(u, U, F::MixedSystem, x, p = nothing) =
    evaluate_and_jacobian!(u, U, F.compiled, x, p)
ModelKit.taylor!(u, v::Val, F::MixedSystem, tx, p = nothing) =
    taylor!(u, v, F.interpreted, tx, p)


text:
"""
Solving the travelling salesman problem.

if you have a costs square symetric matrix with diagonal equal to zero use the function solve_TSP()
if you have the set of (x,y) coordinates then use the function hacematrixdedistancias() first to generate the costs matrix

The function solve_TSP give the minimun costs posible and the elements of the costs matrix that minimize the costs the travel
"""
function solve_TSP(costs_matrix)
    N=size(costs_matrix)[1]
    if N==size(costs_matrix)[2];
        m = JuMP.Model(JuMP.with_optimizer(GLPK.Optimizer))
        @variable(m, x[1:N,1:N], Bin)

        @objective(m, Min, sum(costs_matrix.*x))
        @constraint(m, notself[i=1:N], x[i,i] == 0)
        @constraint(m, goto[i=1:N], sum(x[i,1:N]) == 1)
        @constraint(m, cometo[j=1:N], sum(x[1:N,j]) == 1)
        @variable(m, t[1:N], Int)
        for i=1:N,j=2:N
            if i!=j
                @constraint(m,t[j]>=t[i]+1-N*(1-x[i,j]))
            end
        end

        optimize!(m)
        println("El coste optimo es: ")
        println(objective_value(m))
        for i=1:N
            for j=1:N
                if value(x[i,j])!=0
                    println("d[$i,$j]=",costs_matrix[i,j])
                end
            end
        end
    else println("La matriz debe ser cuadrada, con diagonal principal =0 y simetrica")
    end

end

function hacematrixdedistancias(listadeptos)
    distmatrix=zeros(Float64,(length(listadeptos),length(listadeptos)));
    for i=1:length(listadeptos)
        for j=1:length(listadeptos)
            distmatrix[i,j]=distanciacartesiana(listadeptos[i],listadeptos[j])
        end
    end
    return distmatrix;
end

function distanciacartesiana(pto1,pto2)
    dis=sqrt((pto1[1]-pto2[1])^2+(pto1[2]-pto2[2])^2);
    return dis;
end


text:
using PackageCompiler

create_sysimage(:Poisson,
  sysimage_path=joinpath(@__DIR__,"..","Poisson.so"),
  precompile_execution_file=joinpath(@__DIR__,"..","test","runtests.jl"))


text:
import Tokenize.Lexers: peekchar, prevchar, readchar, iswhitespace, emit, emit_error, backup!, accept_batch, eof

const EmptyWS = Tokens.begin_delimiters
const SemiColonWS = Tokens.end_delimiters
const NewLineWS = Tokens.begin_literal
const WS = Tokens.end_literal
const InvisibleBrackets = Tokens.begin_invisible_keywords
const EmptyWSToken = Token(EmptyWS, (0, 0), (0, 0), -1, -1, "")

"""
    Closer
Struct holding information on the tokens that will close the expression
currently being parsed.
"""
mutable struct Closer
    toplevel::Bool
    newline::Bool
    semicolon::Bool
    tuple::Bool
    comma::Bool
    dot::Bool
    paren::Bool
    quotemode::Bool
    brace::Bool
    inmacro::Bool
    insquare::Bool
    inwhere::Bool
    square::Bool
    block::Bool
    ifelse::Bool
    ifop::Bool
    range::Bool
    trycatch::Bool
    ws::Bool
    wsop::Bool
    precedence::Int
    stop::Int
end
Closer() = Closer(true, true, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, -1, typemax(Int))

"""
    ParseState

The parser's interface with `Tokenize.Lexers.Lexer`. This alters the output of `tokenize` in three ways:
+ Merging of whitespace, comments and semicolons;
+ Skips DOT tokens where they are followed by another operator and marks the trailing operator as dotted.
+ Keeps track of the previous, current and next tokens.

In addition a list of formatting hints is updated as the parsing progresses
and copy of the current `Scope`.
"""
mutable struct ParseState
    l::Lexer
    done::Bool
    lt::Token
    t::Token
    nt::Token
    nnt::Token
    lws::Token
    ws::Token
    nws::Token
    nnws::Token
    dot::Bool
    ndot::Bool
    trackscope::Bool
    formatcheck::Bool
    diagnostics::Vector{Diagnostics.Diagnostic}
    closer::Closer
    errored::Bool
    current_scope
end
function ParseState(str::String)
    ps = ParseState(tokenize(str), false, Token(), Token(), Token(), Token(), Token(), Token(), Token(), Token(), true, true, true, true, [], Closer(), false, Scope{Tokens.TOPLEVEL})
    return next(next(ps))
end

function Base.show(io::IO, ps::ParseState)
    println(io, "ParseState $(ps.done ? "finished " : "")at $(position(ps.l.io))")
    println(io, "last    : ", ps.lt.kind, " ($(ps.lt))", "    ($(wstype(ps.lws)))")
    println(io, "current : ", ps.t.kind, " ($(ps.t))", "    ($(wstype(ps.ws)))")
    println(io, "next    : ", ps.nt.kind, " ($(ps.nt))", "    ($(wstype(ps.nws)))")
end
peekchar(ps::ParseState) = peekchar(ps.l)
wstype(t::Token) = t.kind == EmptyWS ? "empty" :
                   t.kind == NewLineWS ? "ws w/ newline" :
                   t.kind == SemiColonWS ? "ws w/ semicolon" : "ws"

function next(ps::ParseState)
    #  shift old tokens
    ps.lt = ps.t
    ps.t = ps.nt
    ps.nt = ps.nnt
    ps.lws = ps.ws
    ps.ws = ps.nws
    ps.nws = ps.nnws
    ps.dot = ps.ndot

    ps.nnt, ps.done  = next(ps.l, ps.done)
    # Reject new kws for now
    # if ps.nnt.kind == Tokens.STRUCT || ps.nnt.kind == Tokens.STRUCT || ps.nnt.kind == Tokens.MUTABLE
    #     ps.nnt = Token(Tokens.IDENTIFIER, ps.nnt.startpos, ps.nnt.endpos, ps.nnt.startbyte, ps.nnt.endbyte, ps.nnt.val)
    # end

    # Handle dotted operators
    if ps.nt.kind == Tokens.DOT && ps.nws.kind == EmptyWS && isoperator(ps.nnt) && !non_dotted_op(ps.nnt)
        # ps.nt = ps.nnt
        ps.nt = Token(ps.nnt.kind, (ps.nnt.startpos[1], ps.nnt.startpos[2] - 1), ps.nnt.endpos, ps.nnt.startbyte - 1, ps.nnt.endbyte, ps.nnt.val)
        ps.ndot = true
        # combines whitespace, comments and semicolons
        if iswhitespace(peekchar(ps.l)) || peekchar(ps.l) == '#' || peekchar(ps.l) == ';'
            ps.nws = lex_ws_comment(ps.l, readchar(ps.l))
        else
            ps.nws = Token(EmptyWS, (0, 0), (0, 0), ps.nnt.endbyte, ps.nnt.endbyte, "")
        end
        ps.nnt, _ = next(ps.l, ps.done)
    else
        ps.ndot = false
    end
    # combines whitespace, comments and semicolons
    if iswhitespace(peekchar(ps.l)) || peekchar(ps.l) == '#' || peekchar(ps.l) == ';'
        ps.nnws = lex_ws_comment(ps.l, readchar(ps.l))
    else
        ps.nnws = EmptyWSToken
    end
    ps.done = ps.nt.kind == Tokens.ENDMARKER
    return ps
end

"""
    lex_ws_comment(l::Lexer, c)

Having hit an initial whitespace/comment/semicolon continues collecting similar
`Chars` until they end. Returns a WS token with an indication of newlines/ semicolons. Indicating a semicolons takes precedence over line breaks as the former is equivalent to the former in most cases.
"""
function lex_ws_comment(l::Lexer, c)
    newline = c == '\n'
    semicolon = c == ';'
    if c == '#'
        newline = read_comment(l)
    else
        newline, semicolon = read_ws(l, newline, semicolon)
    end
    while iswhitespace(peekchar(l)) || peekchar(l) == '#' || peekchar(l) == ';'
        c = readchar(l)
        if c == '#'
            read_comment(l)
            newline = newline || peekchar(l) == '\n'
            semicolon = semicolon || peekchar(l) == ';'
        elseif c == ';'
            semicolon = true
        else
            newline, semicolon = read_ws(l, newline, semicolon)
        end
    end

    return emit(l, semicolon ? SemiColonWS :
                   newline ? NewLineWS : WS)
end

function read_ws(l::Lexer, newline, semicolon)
    while iswhitespace(peekchar(l))
        c = readchar(l)
        c == '\n' && (newline = true)
        c == ';' && (semicolon = true)
    end
    return newline, semicolon
end

function read_comment(l::Lexer)
    if peekchar(l) != '='
        c = readchar(l)
        if c == '\n' || eof(c)
            backup!(l)
            return true
        end
        while true
            c = readchar(l)
            if c == '\n' || eof(c)
                backup!(l)
                return true
            end
        end
    else
        c = readchar(l) # consume the '='
        n_start, n_end = 1, 0
        while true
            if eof(c)
                emit_error(l, Tokens.EOF_MULTICOMMENT)
                return false
            end
            nc = readchar(l)
            if c == '#' && nc == '='
                n_start += 1
            elseif c == '=' && nc == '#'
                n_end += 1
            end
            if n_start == n_end
                return false
            end
            c = nc
        end
    end
end

isemptyws(t::Token) = t.kind == EmptyWS


text:
using Documenter, BandedMatrices

makedocs(;
    modules = [BandedMatrices],
    format = Documenter.HTML(
        canonical = "https://JuliaMatrices.github.io/BandedMatrices.jl/stable/",
    ),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.com/JuliaMatrices/BandedMatrices.jl/blob/{commit}{path}#L{line}",
    sitename = "BandedMatrices.jl",
    authors = "Sheehan Olver, Mikael Slevinsky, and contributors.",
    assets = String[],
)

deploydocs(;
    repo   = "github.com/JuliaMatrices/BandedMatrices.jl.git"
    )


text:
# # Searching a mesh

# In this tutorial, we will learn
#
#    -  How to import an incidence relation from a Nastran file.
#    -  How to create a mesh from this incidence relation.
#    -  How to print the summary of the mesh.
#    -  How to search for some vertices.
#    -  How to export the selected vertices for visualization.

# The tutorial will produce a visualization in the
# [Paraview](https://www.paraview.org/) format (VTK). One can display this mesh
# by loading the file with `paraview.exe`. When the tutorial is executed in
# `mybinder.org`, the graphics file needs to be downloaded to your desktop, and
# then visualized locally.

# First  we will use a `MeshSteward` function to import a NASTRAN mesh file.
# This file stores a tetrahedral mesh of a hollow cylinder.
using MeshSteward: import_NASTRAN
connectivities = import_NASTRAN("trunc_cyl_shell_0.nas");

# The result, `connectivities`, is an an array of incidence relations (in this
# case there is just one entry in this array). We can create a mesh
# from this incidence relation as
using MeshSteward: Mesh, attach!
mesh = Mesh()
attach!(mesh, connectivities[1]);

# Let us print the summary information for the mesh:
using MeshSteward: summary
println(summary(mesh))

# Here `(3, 0)` means relationship between three-dimensional finite elements and
# the zero dimensional vertices (dimension in the sense of manifolds). It is
# referred to as the "code" of the incidence relation. The final piece of
# information shows the attributes of the vertex shape collection: `geom`
# provides access to the locations of the vertices.

# The imported mesh can be exported, this time in a different format,
# suitable for visualization. The tetrahedral elements are the base incidence
# relation of the mesh.
using MeshSteward: baseincrel
using MeshSteward: vtkwrite
vtkwrite("trunc_cyl_shell_0-elements", baseincrel(mesh))

# Start "Paraview", load the file `"trunc_cyl_shell_0-elements.vtu"` and
# select for instance view as "Surface with Edges". The result will be a view
# of the surface of the tetrahedral mesh.
@async run(`paraview trunc_cyl_shell_0-elements.vtu`)

# Now we will locate some vertices. In particular, let us say we wish to select
# the vertices  in one of the cross-sections of the hollow cylinder, one at $$Z
# = 0.5$$. We shall use the method of a bounding box:  all vertices inside the
# bounding box will be selected. The bounding box will initially start as
# totally flat (zero volume), and in order to capture the vertices which are not
# *precisely* at $$Z = 0.5$$ (given the vagaries of floating-point numbers), we
# shall inflate the box to all sides by a small amount.
using MeshSteward: boundingbox, inflatebox!
box = boundingbox([-Inf -Inf 0.5; Inf Inf 0.5])
box = inflatebox!(box, 0.001)

# The function `vselect` will create an incidence relation that consists of the
# selected vertices.
using MeshSteward: vselect
selectedv = vselect(mesh, box = box);

# This is how many vertices  have been found in that cross-section:
using MeshCore: nshapes
@show nshapes(selectedv.left)

# The incidence relation links some vertices (the shape collection on the left)
# to the shape collection of the original vertices (on the right). We can see it
# from the summary:
summary(selectedv)

# The incidence relation consisting of the selected vertices may be exported for
# viewing with the "Paraview" visualization program. Select to visualize this
# mesh as
# "Points", and select the size of the points as 8.
using MeshSteward: vtkwrite
vtkwrite("trunc_cyl_shell_0-selected-vertices", selectedv)

# Note that a file of this name, `"trunc_cyl_shell_0-selected-vertices.vtu"` is
# provided with the installation. When you run this tutorial  you should get a
# file with the same contents. The selected  vertices may be visualized with the
# point glyphs ("Point Gaussian") in Paraview. Provided the Paraview executable
# is in your path, the line below should start it.
@async run(`paraview trunc_cyl_shell_0-selected-vertices.vtu`)

# It is also possible to start `paraview` and then manually load both
# files, `"trunc_cyl_shell_0-elements.vtu"` and
# `"trunc_cyl_shell_0-selected-vertices.vtu"` into a single visualization.

# The visualization may look like this:
# ![Visualization of the elements and selected vertices](t3_i1.png)


text:
#==========================================================================
  LDA Solvers
==========================================================================#

function class_counts{T<:Integer}(y::Vector{T}, k::T)
    counts = zeros(Int64, k)
    for i = 1:length(y)
        y[i] <= k || error("Index $i out of range.")
        counts[y[i]] += 1
    end
    counts
end

function class_totals{T<:AbstractFloat,U<:Integer}(X::Matrix{T}, y::Vector{U}, k::U)
    n, p = size(X)
    length(y) == n || throw(DimensionMismatch("X and y must have the same number of rows."))
    M = zeros(T, k, p)
    for j = 1:p, i = 1:n
        M[y[i],j] += X[i,j]
    end
    M
end

function center_rows!{T<:AbstractFloat,U<:Integer}(X::Matrix{T}, M::Matrix{T}, y::Vector{U})
    n, p = size(X)
    for j = 1:p, i = 1:n
        X[i,j] -= M[y[i],j]
    end
    X
end

# Takes intermediate matrices such that H_b'H_b = Σ_b, H_w'H_w = Σ_w
# and returns (lda_components, eigenvalues)
function lda_components!{T<:AbstractFloat}(H_b::Matrix{T}, α_b::T, ϵ_b::T, H_w::Matrix{T}, α_w::T, ϵ_w::T)
    Σ_b = syml!(BLAS.syrk('U', 'T', one(T), H_b))  # Σ_b = H_b'*diag(freq)*H_b
    α_b == 0 || regularize!(Σ_b, α_b, trace(Σ_b)/size(Σ_b,1))
    ϵ_b == 0 || perturb!(Σ_b, ϵ_b)
    #tol = eps(T) * maximum(size(H_w)) * maximum(H_w)
    Σ_w = syml!(BLAS.syrk('U', 'T', one(T), H_w))  # Σ_w = H_w'H_w
    α_w == 0 || regularize!(Σ_w, α_w, trace(Σ_w)/size(Σ_w,1))
    ϵ_w == 0 || perturb!(Σ_w, ϵ_w)
    components_geig!(Σ_b, Σ_b)
end

function lda_matrices!{T<:AbstractFloat,U<:Integer}(X::Matrix{T}, y::Vector{U}, k::U, frequency::Vector{T})
    n, p = size(X)
    M = class_totals(X, y, k)
    scale!(one(T) ./ class_counts(y, k), M)
    H_w = center_rows!(X, M, y)  # M = [μ1; μ2; ...]
    scale!(H_w, one(T)/(n-k))  # sampling correction factor for Σ_w
    μ = vec(frequency'M)
    translate!(M,-μ)  # M:= M .- μ'
    H_b = scale!(sqrt(frequency), M)  # M := freq .* M (applies weights to class mean)
    (H_b, H_w, μ)
end

function lda!{T<:AbstractFloat,U<:Integer}(X::Matrix{T}, y::Vector{U}, k::U, frequency::Vector{T},
                                           α_b::T = zero(T), ϵ_b::T = zero(T), α_w::T = zero(T), ϵ_w::T = zero(T))
    H_b, H_w = lda_matrices!(X, y, k, frequency)
    lda_components!(H_b, α_b, ϵ_b, H_w, α_w, ϵ_w)
end

# W is components returned
# Z must be the new matrix to be transformed
transform_lda{T<:AbstractFloat}(W::Matrix{T}, Z::Matrix{T}) = Z * W

# X must be original data matrix
# κ must be original kernel
# W is components returned
# Z must be the new matrix to be transformed
transform_klda{T<:AbstractFloat}(X::Matrix{T}, κ::Kernel{T}, W::Matrix{T}, Z::Matrix{T}) = kernelmatrix(Z, X) * W


text:
@testset "model.jl" begin

end