text:
using OpticalPropagation
using Test
using Unitful
using RecipesBase

@testset "OpticalPropagation.jl" begin
    @testset "MonoLightField2D Constructors" begin
        o=MonoLightField2D([1 2; 3 4],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
        @test begin
            ro=MonoLightField2D(o,data=[3 1; 2 4])
            r=MonoLightField2D([3 1; 2 4],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            ro≈r
        end
        @test begin
            ro=MonoLightField2D(o,wavelength=1.15u"μm")
            r=MonoLightField2D([1 2; 3 4],wavelength=1.15u"μm",size=(1u"mm",1u"mm"))
            ro≈r
        end
        @test begin
            ro=MonoLightField2D(o,size=(1.5u"mm",2.0u"mm"))
            r=MonoLightField2D([1 2; 3 4],wavelength=632.8u"nm",size=(1.5u"mm",2.0u"mm"))
            ro≈r
        end
        @test begin
            ro=MonoLightField2D(o,data=[3 1; 2 4],wavelength=1.15u"μm")
            r=MonoLightField2D([3 1; 2 4],wavelength=1.15u"μm",size=(1u"mm",1u"mm"))
            ro≈r
        end
        @test begin
            ro=MonoLightField2D(o,wavelength=1.15u"μm",size=(1.5u"mm",2.0u"mm"))
            r=MonoLightField2D([1 2; 3 4],wavelength=1.15u"μm",size=(1.5u"mm",2.0u"mm"))
            ro≈r
        end
    end
    @testset "Additive Operator $op" for op in [(+), (-)]
        @test begin
            ad = [1 2; 3 4]
            bd = [1.0 -2; 3im -5]
            a=MonoLightField2D(ad,wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            b=MonoLightField2D(bd,wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            c=MonoLightField2D(op(ad, bd),wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            op(a,b)≈c
        end
        @test_throws DimensionMismatch begin
            a=MonoLightField2D([1 2; 3 4],wavelength=589.6u"nm",size=(1u"mm",1u"mm"))
            b=MonoLightField2D([1.0 -2; 3im -5],wavelength=589u"nm",size=(1u"mm",1u"mm"))
            op(a,b)
        end
        @test_throws DimensionMismatch begin
            a=MonoLightField2D([1 2; 3 4],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            b=MonoLightField2D([1.0 -2; 3im -5],wavelength=632.8u"nm",size=(1.1u"mm",1.1u"mm"))
            op(a,b)
        end
    end
    @testset "Multiplicative Operator" begin
        @test begin
            a=2+3im
            b=MonoLightField2D([1.0 -2; 3im -1+4im],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            c=MonoLightField2D([2.0+3.0im -4.0-6.0im; -9.0+6.0im -14.0+5.0im],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            a*b≈b*a≈c
        end
        @test begin
            a=MonoLightField2D([1.0 -2; 3im -1+4im],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            b=1.5+0.5im
            c=MonoLightField2D([0.6-0.2im -1.2+0.4im; 0.6+1.8im 0.2+2.6im],wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
            a/b≈b\a≈c
        end
    end
    @testset "Array-like" begin
        a=MonoLightField2D(rand(4,5),wavelength=632.8u"nm",size=(1u"mm",1u"mm"))
        @test size(a) == size(a.data)
        @test a[1,1] == a.data[1,1]
        @test begin
            a[2,3] = 10
            a[2,3] == 10
        end
        @test abs.(a) == MonoLightField2D(a, data=abs.(a.data))
    end
    @testset "Plot Recipes" begin
        import RecipesBase.is_key_supported # We hack the RecipesBase to avoid MethodError when apply_recipe
        RecipesBase.is_key_supported(k::Symbol) = false
        l=MonoLightField2D(ones(4,5),wavelength=632.8u"nm",size=(1u"mm",1.2u"mm"))
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), l)
        @test getfield(rec[1],1)[:aspect_ratio] == :equal
        @test getfield(rec[1],1)[:size] == (4, 5)
        @test getfield(rec[1],1)[:xguide] == getfield(rec[1],1)[:yguide] == u"mm"
        @test rec[1].args == (-0.5:0.25:0.25, -0.6:0.24:0.36, ones(4,5))
    end
    @testset "Scale Invariant Propagation: $p" for p in [angularspectrum, fresnel2]
        @test begin # plane wave
            λ=632.8e-9
            lx=1.2e-3
            ly=1e-3
            d=1e-2
            Uin = ones(ComplexF64, 1024, 768)
            Uout = Uin*exp(2im*pi*d/λ)
            p(Uin, d, λ, lx, ly)≈Uout
        end
        @test begin # plane wave
            λ=632.8u"nm"
            lx=1.2u"mm"
            ly=1u"mm"
            d=1u"cm"
            Uin = MonoLightField2D(ones(ComplexF64, 1024, 768),wavelength=λ,size=(lx,ly))
            Uout = Uin*exp(2im*pi*convert(Float64, d/λ))
            p(Uin, d)≈Uout
        end
        @test begin # distance superposition
            λ=632.8u"nm"
            lx=ly=1u"mm"
            d1=1u"cm"
            d2=25u"mm"
            Uin = MonoLightField2D(rand(ComplexF64, 512, 512),wavelength=λ,size=(lx,ly))
            (p(d2)∘p(d1))(Uin)≈(p(d1)∘p(d2))(Uin)≈p(Uin, d1+d2)
        end
    end
    @testset "Fresnel S-FFT" begin
        @test begin # sample rate
            λ=632.8e-9
            lx=1.2e-3
            ly=1e-3
            d=1e-2
            nx = 16
            ny = 9
            Uin = ones(ComplexF64, nx, ny)
            (_, (lxo, lyo)) = fresnel1(Uin, d, λ, lx, ly)
            lxo≈λ*d*nx/lx && lyo≈λ*d*ny/ly
        end
    end
end


text:
export Trend

"""
    Trend(t₀, slope)

Constant `slope` starting at `start`, i.e., `ϕ(t) = (t - t₀) * slope`.
"""
mutable struct Trend <: Sequence
  start
  slope
end

sample(seq::Trend, t::AbstractVector) = (t .- seq.start) * seq.slope

ext: jl


text:
ENV_name = "local"
ENV_root = "output"
println("Environment: Local")

"""Launches the script via qsub, in case dir doesnt exist, create it"""
function launch_job(dir::String, Ncores::Integer, gpu::Bool, julia_script::String, successive_jobs::Integer=1; architecture::String="", hours::Int64=48, fresh::Bool=false, memory::String="", gpu_mem::String="")
  full_dir = "$(ENV["DETERMINATION_DATA"])/$(ENV_root)/$dir"
  run(`mkdir -p $(full_dir)`)
  open("$(full_dir)/job.jl", "w+") do file
    write(file,"""
    $julia_script
    flush(STDOUT) #make sure output is flushed to terminal files
    """)
  end
  cd("$(full_dir)")
  include("$(pwd())/job.jl")
end


text:
###Hastings & Powell 1991
### Reproducing figures 2

u0 = [0.76, 0.16, 9.9]  #initial conditions for x,y,z
tspan = (0.0, 10000.0) #timespan
p = [5.0, 0.1, 3.0, 2.0, 0.4, 0.01] # values for a1, a2, b1, b2, d1 and d2

###Solving the differential equations
prob = ODEProblem(parameterized_model!, u0, tspan, p)
sol = solve(prob)

###Figure 2
fig2a = plot(sol, vars = (0, 1), xlim = (5000, 6500), ylim = (0, 1),
             xlabel = "time", ylabel = "x", leg = false,
             grid = :none, title = "a", titleloc = :left, c = :black
             ) #2a
fig2b = plot(sol, vars = (0, 2), xlim = (5000, 6500), ylim = (0, 0.5),
             xlabel = "time", ylabel = "y", leg = false,
             grid = :none, title="b", titleloc=:left, c=:black
             ) #2b
fig2c = plot(sol, vars = (0, 3), xlim = (5000, 6500), ylim = (7, 10.5),
             xlabel = "time", ylabel = "z", leg = false,
             grid = :none, title = "c", titleloc = :left, c = :black,
             ) #2c

fig2 = plot(fig2a, fig2b, fig2c, layout = (3,1),
            size = (500, 900), titlefontsize = 11,
            tickfontsize = 8, guidefontsize = 10,
            left_margin = 8Plots.mm, right_margin = 5Plots.mm,
            bottom_margin = 2Plots.mm,
            framestyle = :box
            )
## Export figures
savefig(fig2, joinpath("..", "article", "figures", "fig2.pdf"))


text:
import FrankWolfe
using LinearAlgebra
using Test

@testset "Simple benchmark_oracles function" begin
    n = Int(1e3)

    xpi = rand(n)
    total = sum(xpi)
    xp = xpi ./ total

    f(x) = norm(x - xp)^2
    function grad!(storage, x)
        @. storage = 2 * (x - xp)
    end

    lmo_prob = FrankWolfe.ProbabilitySimplexOracle(1)
    x0 = FrankWolfe.compute_extreme_point(lmo_prob, zeros(n))

    FrankWolfe.benchmark_oracles(f, grad!, () -> rand(n), lmo_prob; k=100)
end

@testset "RankOneMatrix" begin
    for n in (1, 2, 5)
        for _ in 1:5
            v = rand(n)
            u = randn(2n)
            M = u * v'
            R = FrankWolfe.RankOneMatrix(u, v)
            for i in 1:2n
                for j in 1:n
                    @test M[i, j] ≈ R[i, j]
                end
            end
            @testset "Right- left-mul" for _ in 1:5
                x = rand(n)
                r1 = R * x
                r2 = M * x
                @test r1 ≈ r2
            end
            @testset "Add and sub" begin
                @test M + R ≈ R + R
                @test M - R ≈ R - R
            end
        end
    end
end

@testset "Line Search methods" begin
    a = [-1.0, -1.0, -1.0]
    b = [1.0, 1.0, 1.0]
    function grad!(storage, x)
        return storage .= 2x
    end
    f(x) = norm(x)^2
    gradient = similar(a)
    grad!(gradient, a)
    @test FrankWolfe.backtrackingLS(f, gradient, a, a - b, 1.0) == (0.5, 1)
    @test abs(FrankWolfe.segment_search(f, grad!, a, a - b, 1.0)[1] - 0.5) < 0.0001

    @inferred FrankWolfe.line_search_wrapper(
        FrankWolfe.Agnostic(),
        3,
        identity,
        nothing,
        [0.3, 0.2],
        ones(2),
        ones(3),
        0.0,
        3,
        nothing,
        1e-6,
        1e-6,
        0.9,
    )
end


text:
#-----------------------------------------------------------------------------------------------------
# Translational Basis
#-----------------------------------------------------------------------------------------------------
struct TranslationalBasis{T <: Number} <: AbstractBasis
    dgt::Vector{Int}
    I::Vector{Int}
    R::Vector{Float64}
    C::T
    B::Int
end
eltype(::TranslationalBasis{T}) where T = T
#-----------------------------------------------------------------------------------------------------
function change!(b::TranslationalBasis, i::Integer)::Float64
    change!(b.dgt, b.I[i], base=b.B)
    1 / b.R[i]
end
#-----------------------------------------------------------------------------------------------------
function index(b::TranslationalBasis)::Tuple{ComplexF64, Int}
    Im, T = indmin(b.dgt, b.B)
    ind = binary_search(b.I, Im)
    if ind == 0
        return 0.0, 1
    else
        N = b.C ^ T * b.R[ind]
        return N, ind
    end
end
#-----------------------------------------------------------------------------------------------------
size(b::TranslationalBasis, i::Integer) = (i == 2 || i == 1) ? length(b.I) : 1
size(b::TranslationalBasis) = (l=length(b.I); (l,l))
#-----------------------------------------------------------------------------------------------------
# Construct
#-----------------------------------------------------------------------------------------------------
export translationalbasis
function translationalbasis(k::Integer, L::Integer; base::Integer=2)
    dgt = zeros(Int, L)
    I, R = Int[], Float64[]
    for i = 1:base^L
        change!(dgt, i, base=base)
        c, r = checkstate(dgt, i, base)
        if c && (k * r % L == 0)
            append!(I, i)
            append!(R, L/sqrt(r))
        end
    end
    expk = exp(-1im*2π/L*k)
    C = imag(expk) == 0 ? real(expk) : expk
    TranslationalBasis(dgt, I, R, C, base)
end
#-----------------------------------------------------------------------------------------------------
function translationalbasis(f, k::Integer, L::Integer; base::Integer=2)
    dgt = zeros(Int, L)
    I, R = Int[], Float64[]
    for i = 1:base^L
        change!(dgt, i, base=base)
        if f(dgt)
            c, r = checkstate(dgt, i, base)
            if c && (k * r % L == 0)
                append!(I, i)
                append!(R, L/sqrt(r))
            end
        end
    end
    expk = exp(-1im*2π/L*k)
    C = imag(expk) == 0 ? real(expk) : expk
    TranslationalBasis(dgt, I, R, C, base)
end

#-----------------------------------------------------------------------------------------------------
# Helper Functions
#-----------------------------------------------------------------------------------------------------
function movebits!(dgt::AbstractVector{<:Integer})
    p = dgt[end]
    for i = 1:length(dgt)
        p, dgt[i] = dgt[i], p
    end
    dgt
end
#-----------------------------------------------------------------------------------------------------
function checkstate(dgt::AbstractVector{<:Integer}, I0::Integer, base::Integer)
    for i=1:length(dgt)-1
        movebits!(dgt)
        In = index(dgt, base=base)
        if In < I0
            change!(dgt, I0, base=base)
            return false, 0
        elseif In == I0
            return true, i
        end
    end
    movebits!(dgt)
    true, length(dgt)
end
#-----------------------------------------------------------------------------------------------------
function indmin(dgt::AbstractVector{<:Integer}, base::Integer)
    I0 = index(dgt, base=base)
    Im, T = I0, 0
    for i=1:length(dgt)-1
        movebits!(dgt)
        In = index(dgt, base=base)
        if In == Im
            change!(dgt, I0, base=base)
            return Im, T
        elseif In < Im
            Im, T = In, i
        end
    end
    movebits!(dgt)
    Im, T
end


text:
using TrainView

using CSV, DataFrames
# using Statistics, LinearAlgebra

function track_to_cabin_movement(input_file,output_file)

    !isfile(input_file) && @error "there is no file: $input_file "

    println("Loading data file from:", input_file)
    println("I hope it's in the right format!")

    lines = open(input_file) do f
        [l for l in eachline(f)]
    end;

    header = lines[1];
    lines = lines[2:end];

    max_v = 720
    max_u = 1280

    trackprop = TrackProperties(track_gauge = 1.435 + 0.065)

    # from a previous calibration exercise
    camera_xyz = [0.0, 0.6771, -2.1754]
    camera_ψθφ = [0.0, -0.2038, -0.08219]

    camera = VideoCamera(camera_xyz;
        ψθφ = camera_ψθφ,
        focal_length = 5.8e-3,
        pixelspermeter = 1 / 5.5e-6
    )

    uv_data = map(lines) do l
        data = parse.(Int,split(l,',')[2:end-1])

        Llen = data[1]
        Rlen = data[2]

        us = data[3:2:end] .- (max_u / 2.0);
        vs = data[4:2:end] .- (max_v / 2.0);

        Lu = us[1:Llen];
        Lv = vs[1:Llen];

        Ru = us[Llen+1:Rlen+Llen];
        Rv = vs[Llen+1:Rlen+Llen];

        return [Lu,Lv,Ru,Rv]
    end;

    # Calculate train car movement from frames
        choose_distortions = [:Y,:Z,:θT,:φT,:α,:β];

        distortions = map(uv_data) do uv

            left_uvs = map(uv[1],uv[2]) do u, v
                [u,v]
            end

            right_uvs = map(uv[3],uv[4]) do u, v
                [u,v]
            end

            distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera, trackprop;
                choose_distortions = choose_distortions,
                iterations = 1)

            return [distortion[k] for k in choose_distortions]
        end;

    distortion_matrix = transpose(hcat(distortions...))

    df = DataFrame(distortion_matrix,choose_distortions)
    rename!(df,choose_distortions)

    CSV.write(output_file,df)

    println("Data saved as a CSV file in $(output_file)")

end

# Code below ensures that the above can be run from the terminal
if length(ARGS) > 1 && typeof(ARGS[1]) == String && typeof(ARGS[2]) == String
    track_to_cabin_movement(ARGS[1],ARGS[2])
end


text:
using ApproxFun, Test
    import ApproxFun: Block, BlockBandedMatrix

## PDEs

S = JacobiWeight(1.,1.,Jacobi(1.,1.))^2
Δ = Laplacian(S)

f = Fun((x,y)->sin(π*x)*sin(π*y),S)

QR1=qr(Δ)
    ApproxFun.resizedata!(QR1,:,400)
    \(QR1,f; tolerance=1E-10)
QR1=qr(Δ)
    @time Δ[Block.(1:40), Block.(1:40)]
    @time ApproxFun.resizedata!(QR1,:,400)
    @time \(QR1,f; tolerance=1E-10)
println("Laplace Dirichlet: should be ~0.015, 0.015, 0.001")

d=ChebyshevInterval()^2
A=[Dirichlet(d); Laplacian(d)]
f=Fun((x,y)->real(exp(x+im*y)),∂(d))

QR=qr(A)
    ApproxFun.resizedata!(QR,:,150)
    \(QR,[f; 0.];tolerance=1E-10)

QR=qr(A)
    @time ApproxFun.resizedata!(QR,:,150)
    @time u=\(QR,[f; 0.];tolerance=1E-10)

println("Laplace: should be ~0.06, 0.001")

d=ChebyshevInterval()^2
u = [Neumann(d); Laplacian(d)+100I] \ [[[1,1],[1,1]],0]
@time u = [Neumann(d); Laplacian(d)+100I] \ [[[1,1],[1,1]],0]
println("Neumann Helmholtz: should be ~0.032")

# dx=Interval(0.,1.);dt=Interval(0.0,0.54)
# d=dx*dt
#
# x,y=Fun(d)
# V=x^2
#
# Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])
#
# ϵ=0.1
#
# u0=Fun(x->exp(-25*(x-.5)^2)*exp(-1.0im/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
# L=1im*ϵ*Dt+.5*ϵ^2*Dx^2-V⊗1
#
# PO=discretize([timedirichlet(d);L],50)
# @time PO=discretize([timedirichlet(d);L],50)
# u=PO\u0
# @time    u=PO\u0
#
# println("Schrodinger: should be ~0.013,0.015")


text:
mutable struct Mesh
    # Original untransformed vertices
    vertices::Array{Point{Float64},1}
    # Transformed vertices
    bucket::Array{Point{Float64},1}

    function Mesh()
        obj = new()

        obj.vertices = Array{Point{Float64},1}[]
        obj.bucket = Array{Point{Float64},1}[]

        obj
    end
end

function add_vertex!(mesh::Mesh, x::T, y::T) where {T <: AbstractFloat}
    push!(mesh.vertices, Point{Float64}(x, y))
end

function build!(mesh::Mesh)
    for p in mesh.vertices
        push!(mesh.bucket, Point{Float64}())
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function sqlite3_errmsg()
    return ccall( (:sqlite3_errmsg, sqlite3_lib),
        Ptr{UInt8}, ()
        )
end
function sqlite3_errmsg(db::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_errmsg, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},),
        db)
end
function sqlite3_open(file::AbstractString,handle)
    return ccall( (:sqlite3_open, sqlite3_lib),
        Cint, (Ptr{UInt8},Ptr{Void}),
        file,handle)
end
function sqlite3_open16(file::UTF16String,handle)
    return ccall( (:sqlite3_open16, sqlite3_lib),
        Cint, (Ptr{UInt16},Ptr{Void}),
        file,handle)
end
function sqlite3_close(handle::Ptr{Void})
    @NULLCHECK handle
    return ccall( (:sqlite3_close, sqlite3_lib),
        Cint, (Ptr{Void},),
        handle)
end
function sqlite3_next_stmt(db::Ptr{Void},stmt::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_next_stmt, sqlite3_lib),
        Ptr{Void}, (Ptr{Void},Ptr{Void}),
        db, stmt)
end
function sqlite3_prepare_v2(handle::Ptr{Void},query::AbstractString,stmt,unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare_v2, sqlite3_lib),
        Cint, (Ptr{Void},Ptr{UInt8},Cint,Ptr{Void},Ptr{Void}),
            handle,query,sizeof(query),stmt,unused)
end
function sqlite3_prepare16_v2(handle::Ptr{Void},query::AbstractString,stmt,unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare16_v2, sqlite3_lib),
        Cint, (Ptr{Void},Ptr{UInt16},Cint,Ptr{Void},Ptr{Void}),
        handle,query,sizeof(query),stmt,unused)
end
function sqlite3_finalize(stmt::Ptr{Void})
    @NULLCHECK stmt
    return ccall( (:sqlite3_finalize, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end

# SQLITE_API int sqlite3_bind_paramter_count(sqlite3_stmt*)
function sqlite3_bind_parameter_count(stmt::Ptr{Void})
    return ccall( (:sqlite3_bind_parameter_count, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end
#SQLITE_API const char* sqlite3_bind_parameter_name(sqlite3_stmt*, int)
function sqlite3_bind_parameter_name(stmt::Ptr{Void}, col::Int)
    return ccall( (:sqlite3_bind_parameter_name, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void}, Cint),
        stmt, col)
end
# SQLITE_API int sqlite3_bind_parameter_index(sqlite3_stmt*, const char *zName);
function sqlite3_bind_parameter_index(stmt::Ptr{Void},value::AbstractString)
    return ccall( (:sqlite3_bind_parameter_index, sqlite3_lib),
        Cint, (Ptr{Void},Ptr{UInt8}),
        stmt,value)
end
# SQLITE_API int sqlite3_bind_double(sqlite3_stmt*, int, double);
function sqlite3_bind_double(stmt::Ptr{Void},col::Int,value::Float64)
    return ccall( (:sqlite3_bind_double, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Float64),
        stmt,col,value)
end
# SQLITE_API int sqlite3_bind_int(sqlite3_stmt*, int, int);
function sqlite3_bind_int(stmt::Ptr{Void},col::Int,value::Int32)
    return ccall( (:sqlite3_bind_int, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Int32),
        stmt,col,value)
end
# SQLITE_API int sqlite3_bind_int64(sqlite3_stmt*, int, sqlite3_int64);
function sqlite3_bind_int64(stmt::Ptr{Void},col::Int,value::Int64)
    return ccall( (:sqlite3_bind_int64, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Int64),
        stmt,col,value)
end
# SQLITE_API int sqlite3_bind_null(sqlite3_stmt*, int);
function sqlite3_bind_null(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_bind_null, sqlite3_lib),
        Cint, (Ptr{Void},Cint),
        stmt,col)
end
# SQLITE_API int sqlite3_bind_text(sqlite3_stmt*, int, const char*, int n, void(*)(void*));
function sqlite3_bind_text(stmt::Ptr{Void},col::Int,value::AbstractString)
    return ccall( (:sqlite3_bind_text, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Ptr{UInt8},Cint,Ptr{Void}),
        stmt,col,value,sizeof(value),C_NULL)
end
function sqlite3_bind_text(stmt::Ptr{Void},col::Int,ptr::Ptr{UInt8},len::Int)
    return ccall( (:sqlite3_bind_text, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Ptr{UInt8},Cint,Ptr{Void}),
        stmt,col,ptr,len,C_NULL)
end
# SQLITE_API int sqlite3_bind_text16(sqlite3_stmt*, int, const void*, int, void(*)(void*));
function sqlite3_bind_text16(stmt::Ptr{Void},col::Int,value::UTF16String)
    return ccall( (:sqlite3_bind_text16, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Ptr{UInt16},Cint,Ptr{Void}),
        stmt,col,value,sizeof(value),C_NULL)
end
function sqlite3_bind_text16(stmt::Ptr{Void},col::Int,ptr::Ptr{UInt16},len::Int)
    return ccall( (:sqlite3_bind_text16, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Ptr{UInt16},Cint,Ptr{Void}),
        stmt,col,ptr,len,C_NULL)
end

# SQLITE_API int sqlite3_bind_blob(sqlite3_stmt*, int, const void*, int n, void(*)(void*));
function sqlite3_bind_blob(stmt::Ptr{Void},col::Int,value)
    return ccall( (:sqlite3_bind_blob, sqlite3_lib),
        Cint, (Ptr{Void},Cint,Ptr{UInt8},Cint,Ptr{Void}),
        stmt,col,value,sizeof(value),SQLITE_STATIC)
end
# SQLITE_API int sqlite3_bind_zeroblob(sqlite3_stmt*, int, int n);
# SQLITE_API int sqlite3_bind_value(sqlite3_stmt*, int, const sqlite3_value*);

function sqlite3_clear_bindings(stmt::Ptr{Void})
    return ccall( (:sqlite3_clear_bindings, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end

function sqlite3_step(stmt::Ptr{Void})
    return ccall( (:sqlite3_step, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end
function sqlite3_column_count(stmt::Ptr{Void})
    return ccall( (:sqlite3_column_count, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end
function sqlite3_column_type(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_type, sqlite3_lib),
        Cint, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_blob(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_blob, sqlite3_lib),
        Ptr{Void}, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_bytes(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_bytes, sqlite3_lib),
        Cint, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_bytes16(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_bytes16, sqlite3_lib),
        Cint, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_double(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_double, sqlite3_lib),
        Cdouble, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_int(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_int, sqlite3_lib),
        Cint, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_int64(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_int64, sqlite3_lib),
        Clonglong, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_text(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_text, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_text16(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_text16, sqlite3_lib),
        Ptr{Void}, (Ptr{Void},Cint),
        stmt,col-1)
end
# function sqlite3_column_value(stmt::Ptr{Void},col::Cint)
#     return ccall( (:sqlite3_column_value, sqlite3_lib),
#             Ptr{Void}, (Ptr{Void},Cint),
#             stmt,col-1)
# end
# SQLITE_API sqlite3_value *sqlite3_column_value(sqlite3_stmt*, int iCol);
function sqlite3_reset(stmt::Ptr{Void})
    return ccall( (:sqlite3_reset, sqlite3_lib),
        Cint, (Ptr{Void},),
        stmt)
end

# SQLITE_API const char *sqlite3_column_name(sqlite3_stmt*, int N);
function sqlite3_column_name(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_name, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},Cint),
        stmt,col-1)
end
function sqlite3_column_name16(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_name16, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},Cint),
        stmt,col-1)
end

function sqlite3_changes(db::Ptr{Void})
    @NULLCHECK db
   return ccall( (:sqlite3_changes, sqlite3_lib),
        Cint, (Ptr{Void},),
        db)
end
function sqlite3_total_changes(db::Ptr{Void})
    @NULLCHECK db
   return ccall( (:sqlite3_changes, sqlite3_lib),
        Cint, (Ptr{Void},),
        db)
end
# SQLITE_API const void *sqlite3_column_name16(sqlite3_stmt*, int N);

# SQLITE_API const char *sqlite3_column_database_name(sqlite3_stmt*,int);
# SQLITE_API const void *sqlite3_column_database_name16(sqlite3_stmt*,int);
# SQLITE_API const char *sqlite3_column_table_name(sqlite3_stmt*,int);
# SQLITE_API const void *sqlite3_column_table_name16(sqlite3_stmt*,int);
# SQLITE_API const char *sqlite3_column_origin_name(sqlite3_stmt*,int);
# SQLITE_API const void *sqlite3_column_origin_name16(sqlite3_stmt*,int);

function sqlite3_column_decltype(stmt::Ptr{Void},col::Int)
    return ccall( (:sqlite3_column_decltype, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},Cint),
        stmt,col-1)
end
# SQLITE_API const char *sqlite3_column_decltype(sqlite3_stmt*,int);
# SQLITE_API const void *sqlite3_column_decltype16(sqlite3_stmt*,int);

# SQLITE_API int sqlite3_data_count(sqlite3_stmt *pStmt);

# SQLITE_API void sqlite3_result_double(sqlite3_context*, double);
function sqlite3_result_double(context::Ptr{Void},value::Float64)
    return ccall( (:sqlite3_result_double, sqlite3_lib),
        Void, (Ptr{Void},Float64),
        context,value)
end
# SQLITE_API void sqlite3_result_error(sqlite3_context*, const char*, int)
function sqlite3_result_error(context::Ptr{Void},msg::AbstractString)
    return ccall( (:sqlite3_result_error, sqlite3_lib),
        Void, (Ptr{Void},Ptr{UInt8},Cint),
        context,value,sizeof(msg)+1)
end
# SQLITE_API void sqlite3_result_error16(sqlite3_context*, const void*, int)
function sqlite3_result_error(context::Ptr{Void},msg::UTF16String)
    return ccall( (:sqlite3_result_error16, sqlite3_lib),
        Void, (Ptr{Void},Ptr{UInt16},Cint),
        context,value,sizeof(msg)+1)
end
# SQLITE_API void sqlite3_result_int(sqlite3_context*, int);
function sqlite3_result_int(context::Ptr{Void},value::Int32)
    return ccall( (:sqlite3_result_int, sqlite3_lib),
        Void, (Ptr{Void},Int32),
        context,value)
end
# SQLITE_API void sqlite3_result_int64(sqlite3_context*, sqlite3_int64);
function sqlite3_result_int64(context::Ptr{Void},value::Int64)
    return ccall( (:sqlite3_result_int64, sqlite3_lib),
        Void, (Ptr{Void},Int64),
        context,value)
end
# SQLITE_API void sqlite3_result_null(sqlite3_context*);
function sqlite3_result_null(context::Ptr{Void})
    return ccall( (:sqlite3_result_null, sqlite3_lib),
        Void, (Ptr{Void},),
        context)
end
# SQLITE_API void sqlite3_result_text(sqlite3_context*, const char*, int n, void(*)(void*));
function sqlite3_result_text(context::Ptr{Void},value::AbstractString)
    return ccall( (:sqlite3_result_text, sqlite3_lib),
        Void, (Ptr{Void},Ptr{UInt8},Cint,Ptr{Void}),
        context,value,sizeof(value)+1,SQLITE_TRANSIENT)
end
# SQLITE_API void sqlite3_result_text16(sqlite3_context*, const void*, int, void(*)(void*));
function sqlite3_result_text16(context::Ptr{Void},value::UTF16String)
    return ccall( (:sqlite3_result_text, sqlite3_lib),
        Void, (Ptr{Void},Ptr{UInt16},Cint,Ptr{Void}),
        context,value,sizeof(value)+1,SQLITE_TRANSIENT)
end
# SQLITE_API void sqlite3_result_blob(sqlite3_context*, const void*, int n, void(*)(void*));
function sqlite3_result_blob(context::Ptr{Void},value)
    return ccall( (:sqlite3_result_blob, sqlite3_lib),
        Void, (Ptr{Void},Ptr{UInt8},Cint,Ptr{Void}),
        context,value,sizeof(value),SQLITE_TRANSIENT)
end
# SQLITE_API void sqlite3_result_zeroblob(sqlite3_context*, int n);
# SQLITE_API void sqlite3_result_value(sqlite3_context*, const sqlite3_value*);
# SQLITE_API void sqlite3_result_error_toobig(sqlite3_context*)
# SQLITE_API void sqlite3_result_error_nomem(sqlite3_context*)
# SQLITE_API void sqlite3_result_error_code(sqlite3_context*, int)

function sqlite3_create_function_v2(db::Ptr{Void},name::AbstractString,nargs::Integer,
                                    enc::Integer,data::Ptr{Void},func::Ptr{Void},
                                    step::Ptr{Void},final::Ptr{Void},
                                    destructor::Ptr{Void})
    @NULLCHECK db
    return ccall(
        (:sqlite3_create_function_v2, sqlite3_lib),
        Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Void},
         Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
        db, name, nargs, enc, data, func, step, final, destructor)
end

# SQLITE_API void* sqlite3_aggregate_context(sqlite3_context*, int nBytes)
function sqlite3_aggregate_context(context::Ptr{Void}, nbytes::Integer)
    return ccall( (:sqlite3_aggregate_context, sqlite3_lib),
        Ptr{Void}, (Ptr{Void}, Cint),
        context, nbytes)
end

# SQLITE_API int sqlite3_value_type(sqlite3_value*)
function sqlite3_value_type(value::Ptr{Void})
    return ccall( (:sqlite3_value_type, sqlite3_lib),
        Cint, (Ptr{Void},),
        value)
end

# SQLITE_API const void* sqlite3_value_blob(sqlite3_value*)
function sqlite3_value_blob(value::Ptr{Void})
    return ccall( (:sqlite3_value_blob, sqlite3_lib),
        Ptr{Void}, (Ptr{Void},),
        value)
end
# SQLITE_API int sqlite3_value_bytes(sqlite3_value*)
function sqlite3_value_bytes(value::Ptr{Void})
    return ccall( (:sqlite3_value_bytes, sqlite3_lib),
        Cint, (Ptr{Void},),
        value)
end
# SQLITE_API int sqlite3_value_bytes16(sqlite3_value*)
function sqlite3_value_bytes16(value::Ptr{Void})
    return ccall( (:sqlite3_value_bytes16, sqlite3_lib),
        Cint, (Ptr{Void},),
        value)
end
# SQLITE_API double sqlite3_value_double(sqlite3_value*)
function sqlite3_value_double(value::Ptr{Void})
    return ccall( (:sqlite3_value_double, sqlite3_lib),
        Cdouble, (Ptr{Void},),
        value)
end
# SQLITE_API int sqlite3_value_int(sqlite3_value*)
function sqlite3_value_int(value::Ptr{Void})
    return ccall( (:sqlite3_value_int, sqlite3_lib),
        Cint, (Ptr{Void},),
        value)
end
# SQLITE_API sqlite_int64 sqlite3_value_int64(sqlite3_value*)
function sqlite3_value_int64(value::Ptr{Void})
    return ccall( (:sqlite3_value_int64, sqlite3_lib),
        Clonglong, (Ptr{Void},),
        value)
end
# SQLITE_API const unsigned char* sqlite3_value_text(sqlite3_value*)
function sqlite3_value_text(value::Ptr{Void})
    return ccall( (:sqlite3_value_text, sqlite3_lib),
        Ptr{UInt8}, (Ptr{Void},),
        value)
end
# SQLITE_API const void* sqlite3_value_text16(sqlite3_value*)
function sqlite3_value_text16(value::Ptr{Void})
    return ccall( (:sqlite3_value_text16, sqlite3_lib),
        Ptr{Void}, (Ptr{Void},),
        value)
end
# SQLITE_API int sqlite3_value_numeric_type(sqlite3_value*)

function sqlite3_initialize()
    return ccall( (:sqlite3_initialize, sqlite3_lib),
        Cint, (),
        )
end
function sqlite3_shutdown()
    return ccall( (:sqlite3_shutdown, sqlite3_lib),
        Cint, (),
        )
end
function sqlite3_os_init()
    return ccall( (:sqlite3_os_init, sqlite3_lib),
        Cint, (),
        )
end
function sqlite3_os_end()
    return ccall( (:sqlite3_os_end, sqlite3_lib),
        Cint, (),
        )
end
function sqlite3_free_table(result::Array{AbstractString,1})
    return ccall( (:sqlite3_free_table, sqlite3_lib),
        Void, (Ptr{Ptr{Void}},),
        result)
end

# SQLITE_API const char *sqlite3_uri_parameter(const char *zFilename, const char *zParam);
# SQLITE_API int sqlite3_uri_boolean(const char *zFile, const char *zParam, int bDefault);
# SQLITE_API sqlite3_int64 sqlite3_uri_int64(const char*, const char*, sqlite3_int64);
function sqlite3_errcode(db::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_errcode, sqlite3_lib),
        Cint, (Ptr{Void},),
        db)
end
function sqlite3_extended_errcode(db::Ptr{Void})
    @NULLCHECK db
    return ccall( (:sqlite3_extended_errcode, sqlite3_lib),
        Cint, (Ptr{Void},),
        db)
end
# SQLITE_API int sqlite3_errcode(sqlite3 *db);
# SQLITE_API int sqlite3_extended_errcode(sqlite3 *db);
function sqlite3_errstr(ret::Cint)
    return ccall( (:sqlite3_errstr, sqlite3_lib),
        Ptr{UInt8}, (Cint,),
        ret)
end
# SQLITE_API const char *sqlite3_errstr(int);

# SQLITE_API int sqlite3_limit(sqlite3*, int id, int newVal);

# SQLITE_API int sqlite3_stmt_readonly(sqlite3_stmt *pStmt);

# SQLITE_API int sqlite3_stmt_busy(sqlite3_stmt*);

# SQLITE_API int sqlite3_table_column_metadata(
#   sqlite3 *db,                /* Connection handle */
#   const char *zDbName,        /* Database name or NULL */
#   const char *zTableName,     /* Table name */
#   const char *zColumnName,    /* Column name */
#   char const **pzDataType,    /* OUTPUT: Declared data type */
#   char const **pzCollSeq,     /* OUTPUT: Collation sequence name */
#   int *pNotNull,              /* OUTPUT: True if NOT NULL constraint exists */
#   int *pPrimaryKey,           /* OUTPUT: True if column part of PK */
#   int *pAutoinc               /* OUTPUT: True if column is auto-increment */
# );

# SQLITE_API int sqlite3_db_status(sqlite3*, int op, int *pCur, int *pHiwtr, int resetFlg);

# Not directly used
function sqlite3_open_v2(file::AbstractString,handle,flags::Cint,vfs::AbstractString)
    return ccall( (:sqlite3_open_v2, sqlite3_lib),
            Cint, (Ptr{UInt8},Ptr{Void},Cint,Ptr{UInt8}),
            file,handle,flags,vfs)
end
function sqlite3_prepare(handle::Ptr{Void},query::AbstractString,stmt,unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare, sqlite3_lib),
        Cint, (Ptr{Void},Ptr{UInt8},Cint,Ptr{Void},Ptr{Void}),
            handle,query,sizeof(query),stmt,unused)
end
function sqlite3_prepare16(handle::Ptr{Void},query::AbstractString,stmt,unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare16, sqlite3_lib),
        Cint, (Ptr{Void},Ptr{UInt8},Cint,Ptr{Void},Ptr{Void}),
            handle,query,sizeof(query),stmt,unused)
end
function sqlite3_close_v2(handle::Ptr{Void})
    @NULLCHECK handle
    try
        return ccall( (:sqlite3_close_v2, sqlite3_lib),
            Cint, (Ptr{Void},),
            handle)
    catch
        # Older versions of the library don't have this, abort to other close
        warn("sqlite3_close_v2 not available.")
        sqlite3_close(handle)
    end
end

max_stars_repo_name: beffaCo/updateSqlite

max_issues_repo_name: beffaCo/updateSqlite

max_forks_repo_name: beffaCo/updateSqlite


text:
"""
The module `PDB` defines types and methods to work with protein structures inside Julia.
It is useful to link structural and sequential information, and
needed for measure the predictive performance at protein contact prediction of mutual information scores.

**Features**

- Read and parse PDF and PDBML files
- Calculate distance and contacts between atoms or residues
- Determine interaction between residues

```julia
using MIToS.PDB
```
"""
module PDB

using RecipesBase           # Plots for PDB Residues
using AutoHashEquals
using StaticArrays
using DataStructures
using PairwiseListMatrices
using NamedArrays
using LightXML
using LinearAlgebra
using Statistics            # mean
using MIToS.Utils
using Formatting

export  # PDBResidues
        PDBResidueIdentifier,
        Coordinates,
        PDBAtom,
        PDBResidue,
        squared_distance,
        distance,
        contact,
        isresidue,
        isatom,
        residues,
        @residues,
        residuesdict,
        @residuesdict,
        atoms,
        @atoms,
        findheavy,
        findatoms,
        findCB,
        selectbestoccupancy,
        bestoccupancy,
        residuepairsmatrix,
        proximitymean,
        # AtomsData
        covalentradius,
        vanderwaalsradius,
        check_atoms_for_interactions,
        # Interaction
        ishydrophobic,
        isaromatic,
        iscationic,
        isanionic,
        ishbonddonor,
        ishbondacceptor,
        hydrogenbond,
        vanderwaals,
        vanderwaalsclash,
        covalent,
        disulphide,
        aromaticsulphur,
        pication,
        aromatic,
        ionic,
        hydrophobic,
        stridehydrogenbond,
        chimerahydrogenbond,
        # PDBParser
        PDBFile,
        # PDBMLParser
        PDBML,
        downloadpdb,
        PDBMLHeader,
        downloadpdbheader,
        getpdbdescription,
        # Kabsch
        kabsch,
        center!,
        rmsd,
        getCA,
        CAmatrix,
        coordinatesmatrix,
        change_coordinates,
        centeredcoordinates,
        centeredresidues,
        superimpose,
        mean_coordinates,
        rmsf,
        # MIToS.Utils
        All,
        # Imported from Base (and exported for docs)
        any,
        parse,
        print,
        angle

include("PDBResidues.jl")
include("AtomsData.jl")
include("Interaction.jl")
include("PDBParser.jl")
include("PDBMLParser.jl")
include("Kabsch.jl")
include("Plots.jl")

end


text:
function expectation(g, distr::ContinuousUnivariateDistribution; kwargs...)
    return first(quadgk(x -> pdf(distr, x) * g(x), extrema(distr)...; kwargs...))
end

## Assuming that discrete distributions only take integer values.
function expectation(g, distr::DiscreteUnivariateDistribution; epsilon::Real=1e-10)
    mindist, maxdist = extrema(distr)
    # We want to avoid taking values up to infinity
    minval = isfinite(mindist) ? mindist : quantile(distr, epsilon)
    maxval = isfinite(maxdist) ? maxdist : quantile(distr, 1 - epsilon)
    return sum(x -> pdf(distr, x) * g(x), minval:maxval)
end

function expectation(g, distr::MultivariateDistribution; nsamples::Int=100, rng::AbstractRNG=GLOBAL_RNG)
    nsamples > 0 || throw(ArgumentError("number of samples should be > 0"))
    # We use a function barrier to work around type instability of `sampler(dist)`
    return mcexpectation(rng, g, sampler(distr), nsamples)
end

mcexpectation(rng, f, sampler, n) = sum(f, rand(rng, sampler) for _ in 1:n) / n

## Leave undefined until we've implemented a numerical integration procedure
# function entropy(distr::UnivariateDistribution)
#     pf = typeof(distr)<:ContinuousDistribution ? pdf : pmf
#     f = x -> pf(distr, x)
#     expectation(distr, x -> -log(f(x)))
# end

function kldivergence(P::Distribution{V}, Q::Distribution{V}; kwargs...) where {V<:VariateForm}
    return expectation(P; kwargs...) do x
        logp = logpdf(P, x)
        return (logp > oftype(logp, -Inf)) * (logp - logpdf(Q, x))
    end
end


text:
# Left division
@inline \(mat::AbstractMatrix, tvec::TransposedVector) = error("Cannot left-divide transposed vector by matrix")

# These four are trivial definitions in Base, resorting to transposition:
#   A_ldiv_Bc
#   A_ldiv_Bt
#   At_ldiv_Bt
#   Ac_ldiv_Bc

# The others need to be overloaded...
@inline At_ldiv_B(mat::AbstractMatrix, tvec::TransposedVector) = error("Cannot left-divide transposed vector by matrix")
@inline Ac_ldiv_B(mat::AbstractMatrix, tvec::TransposedVector) = error("Cannot left-divide transposed vector by matrix")


text:
"""
Complex quantities represent position, velocity etc in a plane.
They are a subset of quantities, although some function methods may
need extending outside T<:Real

"""

# TODO improve this. We want to use promotion to define tuples with
# identical types.
# It might also be nice to use traits, so that we can dispatch on
# 2d quantities.
#=
 dimensionality(T<:Real) = Val{1}()
....etc  d

There are propably examples in Julia base, for n-dimensional arrays.

 foo(x::T) where T = foo(dimensionality(T), x)
 foo(::Val{2}, x) =

 angle(::Val{2}, x::{Quantity{<:Complex} })
 =#

const ComplexQuantity = Quantity{<:Complex}

const RealQuantity = Quantity{<:Real}
const VelocityTuple = Tuple{Velocity, Velocity}
const PositionTuple = Tuple{Length, Length}
const ForceTuple = Tuple{Force, Force}

const QuantityTuple = Tuple{Quantity, Quantity}
ComplexQuantity(p::QuantityTuple) = complex(p[1] , p[2])
QuantityTuple(z::ComplexQuantity) =  reim(z)

function string_polar_form(z::ComplexQuantity)
    strarg_r = string(round(unit(hypot(z)), hypot(z), digits = 3))
    strarg_θ = string(round(angle(z), digits = 3))
    strargument = strarg_r * "∙exp(" * strarg_θ * "im)"
end


text:
export rnn

function rnn(layersNames,layersDims,T,netType=:regression)
  L=length(layersNames)

  global order,sharedParams
  orderGraph=rnnOrderGraph(L,T)
  connectGraph=rnnConnectGraph(L,T)
  order=topoSort(orderGraph)

  layersNames=repmat(layersNames,T,1)
  layersDims=repmat(layersDims,T,1)
  layersDims=correctLayerDims(orderGraph,layersDims)

  net=initNet(layersNames,layersDims,:regression)
  connectNet(net,connectGraph,layersDims)

  for i=1:length(net)
    if i%L==0
      net[i-1].A=linear;net[i-1].∂A=∂linear
    end
  end

  sharedParams=Array{Any}(L-2,2)
  for l=2:L-1
    sharedParams[l-1,1]=collect(l+(1:T-1)*L)
    sharedParams[l-1,2]=[:W,:b]
  end

  return net
end


text:
type DumbFoldAssigner <: FoldAssigner end
let
    @test_throws Exception assign_folds(DumbFoldAssigner(), 5, 2)

    assignment = assign_folds(RandomFoldAssigner(), 4, 2)
    @test length(find(assignment .== 1)) == 2
    @test length(find(assignment .== 2)) == 2

    assignment = assign_folds(RandomFoldAssigner(true), 5, 2)
    @test length(find(assignment .== 1)) == 3
    @test length(find(assignment .== 2)) == 2

    assignment = assign_folds(RandomFoldAssigner(false), 5, 2)
    @test length(find(assignment .== 1)) == 2
    @test length(find(assignment .== 2)) == 2

    @test assign_folds(OrderedFoldAssigner(), 4, 2) == [1, 1, 2, 2]
    @test assign_folds(OrderedFoldAssigner(true), 4, 3) == [1, 1, 2, 3]
    @test assign_folds(OrderedFoldAssigner(true), 5, 2) == [1, 1, 1, 2, 2]
    @test assign_folds(OrderedFoldAssigner(false), 5, 2) == [1, 1, 2, 2, 0]
end


text:
function resolve_ref(x, m::File, state::State, visited_scopes = 0)
    hasref(x) && return true
    return false
end

function resolve_ref(x, m::Nothing, state::State, visited_scopes = 0)
    hasref(x) && return true
    return false
end
function resolve_ref(x, m::T, state::State, visited_scopes = 0) where T
    hasref(x) && return true
    @warn "unhandled $T"
    return false
end
function resolve_ref(x1, m::SymbolServer.ModuleStore, state::State, visited_scopes = 0)
    hasref(x1) && return true
    if isidentifier(x1)
        x = x1
        if valof(x) == m.name
            setref!(x, m)
            return true
        elseif valof(x) in m.exported && haskey(m.vals, valof(x))
            val = m.vals[valof(x)]
            if val isa SymbolServer.PackageRef
                val1 = SymbolServer._lookup(val, getsymbolserver(state.server))
                if val1 !== nothing
                    setref!(x, val1)
                    return true
                else
                    return false
                end
            else
                setref!(x, val)
                return true
            end
        end
    elseif typof(x1) === MacroName
        x = x1.args[2]
        mn = string("@", valof(x))
        if mn in m.exported
            setref!(x, m.vals[mn])
            return true
        end
    elseif typof(x1) === x_Str
        mac = x1
        mn = string("@", valof(mac.args[1]), "_str")
        if mn in m.exported && haskey(m.vals, mn)
            setref!(mac.args[1], m.vals[mn])
            return true
        end
    end
    return false
end

function resolve_ref(x, scope::Scope, state::State, visited_scopes = 0)
    if visited_scopes > 50
        @info "Warning: circular reference found while resolving reference."
        return false
    end
    hasref(x) && return true
    resolved = false
    if typof(x) === BinaryOpCall && kindof(x.args[2]) === CSTParser.Tokens.DOT
        return resolve_getindex(x, scope, state)
    elseif isidentifier(x)
        mn = valof(x)
        x1 = x
        if (mn == "__source__" || mn == "__module__") && _in_macro_def(x)
            setref!(x, NoReference)
            return true
        end
    elseif typof(x) === MacroName
        x1 = x.args[2]
        mn = string("@", valof(x1))
    elseif typof(x) === x_Str
        if typof(x.args[1]) === IDENTIFIER
            x1 = x.args[1]
            mn = string("@", valof(x1), "_str")
        else
            return false
        end
    else
        return false
    end

    if haskey(scope.names, mn)
        setref!(x1, scope.names[mn])
        push!(scope.names[mn].refs, x1)
        resolved = true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for m in scope.modules
            resolved = resolve_ref(x, m[2], state, visited_scopes + 1)
            resolved && break
        end
    end
    if !hasref(x) && !scope.ismodule &&!(parentof(scope) isa EXPR)
        return resolve_ref(x, parentof(scope), state, visited_scopes + 1)
    end
    return resolved
end

function resolve_getindex(x::EXPR, scope::Scope, state::State)
    hasref(x) && return true
    resolved = false
    if typof(x.args[1]) === IDENTIFIER
        resolved = resolve_ref(x.args[1], scope, state)
        if resolved && typof(x.args[3]) === Quotenode && typof(x.args[3].args[1]) === IDENTIFIER
            resolved = resolve_getindex(x.args[3].args[1], refof(x.args[1]), state)
        end
    elseif typof(x.args[1]) === BinaryOpCall && kindof(x.args[1].args[2]) === CSTParser.Tokens.DOT
        resolved = resolve_ref(x.args[1], scope, state)
        if resolved && typof(x.args[3]) === Quotenode && typof(x.args[3].args[1]) === IDENTIFIER
            resolved = resolve_getindex(x.args[3].args[1], refof(x.args[1].args[3].args[1]), state)
        end
    end
    return resolved
end

function resolve_getindex(x::EXPR, b::Binding, state::State)
    hasref(x) && return true
    resolved = false
    if b.t isa Binding
        resolved = resolve_getindex(x, b.t.val, state)
    elseif b.val isa SymbolServer.ModuleStore
        resolved = resolve_getindex(x, b.val, state)
    elseif b.val isa EXPR && typof(b.val) === ModuleH
        resolved = resolve_getindex(x, b.val, state)
    elseif b.val isa Binding && b.val.val isa EXPR && typof(b.val.val) === ModuleH
        resolved = resolve_getindex(x, b.val.val, state)
    end
    return resolved
end

function resolve_getindex(x::EXPR, parent_type, state::State)
    hasref(x) && return true
    return false
end
function resolve_getindex(x::EXPR, parent_type::EXPR, state::State)
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x)
        if (typof(parent_type) === ModuleH || typof(parent_type) === BareModule) && scopeof(parent_type) isa Scope
            resolved = resolve_ref(x, scopeof(parent_type), state)
        elseif CSTParser.defines_struct(parent_type)
            if haskey(scopeof(parent_type).names, valof(x))
                setref!(x, scopeof(parent_type).names[valof(x)])
                push!(scopeof(parent_type).names[valof(x)].refs, x)
                # push!(bindingof(parent_type).refs, x)
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_getindex(x::EXPR, parent::SymbolServer.SymStore, state::State)
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x)
        if parent isa SymbolServer.ModuleStore && haskey(parent.vals, valof(x))
            val = parent.vals[valof(x)]
            if val isa SymbolServer.PackageRef
                val = SymbolServer._lookup(val, getsymbolserver(state.server))
                val === nothing && return false
            end
            setref!(x, val)
            resolved = true
        elseif parent isa SymbolServer.DataTypeStore && valof(x) in parent.fields
        end
    end
    return resolved
end

resolvable_macroname(x) = typof(x) === MacroName && length(x.args) == 2 && isidentifier(x.args[2]) && refof(x.args[2]) === nothing

function hasref(x::EXPR)
    refof(x) !== nothing && refof(x) !== NoReference
end

function _in_macro_def(x)
    if typof(x) === CSTParser.Macro
        return true
    elseif parentof(x) isa EXPR
        return _in_macro_def(parentof(x))
    else
        return false
    end
end


text:
mutable struct WilcoxTestResult<:TestResult
    type::Symbol  # either :paired or :independent
    is_corrected::Bool
    levels::Union{Vector{String}, Nothing}
    n_obs::IntVector
    means::FloatVector
    sds::FloatVector
    statistic::Float64
    pval::Float64
    cohen_d::Union{Float64, Nothing}
end

function show(io::IO, wtr::WilcoxTestResult)
    if wtr.type == :paired
        if wtr.is_corrected
            println(
                io,
                "\n        Wilcoxon signed rank sum test with continuity correction\n"
            )
        else
            println(io, "\n        Wilcoxon signed rank sum exact test\n")
        end
        @printf(io, "Difference (Mean ± SD): %.2f ± %.2f\n",
                wtr.means[1], wtr.sds[1])
    elseif wtr.type == :independent
        if wtr.is_corrected
            println(
                io,
                "\n        Wilcoxon rank sum test with continuity correction\n"
            )
        else
            println(io, "\n        Wilcoxon rank sum exact test\n")
        end
        println(io, "Data (Mean ± SD)")
        @printf(io, "  %s: %.2f ± %.2f\n",
                wtr.levels[1], wtr.means[1], wtr.sds[1])
        @printf(io, "  %s: %.2f ± %.2f\n",
                wtr.levels[2], wtr.means[2], wtr.sds[2])
    end
    println(io, "\nResult:\n", coeftable(wtr))
end

function coeftable(wtr::WilcoxTestResult)
    if (wtr.type == :paired)
        CoefTable(
            [wtr.statistic, wtr.pval],
            ["W", "P-value"],
            [""],
            2, 1
        )
    elseif (wtr.type == :independent)
        CoefTable(
            [wtr.statistic, wtr.pval, wtr.cohen_d],
            ["W", "P-value", "Cohen's d"],
            [""],
            2, 1
        )
    else
        error("invalid Wilcoxon test type $wtr.type")
    end
end

"""
        WilcoxDist(n1, n2)

Distribution of the Wilcox statistic
"""
struct WilcoxDist <: D.DiscreteUnivariateDistribution
    n1::Int
    n2::Int
    w::Array{Float64, 3}
end
D.@distr_support WilcoxDist 0 n1*n2
WilcoxDist(n1::Int, n2::Int) = WilcoxDist(
    n1,
    n2,
    fill(-1.0, (max(n1, n2, 50),
                max(n1, n2, 50),
                ceil(Int, max(n1,50) * max(n2,50) / 2)))
)
Base.show(io::IO, wd::WilcoxDist) = print(io, "Wilcox distribution with (", wd.n1, ", ", wd.n2, ")\n")

# NOTE: Julia starts from 1, not from 0
function cwilcox(w::Array{Float64, 3}, k::Int, m::Int, n::Int)::Float64
    mn = m * n
    if (k < 0) | (k > mn+1)
        return 0.0
    end
    c = floor(Int, mn / 2)
    if k > c
        k = mn - k  # hence k <= floor(mn/2)
    end
    i = min(m, n)
    j = max(m, n)  # hence i <= j
    if j == 0  # hence i == 0
        if k == 0
            return 1.0
        else
            return 0.0
        end
    end

    if (j > 0) & (k < j)
        return cwilcox(w, k, i, k)
    end
    if w[i,j,k] < 0
        w[i,j,k] = cwilcox(w, k-j, i-1, j) + cwilcox(w, k, i, j-1)
    end
    # @printf(stdout, "(i,j,k) = (%d,%d,%d) | c = %.2f\n", i, j, k, w[i,j,k])
    return w[i,j,k]
end

function cdf(d::WilcoxDist, x::Float64)
    n12 = d.n1 * d.n2
    c = binomial(big(d.n1 + d.n2), big(d.n2))
    if x < 0.0
        return 0.0
    elseif x ≥ n12
        return 1.0
    else
        p = 0.0
        if x ≤ (n12/2.0)
            for i = 0:floor(Int, x)
                p += cwilcox(d.w, i, d.n1, d.n2) / c
            end
        else
            x = n12 - x
            for i = 0:(floor(Int, x)-1)
                p += cwilcox(d.w, i, d.n1, d.n2) / c
            end
        end
        return p
    end
end

ccdf(d::WilcoxDist, x::Float64) = 1.0 - cdf(d, x)

"""
        WilcoxSignedDist(n1, n2)

Distribution of the Wilcox statistic
"""
struct WilcoxSignedDist <: D.DiscreteUnivariateDistribution
    n::Int
    w::Array{Float64}
end
D.@distr_support WilcoxSignedDist 0 n
WilcoxSignedDist(n::Int) = WilcoxSignedDist(
    n,
    zeros(ceil(Int, max(n,50)*(max(n,50)+1)/2))
)
Base.show(io::IO, wsd::WilcoxSignedDist) = print(io, "Wilcox signed distribution with ", wsd.n, "\n")

function cwilcoxsigned(w::Array{Float64}, k::Int, n::Int)::Float64
    u = n * (n+1) / 2.0
    c = floor(Int, u/2)

    if (k < 1) | (k > u)
        return 0.0
    end

    if k > c
        k = u - k
    end

    if n == 1
        return 1.0
    end

    if w[1] == 1.0
        return w[k]
    else
        w[1] = 1.0
        w[2] = 1.0
        for j = 3:(n+2)
            i = min(floor(Int, j*(j+1)/2), c)
            while i ≥ j
                w[i] += w[i-j+1]
                i -= 1
            end
        end
        return w[k]
    end
end

function cdf(d::WilcoxSignedDist, x::Float64)
    if x < 0.0
        return 0.0
    elseif x > (d.n*(d.n+1)/2.0)
        return 1.0
    else
        nf = float(d.n)
        f = exp(- nf * log(2.0))
        p = 0.0
        if x ≤ (d.n*(d.n+1)/4.0)
            for i = 1:(ceil(Int, x)+1)
                p += cwilcoxsigned(d.w, i, d.n) * f
            end
        else
            x = d.n*(d.n+1)/2.0 - x
            for i = 1:ceil(Int, x)
                p += cwilcoxsigned(d.w, i, d.n) * f
            end
        end
        return p
    end
end

ccdf(d::WilcoxSignedDist, x::Float64) = 1.0 - cdf(d, x)

"""
    wilcox_test(x1, x2, paired::Bool=false)

Wilcoxon (signed) rank sum test

The arguments `x1` and `x2` can be numeric `Vector`.

Wilcoxon signed rank sum test will be performed if the keyword argument `paired` is `true`. Otherwise, Wilcoxon rank sum test will be performed.
"""
function wilcox_test(
    x1::NumVector,
    x2::NumVector;
    levels::Union{Vector{String}, Nothing}=nothing,
    paired::Bool=false,
    exact::Union{Bool, Nothing}=nothing,
    μ0::AbstractFloat=0.0
)::WilcoxTestResult
    # use temporal sample names
    if !paired & isnothing(levels)
        @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
        levels = ["Sample1", "Sample2"]
    end

    # check types
    if !(x1 isa FloatVector)
        @warn "$levels[1] is converted to Float64"
        x1 = convert(Vector{Float64}, x1)
    end
    if !(x2 isa FloatVector)
        @warn "$levels[2] is converted to Float64"
        x2 = convert(Vector{Float64}, x2)
    end

    # check length
    n1 = length(x1)
    n2 = length(x2)
    n1 > 1 || error("size of $levels[1] is not enough: actual size is $n1")
    n2 > 1 || error("size of $levels[2] is not enough: actual size is $n2")

    # paired or independent
    if paired
        n1 == n2 || error("size of $levels[1] and $levels[2] must be the same: ($n1, $n2)")
        if isnothing(exact)
            exact = n1 < 50
        end
        wilcox_test_signed(x1, x2, exact, μ0)
    else
        if isnothing(exact)
            exact = (n1 < 50) & (n2 < 50)
        end
        wilcox_test_independent(x1, x2, levels, exact, μ0)
    end
end

function wilcox_test_signed(
    x1::FloatVector,
    x2::FloatVector,
    exact::Bool,
    μ0::AbstractFloat
)::WilcoxTestResult
    n = length(x1)
    x = x1 .- x2 .- μ0

    x = filter(a -> a ≠ 0.0, x)
    zeroes = length(x) < n
    n = length(x)  # redefinition

    r = tiedrank(abs.(x))
    w_stat = sum(r[x .> 0])

    n_ties = values(countmap(r))
    ties = false
    if length(n_ties) < n
        ties = true
    end

    pval = 1.0
    is_corrected = false
    if exact & !ties & !zeroes
        p = 1.0
        d = WilcoxSignedDist(n)
        if w_stat > (n*(n+1)/4.0)
            p = cdf(d, w_stat-1.0)
        else
            p = ccdf(d, w_stat)
        end
        pval = min(2.0 * p, 1.0)
    else
        is_corrected = true
        z = w_stat - n*(n+1)/4.0
        σ = sqrt(n*(n+1)*(2*n+1)/24.0 -
            sum(n_ties.^3 .- n_ties)/48.0)
        correction = 0.5 * sign(z)
        z = (z - correction) / σ
        norm_std = Normal()
        pval = 2.0 * min(
            D.cdf(norm_std, z),
            D.ccdf(norm_std, z)
        )
    end

    res = WilcoxTestResult(
        :paired, is_corrected, nothing, [n], [mean(x)],
        [std(x)], w_stat, pval, nothing
    )
    return res
end

function wilcox_test_independent(
    x1::FloatVector,
    x2::FloatVector,
    levels::Vector{String},
    exact::Bool,
    μ0::AbstractFloat
)::WilcoxTestResult
    n1 = length(x1)
    n2 = length(x2)
    r = tiedrank(vcat(x1 .- μ0, x2))
    n_ties = values(countmap(r))

    w_stat = sum(r[1:n1]) - n1 * (n1 + 1) / 2
    ties = false
    if length(n_ties) < (n1+n2)
        ties = true
    end

    pval = 1.0
    is_corrected = false
    if exact & !ties
        p = 1.0
        d = WilcoxDist(n1, n2)
        if w_stat > (n1 * n2 / 2)
            p = cdf(d, w_stat-1.0)
        else
            p = ccdf(d, w_stat)
        end
        pval = min(2.0*p, 1.0)
    else
        is_corrected = true
        z = w_stat - n1 * n2 / 2.0
        σ = sqrt((n1*n2/12.0) *
            ((n1+n2+1.0) - sum(n_ties.^3 .- n_ties) / ((n1+n2) * (n1+n2-1.0))))
        correction = 0.5 * sign(z)
        z = (z - correction) / σ
        norm_std = Normal()
        pval = 2.0 * min(
            D.cdf(norm_std, z),
            D.ccdf(norm_std, z)
        )
    end

    m1 = mean(x1)
    m2 = mean(x2)
    v1 = var(x1)
    v2 = var(x2)
    cohen_d = abs(m1-m2) / sqrt((v1 + v2) / 2.0)

    res = WilcoxTestResult(
        :independent, is_corrected, levels, [n1, n2], [m1, m2],
        [sqrt(v1), sqrt(v2)], w_stat, pval, cohen_d
    )
    return res
end

"""
    WilcoxTestModel

A `WilcoxestModel` type for Wilcoxon (signed) rank sum test

# Members

- `X`: model matrix
- `y`: vector to compare its means
- `type`: either `:paired` or `:independent`
- `levels`: `String` vector of the unique levels (only 1 or 2 elements are allowed)
"""
struct WilcoxTestModel{T<:Real} <: TestModel
    x1::AbstractVector{T}
    x2::AbstractVector{T}
    type::Symbol
    levels::Vector{String}
    function WilcoxTestModel(
        X::Matrix{T},
        y::AbstractVector{T},
        type::Symbol,
        levels::Vector{String};
        μ0::T = 0.0
    ) where T
        dummy = X[:,2]
        x1 = y[dummy .== 0]
        x2 = y[dummy .== 1]
        new{T}(x1, x2, type, levels)
    end
end

function fit(obj::WilcoxTestModel)
    if obj.type == :paired
        wilcox_test(obj.x1, obj.x2, paired=true, levels=obj.levels)
    else
        wilcox_test(obj.x1, obj.x2, paired=false, levels=obj.levels)
    end
end

function fit(
        ::Type{WilcoxTestModel},
        X::AbstractMatrix{<:Real},
        y::AbstractVector{<:Real},
        type::Symbol,
        levels::Vector{String};
        μ0::Real = 0.0
    )
    wilcox_test_obj = WilcoxTestModel(X, y, type, levels, μ0=μ0)
    fit(wilcox_test_obj)
end
function fit(
        ::Type{WilcoxTestModel},
        X::AbstractMatrix{<:Real},
        y::AbstractVector{<:Real},
        type::Symbol
    )
    @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
    fit(WilcoxTestModel, X, y, type, ["Sample1", "Sample2"], μ0=0.0)
end

"""
    wilcox_test(X, y, paired::Bool=false)

An alias for `fit(WilcoxTestModel, X, y, levels)`

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and `DataFrame`.

The argument `type` must be either `:paired` or `:independent`.

The argument `levels` can be a `Vector{String}`.
"""
function wilcox_test(
    X::AbstractMatrix{<:Real},
    y::AbstractVector{<:Real},
    paired::Bool=false
)
    if paired
        fit(WilcoxTestModel, X, y, :paired)
    else
        fit(WilcoxTestModel, X, y, :independent)
    end
end
function wilcox_test(
    formula::FormulaTerm{Term,Term},
    df::DataFrame,
    paired::Bool=false
)
    # obtain levels
    f = apply_schema(formula, schema(formula, df))
    levels = DataAPI.levels(df[!,f.rhs.terms[1].sym])
    if eltype(levels) <: Real
        levels = string.(levels)
    end

    # model matrix
    model_frame = ModelFrame(formula, df)
    X = ModelMatrix(model_frame).m
    y = response(model_frame)

    # perform t-test
    if paired
        fit(WilcoxTestModel, X, y, :paired, levels)
    else
        fit(WilcoxTestModel, X, y, :independent, levels)
    end
end


text:
# ---
# title: 1334. Find the City With the Smallest Number of Neighbors at a Threshold Distance
# id: problem1334
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Graph
# link: <https://leetcode.com/problems/find-the-city-with-the-smallest-number-of-neighbors-at-a-threshold-distance/description/>
# hidden: true
# ---
#
# There are `n` cities numbered from `0` to `n-1`. Given the array `edges` where
# `edges[i] = [fromi, toi, weighti]` represents a bidirectional and weighted
# edge between cities `fromi` and `toi`, and given the integer
# `distanceThreshold`.
#
# Return the city with the smallest number **** of cities that are reachable
# through some path and whose distance is **at most** `distanceThreshold`, If
# there are multiple such cities, return the city with the greatest number.
#
# Notice that the distance of a path connecting cities _**i**_ and _**j**_ is
# equal to the sum of the edges' weights along that path.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/16/find_the_city_01.png)
#
#
#
#     Input: n = 4, edges = [[0,1,3],[1,2,1],[1,3,4],[2,3,1]], distanceThreshold = 4
#     Output: 3
#     Explanation: The figure above describes the graph.
#     The neighboring cities at a distanceThreshold = 4 for each city are:
#     City 0 -> [City 1, City 2]
#     City 1 -> [City 0, City 2, City 3]
#     City 2 -> [City 0, City 1, City 3]
#     City 3 -> [City 1, City 2]
#     Cities 0 and 3 have 2 neighboring cities at a distanceThreshold = 4, but we have to return city 3 since it has the greatest number.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/16/find_the_city_02.png)**
#
#
#
#     Input: n = 5, edges = [[0,1,2],[0,4,8],[1,2,3],[1,4,2],[2,3,1],[3,4,1]], distanceThreshold = 2
#     Output: 0
#     Explanation: The figure above describes the graph.
#     The neighboring cities at a distanceThreshold = 2 for each city are:
#     City 0 -> [City 1]
#     City 1 -> [City 0, City 4]
#     City 2 -> [City 3, City 4]
#     City 3 -> [City 2, City 4]
#     City 4 -> [City 1, City 2, City 3]
#     The city 0 has 1 neighboring city at a distanceThreshold = 2.
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 100`
#   * `1 <= edges.length <= n * (n - 1) / 2`
#   * `edges[i].length == 3`
#   * `0 <= fromi < toi < n`
#   * `1 <= weighti, distanceThreshold <= 10^4`
#   * All pairs `(fromi, toi)` are distinct.
#
#
## @lc code=start
using LeetCode

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


text:
# functions to try to debug cui3.out
# converted to test function

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
    const DEBUG = true
end

if isdefined(:PhyloNetworks)
    PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")
else
    CHECKNET || error("need CHECKNET==true to test snaq in test_correctLik.jl")
end

text = "(Xmayae,((Xhellerii,(((Xclemenciae_F2,Xmonticolus):1.458,(((((Xmontezumae,(Xnezahuacoyotl)#H26:0.247::0.804):0.375,((Xbirchmanni_GARC,Xmalinche_CHIC2):0.997,Xcortezi):0.455):0.63,(#H26:0.0::0.196,((Xcontinens,Xpygmaeus):1.932,(Xnigrensis,Xmultilineatus):1.401):0.042):2.439):2.0)#H7:0.787::0.835,(Xmaculatus,(Xandersi,(Xmilleri,((Xxiphidium,#H7:9.563::0.165):1.409,(Xevelynae,(Xvariatus,(Xcouchianus,(Xgordoni,Xmeyeri):0.263):3.532):0.642):0.411):0.295):0.468):0.654):1.022):0.788):1.917)#H27:0.149::0.572):0.668,Xalvarezi):0.257,(Xsignum,#H27:1.381::0.428):4.669);"

#net = readTopology("cui3.out")
net = readTopology(text)
#printPartitions(net)
#printEdges(net)
cleanBL!(net)
cleanAfterRead!(net,false)
length(net.hybrid)
[n.number for n in net.hybrid] == [7,15,25] || error("wrong hybrid nodes")

# ----------------- update everything but partition ------------------
n=net.hybrid[1];
flag, nocycle, edgesInCycle, nodesInCycle = updateInCycle!(net,n);
flag || error("wrong update in hybrid 1")
[n.number for n in nodesInCycle] == [7,-16,-10,-11,-12] || error("wrong update in hybrid 1")
[n.number for n in edgesInCycle] == [16,24,15,9,8] || error("wrong update in hybrid 1")
flag2, edgesGammaz = updateGammaz!(net,n,false);
flag2 || error("wrong update in hybrid 1")
isempty([n.number for n in edgesGammaz]) || error("wrong update in hybrid 1")
!n.isBadDiamondI || error("wrong update in hybrid 1")
!n.isBadDiamondII || error("wrong update in hybrid 1")
flag3, edgesRoot = updateContainRoot!(net,n);
flag3 || error("wrong update in hybrid 1")
[n.number for n in edgesRoot] == [7] || error("wrong update in hybrid 1")

n=net.hybrid[2];
flag, nocycle, edgesInCycle, nodesInCycle = updateInCycle!(net,n);
flag || error("wrong update in hybrid 2")
[n.number for n in nodesInCycle] == [15,-24,-23,-22,-21,-20,-8] || error("wrong update in hybrid 2")
[n.number for n in edgesInCycle] == [31,32,42,43,44,45,26] || error("wrong update in hybrid 2")
flag2, edgesGammaz = updateGammaz!(net,n,false);
flag2 || error("wrong update in hybrid 2")
isempty([n.number for n in edgesGammaz]) || error("wrong update in hybrid 2")
!n.isBadDiamondI || error("wrong update in hybrid 2")
!n.isBadDiamondII || error("wrong update in hybrid 2")
flag3, edgesRoot = updateContainRoot!(net,n);
flag3 || error("wrong update in hybrid 2")
[n.number for n in edgesRoot] == [25,15,9,6,14,12,10,11,13,24,23,19,17,18,22,20,21] || error("wrong update in hybrid 2")

n=net.hybrid[3];
flag, nocycle, edgesInCycle, nodesInCycle = updateInCycle!(net,n);
flag || error("wrong update in hybrid 3")
[n.number for n in nodesInCycle] == [25,-29,-2,-3,-4] || error("wrong update in hybrid 3")
[n.number for n in edgesInCycle] == [53,54,51,49,48] || error("wrong update in hybrid 3")
flag2, edgesGammaz = updateGammaz!(net,n,false);
flag2 || error("wrong update in hybrid 3")
isempty([n.number for n in edgesGammaz]) || error("wrong update in hybrid 3")
!n.isBadDiamondI || error("wrong update in hybrid 3")
!n.isBadDiamondII || error("wrong update in hybrid 3")
flag3, edgesRoot = updateContainRoot!(net,n);
flag3 || error("wrong update in hybrid 3")
[n.number for n in edgesRoot] == [47,5,3,4,46,45,27,44,28,43,29,42,32,30,41,33,40,34,39,35,38,36,37] || error("wrong update in hybrid 3")

isempty(net.partition) || error("wrong partition after initial update")

# -------------- update partition ---------------------
n=net.hybrid[1];
nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,n);
[n.number for n in edgesInCycle] == [16,24,15,9,8] || error("wrong partition hybrid 1")
[n.number for n in nodesInCycle] == [7,-16,-10,-11,-12] || error("wrong partition hybrid 1")
updatePartition!(net,nodesInCycle)
length(net.partition) == 5 || error("wrong partition hybrid 1")
[n.number for n in net.partition[1].edges] == [7] || error("wrong partition hybrid 1")
[n.number for n in net.partition[2].edges] == [23,19,17,18,22,20,21] || error("wrong partition hybrid 1")
[n.number for n in net.partition[3].edges] == [25] || error("wrong partition hybrid 1")
[n.number for n in net.partition[4].edges] == [14,12,10,11,13] || error("wrong partition hybrid 1")
[n.number for n in net.partition[5].edges] == [6] || error("wrong partition hybrid 1")

n=net.hybrid[2];
nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,n);
[n.number for n in nodesInCycle] == [15,-24,-23,-22,-21,-20,-8] || error("wrong partition in hybrid 2")
[n.number for n in edgesInCycle] == [31,32,42,43,44,45,26] || error("wrong partition in hybrid 2")
updatePartition!(net,nodesInCycle)
length(net.partition) == 11 || error("wrong partition hybrid 2")
[n.number for n in net.partition[6].edges] == [30] || error("wrong partition hybrid 2")
[n.number for n in net.partition[7].edges] == [41,33,40,34,39,35,38,36,37] || error("wrong partition hybrid 2")
[n.number for n in net.partition[8].edges] == [29] || error("wrong partition hybrid 2")
[n.number for n in net.partition[9].edges] == [28] || error("wrong partition hybrid 2")
[n.number for n in net.partition[10].edges] == [27] || error("wrong partition hybrid 2")
[n.number for n in net.partition[11].edges] == [46,5,3,4,47] || error("wrong partition hybrid 2")

n=net.hybrid[3];
nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,n);
[n.number for n in nodesInCycle] == [25,-29,-2,-3,-4] || error("wrong partition in hybrid 3")
[n.number for n in edgesInCycle] == [53,54,51,49,48] || error("wrong partition in hybrid 3")
updatePartition!(net,nodesInCycle)
length(net.partition) == 15 || error("wrong partition hybrid 2")
[n.number for n in net.partition[12].edges] == [52] || error("wrong partition hybrid 3")
[n.number for n in net.partition[13].edges] == [1] || error("wrong partition hybrid 3")
[n.number for n in net.partition[14].edges] == [50] || error("wrong partition hybrid 3")
[n.number for n in net.partition[15].edges] == [2] || error("wrong partition hybrid 3")

#printPartitions(net)


text:
using NonlinearBandits
using Test
using CUDA
using Distributions: Uniform, Normal
using Flux: gpu, cpu
import Base: ==

function mae(x1, x2)
    return sum(abs.(x1 - x2)) / length(x1)
end

@testset "BayesLM.jl" begin
    include("test_BayesLM.jl")
end

@testset "BayesPM.jl" begin
    include("test_BayesPM.jl")
end

@testset "PartitionedBayesPM.jl" begin
    include("test_PartitionedBayesPM.jl")
end

@testset "bandits.jl" begin
    include("test_bandits.jl")
end

@testset "contexts.jl" begin
    include("test_contexts.jl")
end

@testset "rewards.jl" begin
    include("test_rewards.jl")
end

@testset "metrics.jl" begin
    include("test_metrics.jl")
end

@testset "drivers.jl" begin
    include("test_drivers.jl")
end

@testset "PolynomialThompsonSampling." begin
    include("test_PolynomialThompsonSampling.jl")
end

@testset "NeuralLinear.jl" begin
    include("test_NeuralLinear.jl")
end

@testset "NearestNeighboursUCB.jl" begin
    include("test_NearestNeighboursUCB.jl")
end


text:
export pprint, pwrite

################################################################
# I/O
import EzXML: parsexml, parsehtml, readxml, readhtml
# from EzXML
funs = [:parsexml, :parsehtml, :readxml, :readhtml]
for fun in funs

    @eval begin

        # doc
        @doc (@doc $(fun)) $(fun)

        # exporting
        export $(fun)
    end
end
################################################################
"""
    pprint(x)
    pprint(io, x)
    pprint(filename, x)

Pretty prints the xml/html content of a aml type. Also, pretty prints a Node or Document type.
"""
function pprint(io::IO, x::T) where {T<:Union{Document, Node}}
    println("")
    prettyprint(io, x)
end

function pprint(x::T) where {T<:Union{Document, Node}}
    println("")
    prettyprint(x)
end

function pprint(filename::AbstractString, x::T) where {T<:Union{Document, Node}}
    open(io->pprint(io, x), filename, "w")
end

function pprint(io::IO, x)
    println("")
    prettyprint(io, x.aml)
end

function pprint(x)
    println("")
    prettyprint(x.aml)
end

function pprint(filename::AbstractString, x)
    open(io->pprint(io, x.aml), filename, "w")
end

const pwrite = pprint
const write = pprint
const print = pprint


text:
module NodesArrayTests

using Gridap, Test

p = Polytope(HEX_AXIS)
test_hex_polytopes(p)

p = Polytope(HEX_AXIS, HEX_AXIS)
test_hex_polytopes(p)

p = Polytope(HEX_AXIS, HEX_AXIS, HEX_AXIS)
test_hex_polytopes(p)

p = Polytope(HEX_AXIS,TET_AXIS)
# _high_order_test(p,Float64,1)
T = Float64
horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, 2)
@test length(horeffe.shfbasis) == 6
horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, 3)
@test length(horeffe.shfbasis) == 10

p = Polytope(HEX_AXIS,TET_AXIS,TET_AXIS)
# _high_order_test(p,Float64,1)
T = Float64
horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, 2)
@test length(horeffe.shfbasis) == 10
horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, 3)
@test length(horeffe.shfbasis) == 20

# Problems with grid for new machinery
p = Polytope(HEX_AXIS,TET_AXIS,TET_AXIS)
LagrangianRefFE(Float64,p,3)
grid = Grid(p,2)
# writevtk(grid,"grid")

function _high_order_test_scalar(p,T,order)
  horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, order)
  if (Gridap.RefFEs._is_hex(p)) @test length(horeffe.shfbasis) == (order+1)^(length(p.extrusion)) end
  # if (Gridap.RefFEs._is_tet(p)) @test length(horeffe.shfbasis) == sum(1:order+1) end
end

function _high_order_test_vector(p,T,order)
  horeffe = Gridap.RefFEs._high_order_lagrangian_reffe(p, T, order)
  if (Gridap.RefFEs._is_hex(p)) @test length(horeffe.shfbasis) == length(p.extrusion)*(order+1)^(length(p.extrusion)) end
end

function test_hex_polytopes(p)
  _high_order_test_scalar(p,Float64,1)
  _high_order_test_scalar(p,Float64,2)
  _high_order_test_scalar(p,Float64,3)
  _high_order_test_scalar(p,Float64,4)
  T = VectorValue{length(p.extrusion),Float64}
  _high_order_test_vector(p,T,1)
  _high_order_test_vector(p,T,2)
  _high_order_test_vector(p,T,3)
  _high_order_test_vector(p,T,4)
end

end #module


text:
struct CacheDict{T <: Array, N}
    S::Int
    n::Int
    caches::Dict{UInt64, Array}

    function CacheDict{T,N}(S,n) where {T,N}
        new{T,N}(S, n, Dict{UInt64, Array}())
    end
end

const ArrayCache{N} = CacheDict{Array{T,N} where {T <: Number, N}}
const VectorCache = ArrayCache{1}
const MatrixCache = ArrayCache{2}

const VectorArrayCache{N} = CacheDict{Vector{Array{T,N} where T <: Number}, N}
const VectorMatrixCache = VectorArrayCache{2}
const VectorVectorCache = VectorArrayCache{1}

function _new_cache_entry(c::ArrayCache{N}, ST::DataType) where {N}
    ones(ST, Tuple([c.n for _ in 1:N])...)
end

function _new_cache_entry(c::VectorArrayCache{N}, ST::DataType) where {N}
    [ones(ST, Tuple([c.n for _ in 1:N])...) for _ in 1:c.S]
end

_get_type(::ArrayCache{N}, ST) where {N} = Array{ST,N}
_get_type(::VectorArrayCache{N}, ST) where {N} = Vector{Array{ST,N}}

@inline function Base.getindex(c::CacheDict, name::Symbol, ST::DataType)
    key = hash(Threads.threadid(), hash(name, hash(ST)))
    if haskey(c.caches, key)
        c.caches[key]
    else
        c.caches[key] = _new_cache_entry(c, ST)
    end::_get_type(c, ST)
end

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
using SafeTestsets

@safetestset "DiscreteEventsLite" begin
    using DiscreteEventsLite, DataStructures, Test
    f(a; k) = (a, k)
    scheduler = Scheduler()
    register!(scheduler, f, after, .99, "hi "; k=1, description="123")
    event = peek(scheduler.events).first
    @test event.description == "123"
    run!(scheduler, 11)
    @test scheduler.time ≈ 11
    @test length(scheduler.events) == 0

    scheduler = Scheduler()
    register!(scheduler, f, after, .99, "hi "; k=1)
    register!(scheduler, stop!, .5, scheduler)
    run!(scheduler, 11)
    @test scheduler.time ≈ .5
    @test length(scheduler.events) == 2

    scheduler = Scheduler()
    register!(scheduler, f, after, .99, "hi "; k=1)
    event = dequeue!(scheduler.events)
    v1,v2 = event.fun()
    @test v1 == "hi "
    @test v2 == 1

    scheduler = Scheduler()
    register!(scheduler, f, after, .99, "hi "; k=1, id="1")
    register!(scheduler, f, after, 1.99, "hi "; k=1, id="2")
    remove_events!(scheduler, "1")
    event = dequeue!(scheduler.events)
    @test event.id == "2"

    scheduler = Scheduler()
    register!(scheduler, f, every, 1.0, "hi "; k=1)
    run!(scheduler, 10)
    @test scheduler.time == 10

    scheduler = Scheduler()
    register!(scheduler, f, every, 1.0, "hi "; k=1)
    register!(scheduler, ()->(), at, 10.5)
    run!(scheduler, 10)
    @test scheduler.time == 10

    scheduler = Scheduler()
    register!(scheduler, f, every, 1.0, "hi "; k=1)
    register!(scheduler, ()->(), at, 10.5)
    register!(scheduler, stop!, 5.0, scheduler)
    run!(scheduler, 10)
    @test scheduler.time == 5

    scheduler = Scheduler()
    register!(scheduler, f, every, 1.0, "hi "; k=1)
    register!(scheduler, ()->(), at, 10.5)
    register!(scheduler, stop!, 5.0, scheduler)
    DiscreteEventsLite.reset!(scheduler)
    @test scheduler.time == 0
    @test isempty(scheduler.events)

    scheduler = Scheduler(;store=true)
    register!(scheduler, f, every, 1.0, "hi "; k=1)
    register!(scheduler, ()->(), at, 10.5)
    run!(scheduler, 11)
    @test length(scheduler.complete_events) == 13
end


text:
import Base: diff

immutable Equation
    left :: FieldExpr
    right :: FieldExpr
end

for op in (:+, :-, :*, :/, :^, :equal)
    if op == :equal
        @eval ($op)(x::FieldExpr, y::FieldExpr) = Equation(x, y)
    else
        @eval begin
            @eval import Base: $op
            ($op)(x::FieldExpr, y::FieldExpr) = FieldExprNode($(Meta.quot(op)), [x, y])
        end
    end

    # ^ operator is a special case since there is a shortcut for ^(x, y::Integer) in Base,
    # that interferes with our generic definitions.
    # So we have to add specific definitions for Integer to avoid ambiguity.

    for ext_type in (Number, (op == :^ ? (Integer,) : ())...)
        @eval begin
            ($op)(x::FieldExpr, y::$ext_type) = ($op)(x, convert(FieldExpr, y))
            ($op)(x::$ext_type, y::FieldExpr) = ($op)(convert(FieldExpr, x), y)
        end
    end
end

for builtin_function in (:abs, :cosh)
    @eval begin
        import Base: $builtin_function

        function $builtin_function(x::FieldExpr, args...)
            # This check is not necessary, [x, args...] should work anyway,
            # but it currently doesn't because of bug #17003
            if length(args) > 0
                args = [convert(FieldExpr, arg) for arg in args]
                args = [x, args...]
            else
                args = [x]
            end
            FieldExprNode($(Meta.quot(builtin_function)), args)
        end
    end
end

# Taken from StatsBase.jl
function countmap{T}(x::AbstractArray{T,1})
    cm = Dict{T,Int}()
    for v in x
        cm[v] = get(cm, v, 0) + 1
    end
    return cm
end

function diff(f::UnknownField, dimensions...)
    # here we can check that the field actually have the required dimensions
    grouped_dims = countmap(collect(dimensions))
    Diff(f, grouped_dims)
end

function fourier(fexpr::FieldExpr, dimension::Dimension; inverse=false)
    fourier(fexpr, Set{Dimension}([dimension]), inverse=inverse)
end

function fourier(fexpr::FieldExpr, dimension_pairs...; inverse=false)
    if !inverse
        for (dim, fdim) in dimension_pairs
            @assert fdim.base_dimension === dim
        end
        dimension_dict = Dict{Dimension, Dimension}(dimension_pairs)
        Fourier(fexpr, false, dimension_dict)
    else
        for (fdim, dim) in dimension_pairs
            @assert fdim.base_dimension === dim
        end
        dimension_dict = Dict{Dimension, Dimension}(dimension_pairs)
        Fourier(fexpr, true, dimension_dict)
    end
end

function stochastic_mean(fexpr::FieldExpr)
# If the field is not stochastic, should we throw an error?
    if is_stochastic(fexpr)
        StochasticMean(fexpr)
    else
        error("Requires a stochastic field")
    end
end

function stochastic_std(fexpr::FieldExpr)
# If the field is not stochastic, should we throw an error?
    if is_stochastic(fexpr)
        sqrt(stochastic_mean(fexpr^2) - stochastic_mean(fexpr)^2)
    else
        error("Requires a stochastic field")
    end
end

function sample(fexpr::Number, trajectories)
    return Sample(convert(FieldExpr, fexpr), trajectories)
end

function sample(fexpr::FieldExpr, trajectories)
    return Sample(fexpr, trajectories)
end

function random_normal(dimensions, mean, std)
    # TODO: here we need to add the dV factor too, depending on what dimensions are requested
    RandomVariable(:random_normal, dimensions) * std + mean
end


text:
#=
LatinStr/_LatinStr type (ISO Latin1 8-bit subset of Unicode)

Copyright 2017 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on code for ASCIIString that used to be in Julia
=#

## overload methods for efficiency ##

is_ascii(str::MaybeSub{<:Str{<:_LatinCSE}}) = false
is_latin(str::MaybeSub{<:Str{<:LatinCSE}}) = true
is_bmp(str::MS_Latin) = true
is_unicode(str::MS_Latin) = true

const MS_ASCIILatin = MaybeSub{<:Str{<:Union{ASCIICSE, Latin_CSEs}}}

function string(collection::MS_ASCIILatin...)
    length(collection) == 1 && return collection[1]
    len = 0
    @inbounds for str in collection
        len += ncodeunits(str)
    end
    buf, pnt = _allocate(len)
    @inbounds for str in collection
        len = ncodeunits(str)
        _memcpy(pnt, pointer(str), len)
        pnt += len
    end
    Str(LatinCSE, buf)
end

## transcoding to Latin1 ##

function convert(::Type{<:Str{C}}, str::AbstractString) where {C<:Latin_CSEs}
    is_empty(str) && return C === _LatinCSE ? empty_ascii : empty_latin
    # Might want to have invalids_as here
    len, flags, num4byte, num3byte, num2byte = unsafe_check_string(str)
    num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
    buf, pnt = _allocate(UInt8, len)
    @inbounds for ch in str
        set_codeunit!(pnt, ch%UInt8)
        pnt += 1
    end
    Str((C === _LatinCSE && flags == 0) ? ASCIICSE : C, buf)
end

const ASCIISubStr  = SubString{<:Str{ASCIICSE}}
const LatinSubStr  = SubString{<:Str{LatinCSE}}
const _LatinSubStr = SubString{<:Str{_LatinCSE}}

_cpyconvert(::Type{C}, str) where {C} = Str(C, _copysub(str))::Str{C,Nothing,Nothing,Nothing}
_cpysubset(::Type{_LatinCSE}, str) = Str(is_ascii(s) ? ASCIICSE : _LatinCSE, _copysub(s))

convert(::Type{<:Str{LatinCSE}}, s::Str{LatinCSE})  = s
convert(::Type{<:Str{LatinCSE}}, s::Str{ASCIICSE})  = _cpyconvert(LatinCSE, s)
convert(::Type{<:Str{LatinCSE}}, s::Str{_LatinCSE}) = _cpyconvert(LatinCSE, s)

convert(::Type{<:Str{LatinCSE}}, s::LatinSubStr)  = _cpyconvert(LatinCSE, s)
convert(::Type{<:Str{LatinCSE}}, s::ASCIISubStr)  = _cpyconvert(LatinCSE, s)
convert(::Type{<:Str{LatinCSE}}, s::_LatinSubStr) = _cpyconvert(LatinCSE, s)

convert(::Type{<:Str{_LatinCSE}}, s::Str{LatinCSE})  = s
convert(::Type{<:Str{_LatinCSE}}, s::Str{ASCIICSE})  = s
convert(::Type{<:Str{_LatinCSE}}, s::Str{_LatinCSE}) = s

convert(::Type{<:Str{_LatinCSE}}, s::LatinSubStr)  = _cpysubset(_LatinCSE, s)
convert(::Type{<:Str{_LatinCSE}}, s::ASCIISubStr)  = _cpyconvert(LatinCSE, s)
convert(::Type{<:Str{_LatinCSE}}, s::_LatinSubStr) = Str(_LatinCSE, _copysub(s))

# Assumes that has already been checked for validity

# These should be sped up to do chunks at a time, when no bytes > 0x7f
function _utf8_to_latin(pnt::Ptr{UInt8}, len)
    buf, out = _allocate(UInt8, len)
    fin = out + len
    while out < fin
        ch = get_codeunit(pnt)
        # Handle ASCII characters
        ch > 0x7f && (ch = ((ch & 3) << 6) | (get_codeunit(pnt += 1) & 0x3f))
        pnt += 1
        set_codeunit!(out, ch)
        out += 1
    end
    buf
end

function _latin_to_utf8(pnt::Ptr{UInt8}, len)
    buf, out = _allocate(UInt8, len)
    fin = out + len
    while out < fin
        ch = get_codeunit(pnt)
        pnt += 1
        if ch > 0x7f
            set_codeunit!(out, 0xc0 | (ch >>> 6))
            ch = 0x80 | (ch & 0x3f)
            out += 1
        end
        set_codeunit!(out, ch)
        out += 1
    end
    buf
end

# Fast conversion from LatinStr or _LatinStr to RawUTF8Str/UTF8Str
function convert(::Type{<:Str{C}}, str::MS_Latin) where {C<:UTF8_CSEs}
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_str(C)
    @preserve str begin
        pnt = pointer(str)
        Str(C,(cnt = count_latin(len, pnt)) == 0 ? _copysub(str) : _latin_to_utf8(pnt, len + cnt))
    end
end

# Fast conversion from LatinStr or _LatinStr to String
function convert(::Type{String}, str::MS_Latin)
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_string
    @preserve str begin
        pnt = pointer(str)
        (cnt = count_latin(len, pnt)) == 0 ? _copysub(str) : _latin_to_utf8(pnt, len + cnt)
    end
end

function convert(::Type{<:Str{C}}, str::MaybeSub{<:Str{<:Union{Word_CSEs,Quad_CSEs}}}
                 ) where {C<:Latin_CSEs}
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return C === _LatinCSE ? empty_ascii : empty_latin
    @preserve str begin
        pnt = pointer(str)
        # get number of bytes to allocate
        len, flags, num4byte, num3byte, num2byte, latinbyte = fast_check_string(pnt, len)
        num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
        Str((C === _LatinCSE && latinbyte == 0) ? ASCIICSE : C, _cvtsize(UInt8, pnt, len))
    end
end

function convert(::Type{<:Str{LatinCSE}}, str::MS_ByteStr)
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_str(LatinCSE)
    @preserve str begin
        pnt = pointer(str)
        # get number of bytes to allocate
        len, flags, num4byte, num3byte, num2byte, latinbyte = fast_check_string(pnt, siz)
        num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
        Str(LatinCSE, flags == 0 ? _copysub(str) : _utf8_to_latin(pnt, len))
    end
end

convert(::Type{<:Str{_LatinCSE}}, ch::Unsigned) =
    (ch <= 0xff
     ? _convert(ch <= 0x7f ? ASCIICSE : _LatinCSE, ch%UInt8)
     : strerror(StrErrors.LATIN1, ch))

function convert(::Type{<:Str{_LatinCSE}}, str::MS_ByteStr)
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_ascii
    @preserve str begin
        pnt = pointer(str)
        # get number of bytes to allocate
        len, flags, num4byte, num3byte, num2byte, latinbyte = fast_check_string(pnt, siz)
        num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
        Str(latinbyte == 0 ? ASCIICSE : _LatinCSE,
            flags == 0 ? _copysub(str) : _utf8_to_latin(pnt, len))
    end
end

function convert(::Type{<:Str{C}}, vec::Vector{CU}) where {C<:Latin_CSEs,CU<:CodeUnitTypes}
    # handle zero length string quickly
    (len = length(vec)) == 0 && return _empty_str(C)
    @preserve vec begin
        pnt = pointer(vec)
        # get number of bytes to allocate
        len, flags, num4byte, num3byte, num2byte, latinbyte = fast_check_string(pnt, len)
        num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
        Str((C === _LatinCSE && latinbyte == 0) ? ASCIICSE : C,
            ((CU !== UInt8 || flags == 0) ? _cvtsize(UInt8, pnt, len) : _utf8_to_latin(pnt, len)))
    end
end

function convert(::Type{<:Str{C}}, str::MaybeSub{<:Str{<:Union{Text2CSE,Text4CSE}}}
                 ) where {C<:Latin_CSEs}
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_str(C)
    @preserve str begin
        pnt = pointer(str)
        # get number of bytes to allocate
        len, flags, num4byte, num3byte, num2byte, latinbyte = fast_check_string(pnt, len)
        num4byte + num3byte + num2byte == 0 || strerror(StrErrors.INVALID_LATIN1)
        Str((C === _LatinCSE && latinbyte == 0) ? ASCIICSE : C, _cvtsize(UInt8, pnt, len))
    end
end


text:
using XFields
using LinearAlgebra
using FFTransforms
using Test

# mean(x) = sum(x) / length(x)

# @testset "basic" begin
#     include("basic.jl")
# end

# @testset "rFFTimpulses" begin
#     include("rFFTimpulses/test_get_rFFTimpulses_1.jl")
#     include("rFFTimpulses/test_get_rFFTimpulses_2.jl")
# end

@testset "Xfield constructors" begin

	let
		Tf = Float32
		n1, n2 = 100, 256
		p2 = 2π
		ft = 𝕀(n1) ⊗ 𝕎(Tf,n2,p2)

		fx = randn(Tf, size_in(ft))
		fk = plan(ft) * fx
		gk = plan(ft) * randn(Tf, size_in(ft))
		gx = plan(ft) \ gk

		fmap     = @inferred Xmap(ft, fx)
		ffourier = @inferred Xfourier(ft, gk)

		@inferred ft * fmap
		@inferred ft \ fmap
		@inferred ft * ffourier
		@inferred ft \ ffourier

		@test sum(abs2, fmap[:] .- fx) ≈ 0
		@test sum(abs2, fmap[!] .- fk) ≈ 0
		@test sum(abs2, ffourier[!] .- gk) ≈ 0
		@test sum(abs2, ffourier[:] .- gx) ≈ 0
	end

	let
		Tf = Complex{Float64}
		n1, n2 = 111, 256
		ft = 𝕎(Tf,n1,2π) ⊗ 𝕀(n2)

		fx = rand(Bool, size_in(ft))
		fk = plan(ft) * Tf.(fx)
		gk = rand(Bool, size_out(ft))
		Ti = eltype_out(ft)
		gx = plan(ft) \ (Ti.(gk))

		fmap     = @inferred Xmap(ft, fx)
		gfourier = @inferred Xfourier(ft, gk)

		@test sum(abs2, fmap[:] .- fx) ≈ 0
		@test sum(abs2, fmap[!] .- fk) ≈ 0
		@test sum(abs2, gfourier[!] .- gk) ≈ 0
		@test sum(abs2, gfourier[:] .- gx) ≈ 0
	end

	let
		Tf = Complex{Float64}
		n1, n2 = 100, 256
		p1, p2 = 1.0, 2π
		ft = 𝕎(Tf,n1,p1) ⊗ 𝕎(Tf,n2,p2)

		ax = randn(Tf, size_in(ft))
		ak = plan(ft) * ax
		amap = @inferred Xmap(ft, ax)
		bmap = @inferred Xmap(ft,1)
		cmap = @inferred Xmap(ft)

		dk = plan(ft) * randn(Tf, size_in(ft))
		dx = plan(ft) \ dk
		dfourier = @inferred Xfourier(ft, dk)
		efourier = @inferred Xfourier(ft, 1)
		ffourier = @inferred Xfourier(ft)

		@test sum(abs2, amap[:] .- ax) ≈ 0
		@test sum(abs2, amap[!] .- ak) ≈ 0
		@test sum(abs2, bmap[:] .- 1) ≈ 0
		@test sum(abs2, cmap[:] .- 0) ≈ 0

		@test sum(abs2, dfourier[!] .- dk) ≈ 0
		@test sum(abs2, dfourier[:] .- dx) ≈ 0
		@test sum(abs2, efourier[!] .- 1) ≈ 0
		@test sum(abs2, ffourier[!] .- 0) ≈ 0
	end

end

@testset "Xfield DiagOp" begin

	let
		pd  = (1.0, 3.5)
		sz  = (64,28)
		Tf  = Float64
		Ti  = Complex{Tf}
		ft  = 𝕎(Tf,sz,pd)
		F   = typeof(ft)

		Xf = Xfourier{F, Tf, Ti, 2}
		Xm = Xmap{F, Tf, Ti, 2}

		fm = @inferred Xm(ft,rand(Tf,sz))
		ff = @inferred Xf(Xm(ft,rand(Tf,sz)))

		@inferred Xfourier(fm)
		@inferred Xf(fm)
		@inferred Xmap(fm)
		@inferred Xm(fm)

		@inferred Xfourier(ff)
		@inferred Xf(ff)
		@inferred Xmap(ff)
		@inferred Xm(ff)

		@inferred fm + ff
		@inferred - 2 * ff

		L = DiagOp(fm)
		L.f * Xmap(ff)
		ff isa Xf
		L * ff

		L1 = DiagOp(fm+1)
		L2 = DiagOp(ff+1)

		L1 * L2
		L2 * L2
		L1 * L1

		L2 * L2 * fm
		L2 * (L2 * fm)

		L1 * L1 * ff
		L1 * (L1 * ff)

		L1 * L2 * L2

		L1 * L2 \ L2

		for f ∈ (fm, ff), L ∈ (L1, L2), M ∈ (L1, L2)
			@inferred L * f
			@inferred L * M * f
			@inferred L * M \ f
			@inferred L \ M \ f
			@inferred L \ M * f
			@test sum(abs2, fielddata((L * M * f) - (L * (M * f)))) ≈ 0.0 atol=1e-5
			@test sum(abs2, fielddata((L * M \ f) - (inv(M) * (inv(L) * f)))) ≈ 0.0 atol=1e-5
			@test sum(abs2, fielddata((L \ M \ f) - (inv(M) * (inv(inv(L)) * f)))) ≈ 0.0 atol=1e-5
			@test sum(abs2, fielddata((L \ M * f) - (inv(L) * (M * f)))) ≈ 0.0 atol=1e-5
		end

	end
end

@testset "Uniform scaling" begin

	let nᵢ=(256, 50), pᵢ=(1.0, 0.5), Tf  = Float64

		Ti  = Complex{Tf}
		ft  = 𝕎(Tf,nᵢ,pᵢ)

		fmap    = Xmap(ft, rand(Bool, size_in(ft)))
		gfourier = Xfourier(ft, rand(Bool, size_out(ft)))

		Dmap = DiagOp(fmap)
		Dfourier = DiagOp(gfourier)

		Id2 = 2*LinearAlgebra.I

		@inferred Id2 * fmap
		@inferred Id2 * gfourier
		@inferred Id2 \ fmap
		@inferred Id2 \ gfourier

		@inferred Dmap * Id2
		@inferred Dfourier * Id2
		@inferred Dmap + Id2
		@inferred Dfourier + Id2
		@inferred Dmap - Id2
		@inferred Dfourier - Id2

		@inferred Id2 * Dmap
		@inferred Id2 * Dfourier
		@inferred Id2 + Dmap
		@inferred Id2 + Dfourier
		@inferred Id2 - Dmap
		@inferred Id2 - Dfourier

		# TODO add actual tests here ..
	end

end

#=

#' SField Constructors
#' ------------------------------------------

let FT = rFFT(nᵢ=(256, 50), pᵢ=(1.0, 0.5))

	@inferred Rmap{FT}()
	@inferred Rmap{FT}(0)
	@inferred Rfourier{FT}()
	@inferred Rfourier{FT}(0)

	matx = rand(Bool, Grid(FT).nxi...)
	matk = rand(Bool, Grid(FT).nki...)

	@inferred Rmap{FT}(matx)
	@inferred Rfourier{FT}(matk)

end

#' Test using FT for converting between Rmap and Rfourier
#' ------------------------------------------

let FT = rFFT(nᵢ=(256, 50), pᵢ=(1.0, 0.5))

	matx = rand(Float64, Grid(FT).nxi...)
	matk = FT*rand(Float64, Grid(FT).nxi...)
	fx = matx |> Rmap{FT}
	fk = matk |> Rfourier{FT}

	@inferred FT * fx
	@inferred FT * fk
	@inferred FT \ fx
	@inferred FT * fk

	@test all((FT * fx).k .== (FT * matx))
	@test all((FT \ fx).x .== matx)

	@test all((FT * fk).k .== matk)
	@test all((FT \ fk).x .== (FT \ matk))

end

#' ======================================================
pᵢ  = (1.0, 0.5) # periods
nᵢ   = (256, 256) # number of samples (left endpoint included)
FT  = rFFT{Float32,nᵢ,pᵢ,length(nᵢ)}
UFT = rFFTunitary{Float32,nᵢ,pᵢ,length(nᵢ)}

FT_plan = plan(FT)
UFT_plan = plan(UFT)

grid   = Grid(FT)
λ      = wavenumber(FT)
kfulli = frequencies(FT)
xfulli = pixels(FT)

fk = FT * rand(grid.nxi...)
fx = FT \ rand(Complex{Float32}, grid.nki...)

f1 = Rmap{FT}(fx)
f2 = Rfourier{FT}(fk)
f4 = Rfourier{FT}(f1)
f3 = Rmap{FT}(f2)

@test all(f1[!] .== FT_plan * fx)
@test all(f1[:] .== fx)
@test all(f2[!] .== fk)
@test all(f2[:] .== FT_plan \ fk)

@inferred dot(f3, f1)
# real(dot(cplan(FT) * f3[:x], cplan(FT) * f1[:x])) / Grid(FT).Ωk
# real(dot(f3[:x], f1[:x])) / Grid(FT).Ωx

f1 + f2
- 2 * f2

L = DiagOp(f1)
L * f2

L1 = DiagOp(f1)
L2 = DiagOp(f2)
L3 = sqrt(L2)

@test all(L1[!] .== FT_plan * fx)
@test all(L1[:] .== fx)
@test all(L2[!] .== fk)
@test all(L2[:] .== FT_plan \ fk)

L1 * f2
L2 * L1 * f2
L2 * L1

#' ======================================================
pᵢ  = (1.0,)
nᵢ   = (256,)
FT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}

grid   = Grid(FT)
λ      = wavenumber(FT)
kfulli = frequencies(FT)
xfulli = pixels(FT)

fk = FT * rand(grid.nxi...)
fx = FT \ rand(Complex{Float64}, grid.nki...)

f1 = Rmap{FT}(fx)
f2 = Rfourier{FT}(fk)

f4 = Rfourier{FT}(f1)
f3 = Rmap{FT}(f2)

f1 + f2
- 2 * f2

L = DiagOp(f1)
L * f2

#' =========================================

nᵢ  = (256, 50) # number of samples (left endpoint included)
pᵢ  = (1.0, 0.5) # periods
Δxᵢ = tuple((pn[1]/pn[2] for pn in zip(pᵢ,nᵢ))...)

FFT1 = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
FFT2 = rFFT(nᵢ=nᵢ, pᵢ=pᵢ)
FFT3 = rFFT(nᵢ=nᵢ, Δxᵢ=Δxᵢ)
FFTu1 = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
FFTu2 = rFFTunitary(nᵢ=nᵢ, pᵢ=pᵢ)
FFTu3 = rFFTunitary(nᵢ=nᵢ, Δxᵢ=Δxᵢ)

@test FFT1 == FFT2 == FFT3
@test FFTu1 == FFTu2 == FFTu3

grid1 = Grid(FFT1)
grid2 = Grid(FFT2)
grid3 = Grid(FFT3)
gridu1 = Grid(FFTu1)
gridu2 = Grid(FFTu2)
gridu3 = Grid(FFTu3)

@test grid1.ki == grid2.ki == grid3.ki
@test gridu1.ki == gridu2.ki == gridu3.ki

#' =========================================

pᵢ  = (1.0, 0.5) # periods
nᵢ   = (256, 256) # number of samples (left endpoint included)
cFT  = cFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
cUFT  = cFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
rFT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
rUFT  = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}

grid = Grid(cFT)

@inferred Grid(cFT)
@inferred wavenumber(cFT)
@inferred frequencies(cFT)
@inferred pixels(cFT)

cfx = rand(Complex{Float64}, grid.nxi...)
rfx = rand(Float64, grid.nxi...)

@test grid.nxi == grid.nki

@test sum(abs2, plan(cFT)  * cfx .- cplan(rFT)  * cfx) ≈ 0.0
@test sum(abs2, plan(cUFT) * cfx .- cplan(rUFT) * cfx) ≈ 0.0

@test sum(abs2, plan(cFT)  * cfx .- cplan(cFT)  * cfx) ≈ 0.0
@test sum(abs2, plan(cUFT) * cfx .- cplan(cUFT) * cfx) ≈ 0.0

@test sum(abs2, plan(rFT)  * rfx .- rplan(rFT)  * rfx) ≈ 0.0
@test sum(abs2, plan(rUFT) * rfx .- rplan(rUFT) * rfx) ≈ 0.0

@test sum(abs2, plan(rFT)  * rfx .- rplan(cFT)  * rfx) ≈ 0.0
@test sum(abs2, plan(rUFT) * rfx .- rplan(cUFT) * rfx) ≈ 0.0

@inferred adjoint(plan(cFT))
@inferred adjoint(plan(rFT))
@inferred transpose(plan(cFT))
@inferred transpose(plan(rFT))

=#


text:
using IncrementalInference, TransformUtils
using Base.Test

using KernelDensityEstimate

N = 200
fg = emptyFactorGraph()

initCov = eye(6)
[initCov[i,i] = 0.01 for i in 4:6];
odoCov = deepcopy(initCov)

println("Adding PriorPose3 to graph...")
v1 = addNode!(fg,"x1",  0.1*randn(6,N),  N=N)
initPosePrior = PriorPose3(SE3(0), initCov)
f1  = addFactor!(fg,[v1], initPosePrior)

println("Ensure vertex initialized properly")
# start with to tight an initialization
muX1 = Base.mean(getVal(fg,"x1"),2)
stdX1 = Base.std(getVal(fg,"x1"),2)
@test sum(map(Int,abs(muX1) .< 0.1)) == 6
@test sum(map(Int, 0.05 .< stdX1 .< 0.15)) == 6

println("Testing PriorPose3 evaluation...")
priorpts = evalFactor2(fg, fg.g.vertices[2], 1)
means = Base.mean(priorpts,2)
@test sum(map(Int,abs(means[1:3]) .> 0.5)) == 0
@test sum(map(Int,abs(means[4:6]) .> 0.05)) == 0

println("Adding Pose3Pose3 to graph...")
odo = SE3([10;0;0], Quaternion(0))
pts0X2 = projectParticles(getVal(fg,"x1"), odo, odoCov)
odoconstr = Pose3Pose3(odo, odoCov)
v2 = addNode!(fg,"x2",  pts0X2, N=N)
addFactor!(fg,[v1;v2],odoconstr)

println("Testing Pose3Pose3 evaluation...")
X1pts = evalFactor2(fg, fg.g.vertices[4], 1)
X2pts = evalFactor2(fg, fg.g.vertices[4], 3)
X2ptsMean = Base.mean(X2pts,2)
X1ptsMean = Base.mean(X1pts,2)
@test  sum(map(Int, abs(X1ptsMean) .< 0.5 )) == 6
@test  sum(map(Int, abs(X2ptsMean - [10.0;0;0;0;0;0]) .< 0.5 )) == 6

println("Construct Bayes tree and perform inference...")
tree = prepBatchTree!(fg);
inferOverTree!(fg, tree)

println("Ensure basic parameters on x1,x2 after inference...")
# check mean and covariances after one up and down pass over the tree
muX1 = Base.mean(getVal(fg,"x1"),2)
stdX1 = Base.std(getVal(fg,"x1"),2)
@test sum(map(Int,abs(muX1[1:3]) .< 1.0)) == 3
@test sum(map(Int,abs(muX1[4:6]) .< 0.1)) == 3
@test sum(map(Int, 0.5 .< stdX1[1:3] .< 1.5)) == 3
@test sum(map(Int, 0.05 .< stdX1[4:6] .< 0.15)) == 3
muX2 = Base.mean(getVal(fg,"x2"),2)
stdX2 = Base.std(getVal(fg,"x2"),2)
@test sum(map(Int, abs(muX2[1:3]-[10.0;0;0]) .< 1.0)) == 3
@test sum(map(Int, abs(muX2[4:6]) .< 0.1)) == 3
@test sum(map(Int, 1.0 .< stdX2[1:3] .< 2.0)) == 3
@test sum(map(Int, 0.05 .< stdX2[4:6] .< 0.15)) == 3

# println("Plot marginals to see what is happening")
# plotKDE(marginal(getVertKDE(fg,"x1"),[1]))
# plotKDE(marginal(getVertKDE(fg,"x2"),[1]))

println("Modify factor graph, adding Pose3Pose3's X3,X4 to graph...")
#X3, yaw 90 degrees
odo2 = SE3([0;0;0], AngleAxis(pi/2,[0;0;1]))
odoconstr2 = Pose3Pose3(odo2, odoCov)
pts0x3 = projectParticles(getVal(fg,"x2"), odo2, odoCov) # init new pose location
v3 = addNode!(fg,"x3", pts0x3, N=N)
addFactor!(fg,[v2;v3],odoconstr2)

# X4, drive forward 10m
odo3 = SE3([10;0;0], SO3(0))
odoconstr3 = Pose3Pose3(odo3, odoCov)
pts0X4 = projectParticles(getVal(fg,"x3"),odo3, odoCov)
v4 = addNode!(fg,"x4", pts0X4, N=N)
addFactor!(fg,[v3;v4],odoconstr3)

# println("Plot marginals to see what is happening")
# plotKDE(marginal(getVertKDE(fg,"x3"),[1]))
# plotKDE(marginal(getVertKDE(fg,"x4"),[1]))

println("Reconstruct Bayes tree and perform inference...")
tree = wipeBuildNewTree!(fg)
inferOverTree!(fg, tree)

#