text:
# functions which work on all operating systems

include("ft_createdeviceinfolist.jl")
include("ft_getdeviceinfolist.jl")
include("ft_open.jl")
include("ft_openex.jl")
include("ft_close.jl")
include("ft_read!.jl")
include("ft_write.jl")
include("ft_setbaudrate.jl")
include("ft_setdatacharacteristics.jl")
include("ft_settimeouts.jl")
include("ft_setflowcontrol.jl")
include("ft_setdtr.jl")
include("ft_clrdtr.jl")
include("ft_setrts.jl")
include("ft_clrrts.jl")
include("ft_getmodemstatus.jl")
include("ft_getqueuestatus.jl")
include("ft_getdeviceinfo.jl")
include("ft_seteventnotification.jl")
include("ft_getstatus.jl")
include("ft_setchars.jl")
include("ft_setbreakon.jl")
include("ft_setbreakoff.jl")
include("ft_purge.jl")
include("ft_resetdevice.jl")
include("ft_stopintask.jl")
include("ft_restartintask.jl")
include("ft_setdeadmantimeout.jl")
include("ft_readee.jl")
include("ft_writeee.jl")
include("ft_eraseee.jl")
include("ft_ee_read__ft_ee_program.jl")
include("ft_ee_uasize.jl")
include("ft_ee_uaread!.jl")
include("ft_ee_uawrite.jl")
include("ft_setlatencytimer.jl")
include("ft_getlatencytimer.jl")
include("ft_setbitmode.jl")
include("ft_getbitmode.jl")
include("ft_setusbparameters.jl")


text:
"""Initialises the breaks array based on various methods:

    maxentropy: quantile / maximum entropy.
    rand: random breaks without empty classes."""
function BreaksInit(n::Int,X::Array{T,1},method::String="maxentropy") where {T<:AbstractFloat}

    if method == "maxentropy"

        np = length(X) ÷ n       # number of data points per chunk

        # breaks are determind by the lower bounds index of the data array X
        breaks = Array{Int,1}(undef,n+1)

        # the lower bounds index of the first interval is always 1
        # quantile-approach / maximum entropy method
        breaks[1:end-1] = collect(1:np:np*n)

        # this is redundant but makes class collection easier
        breaks[end] = length(X)+1

    elseif method == "rand"     # random initialisation

        # breaks are determind by the lower bounds index of the data array X
        breaks = Array{Int,1}(undef,n+1)
        breaks[1] = 1
        breaks[end] = length(X)+1

        breaks[2:end-1] = randintNDNC(3,length(X)-1,n-1)

    else
        throw(error("Method '$method' not defined."))
    end

    return breaks
end

"""randint with no duplicates and no consecutive numbers"""
function randintNDNC(N0::Int,N1::Int,n::Int)

    ratio = n/(N1-N0)

    if ratio > 0.33
        throw(error("Subset-set size ratio is $ratio, too likely that no solution can be found. Increase n."))
    elseif ratio > 0.2 && n > 10000
        @warn "Finding $n random numbers for a subset-set size ratio of $ratio may take a long time."
    end

    v = Array{Int,1}(undef,n)
    v[1] = rand(N0:N1)

    for i in 2:n
        r = rand(N0:N1)     # pick a random number
        while any(abs.(r .- v[1:i-1]) .<= 1)    # check that it's not duplicate or a consecutive
            r = rand(N0:N1) # if that's the case try another one.
        end
        v[i] = r
    end

    sort!(v)
    return v
end

ext: jl


text:
# Create sticky REPL functions for shell and help modes
const REPL_MODES_METADATA = @compat Dict(
    :shell => 2,
    :help => 3
)

for (repl_mode, mode_number) in REPL_MODES_METADATA
    sticky_mode = @compat Symbol(:sticky_, repl_mode, :_mode)
    @eval function $(sticky_mode)(enable::Bool = true)
        Base.active_repl.interface.modes[$mode_number].sticky = enable
    end
end

repl_prompt(prompt::String = "julia> ") = Base.active_repl.interface.modes[1].prompt = prompt


text:
#!/usr/bin/env julia

using WriteVTK
using Test
import Compat: UTF8String, readstring
using SHA: sha1

const tests = ["multiblock.jl",
               "rectilinear.jl",
               "imagedata.jl",
               "structured.jl",
               "unstructured.jl",
               "pvdCollection.jl",
               "array.jl"]

# Only toggle to generate new checksums, if new tests are added.
const OVERWRITE_CHECKSUMS = true
const checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
const checksum_list = read(checksums_file, String)

if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
end

ProjDir = dirname(@__FILE__)
cd(ProjDir*"/output") do

  # Run the test scripts.
  for test in tests
      println("TEST (first run): ", test)
      outfiles = evalfile(test)::Vector{UTF8String}

      # Check that the generated files match the stored checksums.
      for file in outfiles
          sha_str = bytes2hex(open(sha1, file)) * "  $file\n"
          if OVERWRITE_CHECKSUMS
              write(csio, sha_str)
          else
              # Returns 0:-1 if string is not found.
              cmp = search(checksum_list, sha_str)
              @test cmp != 0:-1
          end
      end
      println()
  end

  OVERWRITE_CHECKSUMS && close(csio)

  println("="^60, "\n")

  # Run the tests again, just to measure the time and allocations once all the
  # functions have already been compiled.
  for test in tests
      println("TEST (second run): ", test)
      outfiles = evalfile(test)::Vector{UTF8String}
      println()
  end

end


text:
module FlexibleIntervals

import Base: hash, convert, promote_rule,
    string, show, showcompact, showall,
    isnull,
    isnan, isinf, isfinite, issubnormal,
    sign, signbit, copysign, flipsign,
    ispow2, prevpow2, nextpow2, frexp, ldexp,
    abs, sqrt, abs2, hypot,
    isequal, isless,
    (==), (!=), (<), (<=), (>=), (>),
    (+),(-),(*),(/),(\),(%),(^),(//),
    fma,
    log, log1p, log2, log10, exp, expm1,
    sin, cos, tan, csc, sec, cot,
    sinh, cosh, tanh, csch, sech, coth,
    asin, acos, atan, acsc, asec, acot, atan2,
    asinh, acosh, atanh, acsch, asech, acoth

import Genera: Float, StdFloat

export show, showcompact, showall,  # sometimes necessary (q.v. Nemo/arb)
    isnull,
    ispos, isneg, isnormal, ntneg, ntnan, ntnormal,
    ispow2, prevpow2, nextpow2, frexp, ldexp,
    abs, sqrt, abs2, hypot,
    sqr, fma, fms, fam, fsm

include("type/FlexibleInterval.jl")
include("type/io.jl")
include("type/generous.jl")
include("type/predicates.jl")
include("type/ordering.jl")

end


text:
function testsameness(a, b)
    try
        @test is(a, b)
    catch
        println(a, "::", typeof(a), " is not ", b, "::", typeof(b))
        rethrow()
    end
end

function testdberror(dbobj, expected)
    try
        @test errcode(dbobj) == expected
    catch
        println(errstring(dbobj))
        println(bytestring(PostgreSQL.PQresStatus(expected)))
        rethrow()
    end
end

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using MyModule
using Plots
gr()

a = 1.0
b = 2.0
c = [1. ,2. ,3.]
Data = [1., 2., 3, 4., 5., 6.]
Num =0.
u = 0.0:0.1:1.0
v = 8. .*u

function mul(x,y)
    z = x*y
    return z
end

##別ファイルで定義した関数が使えるか
println(plus(a, b))
### →使える

##"."で要素ごとの演算ができる性質は維持されるのか
println(plus.(c, b))
### →維持される

##同じファイルで複数の関数を定義した際にそれぞれを使えるか
println(minus(a,b))
### →使える

##元ファイルで定義した変数を関数ファイルで使えるか
Print()
### →使える

##関数内で、元ファイルの変数を上書きできるか
println(Num)
OverWrite()
println(Num)
### →できない

##元ファイルで定義した関数を関数ファイルで使えるか
println(Pow(b, 5))
### →使える

##元ファイルで導入したモジュールが関数ファイルで使えるか
Pl(u,v)
### →使えない(多分?)

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
using CUDA.CUTENSOR
using CUDA
using LinearAlgebra

@test has_cutensor()
@test CUTENSOR.version() isa VersionNumber

@testset "CuTensor type basics" begin
    N = 2
    dmax = 2^div(18,N)
    dims = rand(2:dmax, N)
    p = randperm(N)
    indsA = collect(('a':'z')[1:N])
    dimsA = dims
    A = rand(Float64, dimsA...)
    dA = CuArray(A)
    p = randperm(N)
    indsA = collect(('a':'z')[1:N])
    ctA = CuTensor(dA, indsA)
    @test length(ctA) == length(A)
    @test size(ctA) == size(A)
    @test size(ctA, 1) == size(A, 1)
    @test ndims(ctA) == ndims(A)
    @test strides(ctA) == strides(A)
    @test eltype(ctA) == eltype(A)
end


text:
function getproposal(p::Prior, nparams)
  newparams = zeros(Float64, nparams)
  update_newparams!(newparams, p)
  return newparams
end

update_newparams!(newparams,p::Prior) = update_newparams!(newparams, 1, p.distribution...)

@inline function update_newparams!(newparams, i, x, y...)
  newparams[i] = rand(x)
  update_newparams!(newparams, i + 1, y...)
end
@inline function update_newparams!(newparams,i,x)
  newparams[i] = rand(x)
end

function perturbparticle(particle, kernel::Kernel)
  newparticle = copyparticle(particle)
  newparams = zeros(Float64, length(newparticle.params))
  for i in 1:length(newparams)
    newparams[i] = kernel.perturbation_function(newparticle.params[i], kernel.kernel_parameters[i])
  end
  newparticle.params = newparams
  return newparticle
end

function perturbmodel(ABCsetup, mstar, modelprob)
    prob = ABCsetup.modelkern
    mprob = ones(Float64, length(modelprob))
    mprob[modelprob.==0.0] .= 0.0
    nsurvivingmodels = sum(mprob)
    mprob[mprob.> 0.0] .= (1 - prob) / (nsurvivingmodels - 1)
    mprob[mstar] = prob
    wsample(1:ABCsetup.nmodels, mprob)
end


text:
struct Queen
    x::Int
    y::Int
end

qhorz(qa, qb) = qa.x == qb.x
qvert(qa, qb) = qa.y == qb.y
qdiag(qa, qb) = abs(qa.x - qb.x) == abs(qa.y - qb.y)

qhvd(qa, qb) = qhorz(qa, qb) || qvert(qa, qb) || qdiag(qa, qb)
qany(testq, qs) = any(q -> qhvd(testq, q), qs)

function qsolve(nsqsx, nsqsy, nqs, presqs = ())
    nqs == 0 && return presqs
    for xsq in 1:nsqsx
        for ysq in 1:nsqsy
            testq = Queen(xsq, ysq)
            if !qany(testq, presqs)
                tryqs = (presqs..., testq)
                maybe = qsolve(nsqsx, nsqsy, nqs - 1, tryqs)
                maybe !== nothing && return maybe
            end
        end
    end
    return nothing
end

qsolve(nqs) = qsolve(nqs, nqs, nqs, presqs = ())


text:
#!/usr/bin/env julia

using Luxor, Test, Random

Random.seed!(42)

function test_styles(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    background("antiquewhite")
    ds = Style()
    dark = Style(colorant"orange")
    jm = Style(colorant"green", "JuliaMono-Black", 20)
    applystyle(ds)
    box(O, 140, 140, :fill)
    applystyle(dark)
    box(O, 40, 240, :fill)
    text("hello world")
    applystyle(jm)
    text("hello world again", O + (20, 30))
    @test finish() == true
    println("...finished styles test, saved in $(fname)")
end
fname = "styles-test.png"
test_styles(fname)


text:
#=
Julia port of Andrew Patterson's PyExpUtils file permutations.py

https://github.com/andnp/PyExpUtils/blob/master/PyExpUtils/utils/permute.py
=#

function flattenToArray(thing::Dict)
    accum = []

    function inner(thing, path)
        if isa(thing, Array)
            if isa(thing[1], Dict)
                i=1
                for sub in thing
                    inner(sub, string(path,".[", i, "]"))
                    i+=1
                end
                return
            end

            push!(accum, (path, thing))
            return
        end

        if isa(thing, Dict)
            for key in keys(thing)
                new_path = isempty(path) ? key : string(path,".",key)
                inner(thing[key], new_path)
            end
            return
        end
        push!(accum, (path, [thing]))
        return
    end

    inner(thing, "")
    return accum
end

function getTotalCombinations(pairs::Vector{Any})
    accum = 1
    for (path, values) in pairs
        num = length(values) > 0 ? length(values) : 1
        accum *= num
    end
    return accum
end

getTotalCombinations(cfg::Dict) = getTotalCombinations(flattenToArray(cfg))

function getPermutation(pairs::Vector{Any}, idx::Int)
    perm = Dict()
    accum = 1

    for (key, values) in pairs
        num = length(values)

        if num == 0
            perm[key] = []
            continue
        end

        value_idx = Int(floor((idx-1)/accum)) % num
        perm[key] = values[value_idx + 1]
        accum *= num
    end

    return reconstructParameters(perm)
end

getPermutation(cfg::Dict, idx::Int) = getPermutation(flattenToArray(cfg), idx)

function reconstructParameters(perm::Dict)
    res = Dict()
    for key in keys(perm)
        set_at_path(res, key, perm[key])
    end
    return res
end

function set_at_path(d::Dict, path::String, val::Any)
    function inner(d::Dict, path::String, val, last)
        if length(path) == 0
            return d
        end
        sp = split(path, "."; limit=2)
        sp = map(s->String(s), sp)

        part, rest = length(sp) > 1 ? sp : (sp[1], "")
        nxt = split(rest, ".")[1]

        # lists
        if startswith(part, "[")
            num = parse(Int, replace(part, r"[\[,\]]" =>s""))

            if length(d[last]) > num
                piece = length(rest) > 0 ? inner(d[last][num], rest, val, "") : val
                d[last][num] = piece
            else
                piece = length(rest) > 0 ? inner(Dict(), rest, val, "") : val
                push!(d[last], piece)
            end
            return d

        # objects
        elseif length(rest) > 0
            if startswith(nxt, "[")
                piece = get!(d, part, [])
                return inner(d, rest, val, part)
            else
                piece = get!(d, part, Dict())
                return inner(piece, rest, val, part)
            end
        # everything else
        else
            get!(d, part, val)
            return d
        end
    end

    inner(d, path, val, "")
    return d
end

function getDiff!(dict::Dict, cfg_dict::Dict)
    dkeys = []
    for (k,v) in cfg_dict
        if typeof(dict[k]) <: Dict
            subdiff = getDiff(dict[k], cfg_dict[k])
            if isempty(subdiff)
                push!(dkeys, k)
            end
        elseif all(v.==dict[k])
            push!(dkeys, k)
        end
    end

    map(k->delete!(cfg_dict, k), dkeys)

    return cfg_dict
end

getDiff(dict::Dict, cfg_dict::Dict) = getDiff!(dict, copy(cfg_dict))

function _ensure(dirs::Vector{String})
    for d in dirs
        if !isdir(d)
            mkpath(d)
        end
    end
end


text:
# real U
@inline function subspace_mul_kernel_complex!(S::AbstractMatrix{T}, (y_re, y_im), indices, U::AbstractMatrix{T}, k::Int, offset::Int) where {T <: Base.HWReal}
    @avx for i in axes(U, 1)
        y_re_i = zero(T)
        y_im_i = zero(T)
        for j in axes(U, 2)
            idx_j = indices[j] + k + offset
            y_re_i += U[i, j] * S[1, idx_j]
            y_im_i += U[i, j] * S[2, idx_j]
        end
        y_re[i] = y_re_i
        y_im[i] = y_im_i
    end

    @avx for i in axes(U, 1)
        idx_i = indices[i] + k + offset
        S[1, idx_i] = y_re[i]
        S[2, idx_i] = y_im[i]
    end
end

# Diagonal real U
@inline function subspace_mul_kernel_complex!(S::AbstractMatrix{T}, (y_re, y_im), indices,
    U::Diagonal{T}, k::Int, offset::Int) where {T <: Base.HWReal}

    @avx for i in axes(U, 1)
        idx_i = indices[i] + k + offset
        S[1, idx_i] = U[i, i] * S[1, idx_i]
        S[2, idx_i] = U[i, i] * S[2, idx_i]
    end
end

# complex U
@inline function subspace_mul_kernel_complex!(S::AbstractMatrix{T}, (y_re, y_im), indices, (U_re, U_im), k::Int, offset::Int) where {T <: Base.HWReal}
    @avx for i in axes(U_re, 1)
        y_re_i = zero(T)
        y_im_i = zero(T)
        for j in axes(U_re, 2)
            idx_j = indices[j] + k + offset
            y_re_i += U_re[i, j] * S[1, idx_j] - U_im[i, j] * S[2, idx_j]
            y_im_i += U_re[i, j] * S[2, idx_j] + U_im[i, j] * S[1, idx_j]
        end
        y_re[i] = y_re_i
        y_im[i] = y_im_i
    end

    @avx for i in axes(U_re, 1)
        idx_i = indices[i] + k + offset
        S[1, idx_i] = y_re[i]
        S[2, idx_i] = y_im[i]
    end
end

# Diagonal complex U (not split up by split_op)
@inline function subspace_mul_kernel_complex!(S::AbstractMatrix{T}, (y_re, y_im), indices,
    U::Diagonal{Complex{T}}, k::Int, offset::Int) where {T <: Base.HWReal}

    @avx for i in axes(U, 1)
        idx_i = indices[i] + k + offset
        S[1, idx_i], S[2, idx_i] =
            real(U[i, i]) * S[1, idx_i] - imag(U[i, i]) * S[2, idx_i],
            real(U[i, i]) * S[2, idx_i] + imag(U[i, i]) * S[1, idx_i]
    end
end

# real U, multiple states in S
@inline function subspace_mul_kernel_complex!(S::AbstractArray{T, 3}, (C_re, C_im), indices, U::AbstractMatrix{T}, k::Int, _b::Int, Bmax::Int, offset::Int) where {T <: Base.HWReal}
    b = _b << 3;
    bmax = b + 8
    if bmax ≤ Bmax # full static block
        @avx for n ∈ 1:8, m ∈ axes(U, 1)
            C_re_m_n = zero(T)
            C_im_m_n = zero(T)
            for i ∈ axes(U, 2)
                j = k + indices[i] + offset
                C_re_m_n += U[m,i] * S[1,n+b,j]
                C_im_m_n += U[m,i] * S[2,n+b,j]
            end
            C_re[m,n] = C_re_m_n
            C_im[m,n] = C_im_m_n
        end
        @avx for n ∈ 1:8, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1,n+b,S_m] = C_re[m,n]
            S[2,n+b,S_m] = C_im[m,n]
        end
        # AmulB!(C_re, C_im, U, U_im,
    else # dynamic block
        Nmax = 8 + Bmax - bmax
        @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
            C_re_m_n = zero(T)
            C_im_m_n = zero(T)
            for i ∈ axes(U, 2)
                j = k + indices[i] + offset
                C_re_m_n += U[m,i] * S[1,n+b,j]
                C_im_m_n += U[m,i] * S[2,n+b,j]
            end
            C_re[m,n] = C_re_m_n
            C_im[m,n] = C_im_m_n
        end
        @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1,n+b,S_m] = C_re[m,n]
            S[2,n+b,S_m] = C_im[m,n]
        end
    end
    return
end

# Diagonal real U, multiple states in S
@inline function subspace_mul_kernel_complex!(S::AbstractArray{T, 3}, (C_re, C_im), indices,
    U::Diagonal{T}, k::Int, _b::Int, Bmax::Int, offset::Int) where {T <: Base.HWReal}

    b = _b << 3;
    bmax = b + 8
    if bmax ≤ Bmax # full static block
        @avx for n ∈ 1:8, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1, n+b, S_m] = U[m, m] * S[1, n+b, S_m]
            S[2, n+b, S_m] = U[m, m] * S[2, n+b, S_m]
        end
        # AmulB!(C_re, C_im, U, U_im,
    else # dynamic block
        Nmax = 8 + Bmax - bmax
        @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1, n+b, S_m] = U[m, m] * S[1, n+b, S_m]
            S[2, n+b, S_m] = U[m, m] * S[2, n+b, S_m]
        end
    end
    return
end

# complex U, multiple states in S
@inline function subspace_mul_kernel_complex!(S::AbstractArray{T, 3}, (C_re, C_im), indices, (U_re, U_im), k::Int, _b::Int, Bmax::Int, offset::Int) where {T <: Base.HWReal}
    b = _b << 3;
    bmax = b + 8
    if bmax ≤ Bmax # full static block
        @avx for n ∈ 1:8, m ∈ axes(U_re, 1)
            C_re_m_n = zero(T)
            C_im_m_n = zero(T)
            for i ∈ axes(U_re, 2)
                j = k + indices[i] + offset
                C_re_m_n += U_re[m,i] * S[1,n+b,j] - U_im[m,i] * S[2,n+b,j]
                C_im_m_n += U_re[m,i] * S[2,n+b,j] + U_im[m,i] * S[1,n+b,j]
            end
            C_re[m,n] = C_re_m_n
            C_im[m,n] = C_im_m_n
        end
        @avx for n ∈ 1:8, m ∈ axes(U_re, 1)
            S_m = k + indices[m] + offset
            S[1,n+b,S_m] = C_re[m,n]
            S[2,n+b,S_m] = C_im[m,n]
        end
        # AmulB!(C_re, C_im, U_re, U_im,
    else # dynamic block
        Nmax = 8 + Bmax - bmax
        @avx for n ∈ 1:Nmax, m ∈ axes(U_re, 1)
            C_re_m_n = zero(T)
            C_im_m_n = zero(T)
            for i ∈ axes(U_re, 2)
                j = k + indices[i] + offset
                C_re_m_n += U_re[m,i] * S[1,n+b,j] - U_im[m,i] * S[2,n+b,j]
                C_im_m_n += U_re[m,i] * S[2,n+b,j] + U_im[m,i] * S[1,n+b,j]
            end
            C_re[m,n] = C_re_m_n
            C_im[m,n] = C_im_m_n
        end
        @avx for n ∈ 1:Nmax, m ∈ axes(U_re, 1)
            S_m = k + indices[m] + offset
            S[1,n+b,S_m] = C_re[m,n]
            S[2,n+b,S_m] = C_im[m,n]
        end
    end
    return
end

# Diagonal complex U, multiple states in S
@inline function subspace_mul_kernel_complex!(S::AbstractArray{T, 3}, (C_re, C_im), indices,
    U::Diagonal{Complex{T}}, k::Int, _b::Int, Bmax::Int, offset::Int) where {T <: Base.HWReal}

    b = _b << 3;
    bmax = b + 8
    if bmax ≤ Bmax # full static block
        @avx for n ∈ 1:8, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1, n+b, S_m], S[2, n+b, S_m] =
                real(U[m, m]) * S[1, n+b, S_m] - imag(U[m, m]) * S[2, n+b, S_m],
                real(U[m, m]) * S[2, n+b, S_m] + imag(U[m, m]) * S[1, n+b, S_m]
        end
        # AmulB!(C_re, C_im, U, U_im,
    else # dynamic block
        Nmax = 8 + Bmax - bmax
        @avx for n ∈ 1:Nmax, m ∈ axes(U, 1)
            S_m = k + indices[m] + offset
            S[1, n+b, S_m], S[2, n+b, S_m] =
                real(U[m, m]) * S[1, n+b, S_m] - imag(U[m, m]) * S[2, n+b, S_m],
                real(U[m, m]) * S[2, n+b, S_m] + imag(U[m, m]) * S[1, n+b, S_m]
        end
    end
    return
end


text:
using Random
using SHTOOLS
using Test

include("test_legendre.jl")
include("test_transforms.jl")
include("test_storage.jl")


text:
const FN = typeof(identity)
const FC = typeof(conj)
const FA = typeof(adjoint)
const FT = typeof(transpose)
_conj(::FN) = conj
_conj(::FC) = identity
_conj(::FA) = transpose
_conj(::FT) = adjoint
_transpose(::FN) = transpose
_transpose(::FC) = adjoint
_transpose(::FA) = conj
_transpose(::FT) = identity
_adjoint(::FN) = adjoint
_adjoint(::FC) = transpose
_adjoint(::FA) = identity
_adjoint(::FT) = conj

abstract type AbstractStridedView{T,N,F<:Union{FN,FC,FA,FT}} <: AbstractArray{T,N} end

Base.elsize(::Type{<:AbstractStridedView{T}}) where {T} =
    Base.isbitstype(T) ? sizeof(T) :
        (Base.isbitsunion(T) ? Base.bitsunionsize(T) : sizeof(Ptr))

# Converting back to other DenseArray type:
function Base.convert(T::Type{<:DenseArray}, a::AbstractStridedView)
    b = T(undef, size(a))
    copyto!(StridedView(b), a)
    return b
end
# following method because of ambiguity warning
function Base.convert(::Type{T}, a::AbstractStridedView) where {T<:Array}
    b = T(undef, size(a))
    copyto!(StridedView(b), a)
    return b
end
function Base.convert(::Type{Array}, a::AbstractStridedView{T}) where {T}
    b = Array{T}(undef, size(a))
    copyto!(StridedView(b), a)
    return b
end

# Methods based on map!
Base.copyto!(dst::AbstractStridedView{<:Any,N}, src::AbstractStridedView{<:Any,N}) where {N} =
    map!(identity, dst, src)
Base.conj!(a::AbstractStridedView{<:Real}) = a
Base.conj!(a::AbstractStridedView) = map!(conj, a, a)
LinearAlgebra.adjoint!(dst::AbstractStridedView{<:Any,N}, src::AbstractStridedView{<:Any,N}) where {N} =
    copyto!(dst, adjoint(src))
Base.permutedims!(dst::AbstractStridedView{<:Any,N}, src::AbstractStridedView{<:Any,N}, p) where {N} =
    copyto!(dst, permutedims(src, p))

# linear algebra
LinearAlgebra.rmul!(dst::AbstractStridedView, α::Number) = mul!(dst, dst, α)
LinearAlgebra.lmul!(α::Number, dst::AbstractStridedView) = mul!(dst, α, dst)

function LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N}, α::Number, src::AbstractStridedView{<:Number,N}) where {N}
    if α == 1
        copyto!(dst, src)
    else
        dst .= α .* src
    end
    return dst
end
function LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N}, src::AbstractStridedView{<:Number,N}, α::Number) where {N}
    if α == 1
        copyto!(dst, src)
    else
        dst .= src .* α
    end
    return dst
end
function LinearAlgebra.axpy!(a::Number, X::AbstractStridedView{<:Number,N}, Y::AbstractStridedView{<:Number,N}) where {N}
    if a == 1
        Y .= X .+ Y
    else
        Y .= a .* X .+ Y
    end
    return Y
end
function LinearAlgebra.axpby!(a::Number, X::AbstractStridedView{<:Number,N}, b::Number, Y::AbstractStridedView{<:Number,N}) where {N}
    if b == 1
        axpy!(a, X, Y)
    elseif b == 0
        mul!(Y, a, X)
    else
        Y .= a .* X .+ b .* Y
    end
    return Y
end

function LinearAlgebra.mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{<:Any,2}, B::AbstractStridedView{<:Any,2}, α::Number = true, β::Number = false) where {T}
    if !(eltype(C) <: LinearAlgebra.BlasFloat && eltype(A) == eltype(B) == eltype(C))
        return __mul!(C, A, B, α, β)
    end
    # C.op is identity or conj
    if C.op == conj
        if stride(C,1) < stride(C,2)
            _mul!(conj(C), conj(A), conj(B), conj(α), conj(β))
        else
            _mul!(C', B', A', conj(α), conj(β))
        end
    elseif stride(C,1) > stride(C,2)
        _mul!(transpose(C), transpose(B), transpose(A), α, β)
    else
        _mul!(C, A, B, α, β)
    end
    return C
end

function isblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}
    if A.op == identity
        return stride(A,1) == 1 || stride(A,2) == 1
    elseif A.op == conj
        return stride(A, 2) == 1
    else # should never happen
        return false
    end
end
function getblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}
    if A.op == identity
        if stride(A,1) == 1
            return blasstrides(A), 'N'
        else
            return blasstrides(transpose(A)), 'T'
        end
    else
        return blasstrides(adjoint(A)), 'C'
    end
end

# here we will have C.op == :identity && stride(C,1) < stride(C,2)
function _mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{T,2}, B::AbstractStridedView{T,2}, α, β) where {T<:LinearAlgebra.BlasFloat}
    if stride(C,1) == 1 && isblasmatrix(A) && isblasmatrix(B)
        nthreads = use_threaded_mul() ? get_num_threads() : 1
        _threaded_blas_mul!(C, A, B, α, β, nthreads)
    else
        return __mul!(C, A, B, α, β)
    end
end

function _threaded_blas_mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{T,2}, B::AbstractStridedView{T,2}, α, β, nthreads) where {T<:LinearAlgebra.BlasFloat}
    m, n = size(C)
    m == size(A, 1) && n == size(B, 2) || throw(DimensionMismatch())
    if nthreads == 1 || m*n < 1024
        A2, CA = getblasmatrix(A)
        B2, CB = getblasmatrix(B)
        C2 = blasstrides(C)
        LinearAlgebra.BLAS.gemm!(CA, CB, convert(T, α), A2, B2, convert(T, β), C2)
    else
        if m > n
            m2 = round(Int, m/16)*8
            nthreads2 = nthreads >> 1
            t = Threads.@spawn _threaded_blas_mul!(C[1:m2, :], A[1:m2, :], B, α, β, nthreads2)
            _threaded_blas_mul!(C[m2+1:m, :], A[m2+1:m, :], B, α, β, nthreads - nthreads2)
            wait(t)
            return C
        else
            n2 = round(Int, n/16)*8
            nthreads2 = nthreads >> 1
            t = Threads.@spawn _threaded_blas_mul!(C[:, 1:n2], A, B[:, 1:n2], α, β, nthreads2)
            _threaded_blas_mul!(C[:, n2+1:n], A, B[:, n2+1:n], α, β, nthreads - nthreads2)
            wait(t)
            return C
        end
    end
end

# This implementation is faster than LinearAlgebra.generic_matmatmul
function __mul!(C::AbstractStridedView{<:Any,2}, A::AbstractStridedView{<:Any,2}, B::AbstractStridedView{<:Any,2}, α, β)
    (size(C,1) == size(A,1) && size(C,2) == size(B,2) && size(A,2) == size(B,1)) ||
        throw(DimensionMatch("A has size $(size(A)), B has size $(size(B)), C has size $(size(C))"))

    m,n = size(C)
    k = size(A,2)
    A2 = sreshape(A, (m, 1, k))
    B2 = sreshape(permutedims(B,(2,1)), (1, n, k))
    C2 = sreshape(C, (m, n, 1))

    if α == 0 || k == 0
        rmul!(C, β)
    elseif α == 1
        if β == 0
            _mapreducedim!(*, +, zero, (m,n,k), (C2,A2,B2))
        elseif β == 1
            _mapreducedim!(*, +, nothing, (m,n,k), (C2,A2,B2))
        else
            _mapreducedim!(*, +, x->x*β, (m,n,k), (C2,A2,B2))
        end
    else
        f = (x,y)->(x*y*α)
        if β == 0
            _mapreducedim!(f, +, zero, (m,n,k), (C2,A2,B2))
        elseif β == 1
            _mapreducedim!(f, +, nothing, (m,n,k), (C2,A2,B2))
        else
            _mapreducedim!(f, +, x->x*β, (m,n,k), (C2,A2,B2))
        end
    end
    return C
end

# ParentIndex: index directly into parent array
struct ParentIndex
    i::Int
end

function sreshape end

sreshape(a::AbstractStridedView,args::Vararg{Int}) = sreshape(a, args)

struct ReshapeException <: Exception
end
Base.show(io::IO, e::ReshapeException) = print(io, "Cannot produce a reshaped StridedView without allocating, try sreshape(copy(array), newsize) or fall back to reshape(array, newsize)")

@inline sview(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(a, I...)
@inline sview(a::AbstractStridedView, I::SliceIndex) =
    getindex(sreshape(a, (length(a),)), I...)

@inline Base.view(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(a, I...)

@inline sview(a::DenseArray{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(StridedView(a), I...)
@inline sview(a::DenseArray, I::SliceIndex) =
    getindex(sreshape(StridedView(a), (length(a),)), I...)

# Auxiliary routines
@inline _computeind(indices::Tuple{}, strides::Tuple{}) = 1
@inline _computeind(indices::NTuple{N,Int}, strides::NTuple{N,Int}) where {N} =
    (indices[1]-1)*strides[1] + _computeind(tail(indices), tail(strides))

@inline function _computeviewsize(oldsize::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Int)
        return _computeviewsize(tail(oldsize), tail(I))
    elseif isa(I[1], Colon)
        return (oldsize[1], _computeviewsize(tail(oldsize), tail(I))...)
    else
        return (length(I[1]), _computeviewsize(tail(oldsize), tail(I))...)
    end
end
_computeviewsize(::Tuple{}, ::Tuple{}) = ()

@inline function _computeviewstrides(oldstrides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Int)
        return _computeviewstrides(tail(oldstrides), tail(I))
    elseif isa(I[1], Colon)
        return (oldstrides[1], _computeviewstrides(tail(oldstrides), tail(I))...)
    else
        return (oldstrides[1]*step(I[1]), _computeviewstrides(tail(oldstrides), tail(I))...)
    end
end
_computeviewstrides(::Tuple{}, ::Tuple{}) = ()

@inline function _computeviewoffset(strides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Colon)
        return _computeviewoffset(tail(strides), tail(I))
    else
        return (first(I[1])-1)*strides[1]+_computeviewoffset(tail(strides), tail(I))
    end
end
_computeviewoffset(::Tuple{}, ::Tuple{}) = 0

_simplify(size::Tuple{}, strides::Tuple{}) = size, strides
_simplify(size::Dims{1}, strides::Dims{1}) = size, strides
function _simplify(size::Dims{N}, strides::Dims{N}) where {N}
    tailsize, tailstrides = _simplify(tail(size), tail(strides))
    if size[1] == 1
        return (tailsize..., 1), (tailstrides..., strides[1])
    elseif size[1]*strides[1] == tailstrides[1]
        return (size[1]*tailsize[1], tail(tailsize)..., 1),
            (strides[1], tail(tailstrides)..., tailsize[1]*tailstrides[1])
    else
        return (size[1], tailsize...), (strides[1], tailstrides...)
    end
end

_normalizestrides(size::Tuple{}, strides::Tuple{}) = strides
function _normalizestrides(size::Dims{N}, strides::Dims{N}) where {N}
    if size[1] == 1
        strides = Base.setindex(strides, 1, 1)
    elseif size[1] == 0
        strides = one.(strides)
    end
    for i = 2:N
        if size[i] == 1
            strides = Base.setindex(strides, size[i-1]*strides[i-1], i)
        elseif size[i] == 0
            for j = i:N
                strides = Base.setindex(strides, size[i-1]*strides[i-1], j)
            end
        end
    end
    return strides
end

_computereshapestrides(newsize::Tuple{}, oldsize::Tuple{}, strides::Tuple{}) = strides
function _computereshapestrides(newsize::Tuple{}, oldsize::Dims{N}, strides::Dims{N}) where {N}
    all(isequal(1), oldsize) || throw(DimensionMismatch())
    return ()
end

function _computereshapestrides(newsize::Dims, oldsize::Tuple{}, strides::Tuple{})
    all(isequal(1), newsize) || throw(DimensionMismatch())
    return newsize
end
function _computereshapestrides(newsize::Dims, oldsize::Dims{N}, strides::Dims{N}) where {N}
    d,r = divrem(oldsize[1], newsize[1])
    if r == 0
        s1 = strides[1]
        if d == 1
            oldsize = (tail(oldsize)..., 1)
            strides = (tail(strides)..., newsize[1]*s1)
            return (s1, _computereshapestrides(tail(newsize), oldsize, strides)...)
        else
            oldsize = (d, tail(oldsize)...)
            strides = (newsize[1]*s1, tail(strides)...)
            return (s1, _computereshapestrides(tail(newsize), oldsize, strides)...)
        end
    else
        if prod(newsize) != prod(oldsize)
            throw(DimensionMismatch())
        else
            throw(ReshapeException())
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractProblem

A generic problem in geostatistics.
"""
abstract type AbstractProblem end

#------------------
# IMPLEMENTATIONS
#------------------
include("problems/estimation.jl")
include("problems/simulation.jl")
include("problems/learning.jl")


text:
# # Using different kernels
#
# Starting with version v0.6, two new kernels can be chosen: a **linear** kernel and a
# **polynomial** kernel.
# In this document we will see how to handle choosing a different kernel, and we'll
# showcase an example.
#
# First, we need to import all the necessary packages.
using LeastSquaresSVM
using MLJ, MLJBase
using DataFrames, CSV
using CategoricalArrays, Random
using Plots
gr();
rng = MersenneTwister(812);

# For this example, we will create a large classification problem. It is actually
# inspired from a similar classification problem from [`scikit-learn`](https://scikit-learn.org/stable/auto_examples/svm/plot_separating_hyperplane_unbalanced.html#sphx-glr-auto-examples-svm-plot-separating-hyperplane-unbalanced-py).
#
# The idea is to have a very large number of features (5000), and a small number of
# instances.
# This has been [reported](https://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf) as
# as being a good use case or rule of thumb
#
# > Whenever the number of features is _larger_ than the number of instances, use a
# > linear kernel.
#
X, y = MLJ.make_blobs(500, 2_000; centers=2, cluster_std=[1.5, 0.5]);

# Of course, this is just to showcase the implementation within `LeastSquaresSVM`. There are
# actually better ways to handle this kind of problem, e.g. dimensionality-reduction
# algorithms.
#
# The `make_blobs` function is very similar to that of [`scikit-learn`s](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_blobs.html#sklearn.datasets.make_blobs).
# The idea is to create circles, or clusters, and to classify between them.
#
# We'll use them to test our _linear_ kernel.

# We need to construct a `DataFrame` with the arrays created to better handle the data,
# as well as a better integration with `MLJ`.
df = DataFrame(X);
df.y = y;

# Recall that we need to change the primitive types of `Julia` to `scitypes`.
dfnew = coerce(df, autotype(df));

# We can then observe the first three columns, together with their new types.
# We'll just look at the first 8 features to avoid cluttering the space.
first(dfnew[:, 1:8], 3) |> pretty

# We should also check out the basic statistics of the dataset. We'll only use a small
# subset as the data frame it quite large.
describe(dfnew[1:20, 1:10], :mean, :std, :eltype)

# Recall that we also need to standardize the dataset, we can see here that the mean is
# close to zero, but not quite, and we also need an unitary standard deviation.

# Split the dataset into training and testing sets.
y, X = unpack(dfnew, ==(:y), colname -> true);
train, test = partition(eachindex(y), 0.75, shuffle=true, rng=rng);

# In this document, we will use a _pipeline_ to integrate both the standardization and
# the LS-SVM classifier into one step. This should make it easier to train and to also
# showcase the ease of use of the `@pipeline` macro from `MLJ.jl`.
#
# First, we define our classifier and the hyperparameter range. The `self_tuning_model`
# variable will hold the model with the best performance.
#
# For the case of a _linear_ kernel, no hyperparameter is needed. Instead, the only
# hyperparameter that needs to be adjusted is the ``\sigma`` value that is intrinsic
# of the least-squares formulation. We will search for a good hyperparameter now.
#
# We will use the `accuracy` as a metric. The accuracy is simply defined as
#
# ```math
# \text{accuracy} = \frac{\text{number of correct predictions}}{\text{total number of predictions}}
# ```
#
# Note that the accuracy is not always a good measure of classification, but it will do
# fine on this dataset.
#
model = LSSVClassifier(kernel=:linear);
r1 = range(model, :σ, lower=1.0, upper=1000.0);
self_tuning_model = TunedModel(
    model=model,
    tuning=Grid(goal=400, rng=rng),
    resampling=StratifiedCV(nfolds=5),
    range=[r1],
    measure=accuracy,
    acceleration=CPUThreads(), # We use this to enable multithreading
);

# Then, we will build the pipeline. The first step is to standardize the inputs and then
# pass it to the classifier.
pipe = @pipeline(Standardizer(), self_tuning_model);

# !!! warning
#     Remember that the least-squares formulation uses **all** the data samples, so the
#     following will actually consume at least 6 GB of RAM or more. Do not run this on
#     your hardware if you are not sure you have this kind of resources available.
#     At the very least, replace `CPUThreads()` with `CPU1()` to disable multithreading.
#     Methods to handle memory more efficiently will be available in future
#     versions.
#
# And now we proceed to train all the models and find the best one!
#
mach = machine(pipe, X, y);
fit!(mach, rows=train, verbosity=0);
fitted_params(mach).deterministic_tuned_model.best_model

# Having found the best hyperparameters for the regressor model we proceed to check how
# the model generalizes.
# To do this we use the test set to check the performance.
ŷ = MLJBase.predict(mach, rows=test);
result = accuracy(ŷ, y[test]);
result

# We can see that we did quite well. A value of 1, or close enough, means the classifier
# is _perfect._ That is, it can classify correctly between each class.
#
# Finally, let us look at the so-called _confusion matrix._ This table shows us useful
# information about the performance of our classifier.
#
# Let us compute it first, and then we'll analyse it. Notice, however, that we need to
# first coerce the types to `OrderedFactor` _scitypes_ in order for the confusion matrix
# to be computed correctly.
ŷ = coerce(ŷ, OrderedFactor);
y_ordered = coerce(y[test], OrderedFactor);
confusion_matrix(ŷ, y_ordered)

# The way you read the confusion matrix is the following. The main diagonal tells us how
# many correct predictions were obtained by the classifier for both classes.
# On the other hand, the other values are the following
#
# - The _upper right_ value is known as the **false positive**. This is the number of instances that were classified as belonging to a given class, when actually they were instances of the other one. An example would be if we had an instance ``x_1`` which belonged to the class `b`, but the classifier would have predicted it actually belonged to class `a`.
# - The _lower left_ value is known as the **false negative**. This is the number of instances classified as _not_ belonging to a given class, when they actually belonged to a class. An example would be if we had an instance ``x_2`` which belonged to class `a`, and the classifier actually predicted it belonged to class `b`.
#
# It might be a little bit confusing, so a good starting point for more information on the
# subject is the excellent [Wikipedia article.](https://en.wikipedia.org/wiki/Confusion_matrix)
# You might also be interested in the following [document](https://developers.google.com/machine-learning/crash-course/classification/accuracy) from a Google's Machine Learning Crash Course.
#
# In this case, we can see that no false negative or positive values were found, which
# means that the classifier did outstandingly good.
# Normally, we can expect to have at least some percentage of false negative or positives.
#

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
# ## Prelims
#
# This tutorial builds upon the previous ensemble tutorial with a home-made Random Forest regressor on the "boston" dataset.
#

using MLJ
using PyPlot
using PrettyPrinting
using StableRNGs
import DataFrames

MLJ.color_off() # hide
ioff() # hide
X, y = @load_boston
sch = schema(X)
p = length(sch.names)
n = sch.nrows
@show (n, p)
DataFrames.describe(y)

# Let's load the decision tree regressor

@load DecisionTreeRegressor

# Let's first check the performances of just a single Decision Tree Regressor (DTR for short):

tree = machine(DecisionTreeRegressor(), X, y)
e = evaluate!(tree, resampling=Holdout(fraction_train=0.8),
              measure=[rms, rmslp1])
e |> pprint # use PrettyPrinting

# Note that multiple measures can be reported simultaneously.

# ## Random forest
#
# Let's create an ensemble of DTR and fix the number of subfeatures to 3 for now.

forest = EnsembleModel(atom=DecisionTreeRegressor())
forest.atom.n_subfeatures = 3

# (**NB**: we could have fixed `n_subfeatures` in the DTR constructor too).
#
# To get an idea of how many trees are needed, we can follow the evaluation of the error (say the `rms`) for an increasing number of tree over several sampling round.

rng = StableRNG(5123) # for reproducibility
m = machine(forest, X, y)
r = range(forest, :n, lower=10, upper=1000)
curves = learning_curve!(m, resampling=Holdout(fraction_train=0.8, rng=rng),
                         range=r, measure=rms);

# let's plot the curves

figure(figsize=(8,6))
plot(curves.parameter_values, curves.measurements)
ylabel("Root Mean Squared error", fontsize=16)
xlabel("Number of trees", fontsize=16)
xticks([10, 250, 500, 750, 1000], fontsize=14)
yticks(fontsize=14)

savefig(joinpath(@OUTPUT, "A-ensembles-2-curves.svg")) # hide

# \figalt{RMS vs number of trees}{A-ensembles-2-curves.svg}
#
# The curve is pretty noisy but let's just go for 150 trees:

forest.n = 150;

# ### Tuning
#
# As `forest` is a composite model, it has nested hyperparameters:

params(forest) |> pprint

# Let's define a range for the number of subfeatures and for the bagging fraction:

r_sf = range(forest, :(atom.n_subfeatures), lower=1, upper=12)
r_bf = range(forest, :bagging_fraction, lower=0.4, upper=1.0);

# And build a tuned model as usual that we fit on a 80/20 split.
# We use a low-resolution grid here to make this tutorial faster but you could of course use a finer grid.

tuned_forest = TunedModel(model=forest,
                          tuning=Grid(resolution=3),
                          resampling=CV(nfolds=6, rng=StableRNG(32)),
                          ranges=[r_sf, r_bf],
                          measure=rms)
m = machine(tuned_forest, X, y)
e = evaluate!(m, resampling=Holdout(fraction_train=0.8),
              measure=[rms, rmslp1])
e |> pprint

# ### Reporting
# Again, you could show a 2D heatmap of the hyperparameters

r = report(m)

figure(figsize=(8,6))

res = r.plotting

vals_sf = res.parameter_values[:, 1]
vals_bf = res.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, res.measurements)
xticks(1:3:12, fontsize=12)
xlabel("Number of sub-features", fontsize=14)
yticks(0.4:0.2:1, fontsize=12)
ylabel("Bagging fraction", fontsize=14)

savefig(joinpath(@OUTPUT, "A-ensembles-2-heatmap.svg")) # hide

# \fig{A-ensembles-2-heatmap.svg}
#
# Even though we've only done a very rough search, it seems that around 7 sub-features and a bagging fraction of around `0.75` work well.
#
# Now that the machine `m` is trained, you can use use it for predictions (implicitly, this will use the best model).
# For instance we could look at predictions on the whole dataset:

ŷ = predict(m, X)
rms(ŷ, y)

PyPlot.close_figs() # hide


text:
# Solve reinitialisation equation to reset ϕ as signed-distance function
# Alex Tam, 09/11/2021

function reinitialisation(ϕ, par, dx, dy, nI)
    dτ = 0.5*(dx/5 + dy/5)
    for i = 1:nI
        ϕ = reinit_lt(ϕ, par, dx, dy, dτ)
    end
    return ϕ
end

"Compute one τ-step of the reinitialisation equation using Lin-Tadmor"
function reinit_lt(ϕ, par, dx, dy, dt)
    ϕn = Array{Float64}(undef, par.Nx, par.Ny) # Pre-allocate array of ϕ
    dxϕ = Array{Float64}(undef, par.Nx, par.Ny) # Pre-allocate x-derivatives at non-staggered grid points
    dyϕ = Array{Float64}(undef, par.Nx, par.Ny) # Pre-allocate y-derivatives at non-staggered grid points
    S = Array{Float64}(undef, par.Nx, par.Ny) # Pre-allocate array of S
    # Compute derivatives at non-staggered points at τ_n
    for i = 1:par.Nx
        for j = 1:par.Ny # Loop over grid points
            # dϕ/dx
            if i == 1 # Boundary condition
                dxϕ[i,j] = (ϕ[i+1,j] - ϕ[i,j])/dx
            elseif i == par.Nx # Boundary condition
                dxϕ[i,j] = (ϕ[i,j] - ϕ[i-1,j])/dx
            else # Apply min-mod at interior grid points
                dxϕ[i,j] = minmod(par.θ*(ϕ[i+1,j]-ϕ[i,j])/dx, (ϕ[i+1,j]-ϕ[i-1,j])/(2*dx), par.θ*(ϕ[i,j]-ϕ[i-1,j])/dx)
            end
            # dϕ/dy
            if j == 1 # Boundary condition
                dyϕ[i,j] = (ϕ[i,j+1] - ϕ[i,j])/dy
            elseif j == par.Ny # Boundary condition
                dyϕ[i,j] = (ϕ[i,j] - ϕ[i,j-1])/dy
            else # Apply min-mod at interior grid points
                dyϕ[i,j] = minmod(par.θ*(ϕ[i,j+1]-ϕ[i,j])/dy, (ϕ[i,j+1]-ϕ[i,j-1])/(2*dy), par.θ*(ϕ[i,j]-ϕ[i,j-1])/dy)
            end
        end
    end
    # Compute sign function on grid points at τ_n
    for i = 1:par.Nx
        for j = 1:par.Ny
            S[i,j] = mod_sign(ϕ[i,j], dxϕ[i,j], dyϕ[i,j], dx)
        end
    end
    # Compute solutions on staggered grid
    ϕs = lt_staggered_reinit(S, ϕ, dxϕ, dyϕ, par, dx, dy, dt)
    # Compute derivatives on staggered grid
    dxϕs = Array{Float64}(undef, par.Nx-1, par.Ny-1) # Pre-allocate x-derivatives at staggered grid points
    dyϕs = Array{Float64}(undef, par.Nx-1, par.Ny-1) # Pre-allocate y-derivatives at staggered grid points
    for i = 1:par.Nx-1
        for j = 1:par.Ny-1
            if i == 1
                dxϕs[i,j] = (ϕs[i+1,j]-ϕs[i,j])/dx # ϕ'_{i+1/2,j+1/2}/Δx
            elseif i == par.Nx-1
                dxϕs[i,j] = (ϕs[i,j]-ϕs[i-1,j])/dx # ϕ'_{i+1/2,j+1/2}/Δx
            else
                dxϕs[i,j] = minmod(par.θ*(ϕs[i+1,j]-ϕs[i,j])/dx, (ϕs[i+1,j]-ϕs[i-1,j])/(2*dx), par.θ*(ϕs[i,j]-ϕs[i-1,j])/dx) # ϕ'_{i+1/2,j+1/2}/Δx
            end
            if j == 1
                dyϕs[i,j] = (ϕs[i,j+1]-ϕs[i,j])/dy # ϕ`_{i+1/2,j+1/2}/Δy
            elseif j == par.Ny-1
                dyϕs[i,j] = (ϕs[i,j]-ϕs[i,j-1])/dy # ϕ`_{i+1/2,j+1/2}/Δy
            else
                dyϕs[i,j] = minmod(par.θ*(ϕs[i,j+1]-ϕs[i,j])/dy, (ϕs[i,j+1]-ϕs[i,j-1])/(2*dy), par.θ*(ϕs[i,j]-ϕs[i,j-1])/dy) # ϕ`_{i+1/2,j+1/2}/Δy
            end
        end
    end
    # Perform cell-averaging at non-staggered points
    for i = 1:par.Nx
        for j = 1:par.Ny
            if (i == 1) || (j == 1) || (i == par.Nx) || (j == par.Ny) # Ignore boundaries
                ϕn[i,j] = ϕ[i,j]
            else
                ϕn[i,j] = (ϕs[i,j] + ϕs[i-1,j] + ϕs[i,j-1] + ϕs[i-1,j-1])/4 +
                (dx*(dxϕs[i-1,j-1]-dxϕs[i,j-1]) + dx*(dxϕs[i-1,j]-dxϕs[i,j]) + dy*(dyϕs[i-1,j-1]-dyϕs[i-1,j]) + dy*(dyϕs[i,j-1]-dyϕs[i,j]))/16 # Compute non-staggered cell average
            end
        end
    end
    return ϕn
end

"Implement staggered Lin-Tadmor scheme"
function lt_staggered_reinit(S, ϕ, dxϕ, dyϕ, par, dx, dy, dt)
    ϕ_stag = Array{Float64}(undef, par.Nx-1, par.Ny-1) # Pre-allocate solution on staggered grid ϕ_stag[i,j] = ϕ[i+1/2, j+1/2], one-based indexing
    S_stag = Array{Float64}(undef, par.Nx-1, par.Ny-1) # Pre-allocate velocity on staggered grid V_stag[i,j] = ϕ[i+1/2, j+1/2], one-based indexing
    ϕ_mid = Array{Float64}(undef, par.Nx, par.Ny) # Pre-allocate mid-values at non-staggered grid points
    # Evaluate missing mid-values on full non-staggered grid
    for i = 1:par.Nx
        for j = 1:par.Ny
            ϕ_mid[i,j] = ϕ[i,j] - dt/2*hamiltonian_reinit(S[i,j], dxϕ[i,j], dyϕ[i,j]) # Missing mid-values
        end
    end
    # Obtain sign function on staggered grid
    for i = 1:par.Nx-1
        for j = 1:par.Ny-1
            S_stag[i,j] = (S[i,j] + S[i+1,j] + S[i,j+1] + S[i+1,j+1])/4
        end
    end
    # Evaluate solution on staggered grid
    for i = 1:par.Nx-1
        for j = 1:par.Ny-1
            ϕ_stag[i,j] = (ϕ[i,j] + ϕ[i+1,j] + ϕ[i,j+1] + ϕ[i+1,j+1])/4 +
            (dxϕ[i,j]-dxϕ[i+1,j]+dxϕ[i,j+1]-dxϕ[i+1,j+1])*dx/16 +
            (dyϕ[i,j]-dyϕ[i,j+1]+dyϕ[i+1,j]-dyϕ[i+1,j+1])*dy/16 -
            dt/2*(hamiltonian_reinit(S_stag[i,j], (ϕ_mid[i+1,j]-ϕ_mid[i,j])/dx, (ϕ_mid[i+1,j+1]-ϕ_mid[i+1,j])/dy) +
            hamiltonian_reinit(S_stag[i,j], (ϕ_mid[i+1,j+1]-ϕ_mid[i,j+1])/dx, (ϕ_mid[i,j+1]-ϕ_mid[i,j])/dy))
        end
    end
    return ϕ_stag
end

"Implement the modified sign function"
function mod_sign(ϕ, ϕx, ϕy, dx)
    return ϕ/sqrt(ϕ^2 + dx^2*(ϕx^2 + ϕy^2))
end

"Hamiltonian for reinitialisation equation"
function hamiltonian_reinit(s, ϕx, ϕy)
    return s*(sqrt(ϕx^2 + ϕy^2) - 1)
end


text:
@testset "Polygon2D" begin
    verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2, 2)]

    # cases for closed = true
    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    # cases for closed = false
    img = zeros(RGB, 7, 7)
    expected = ones(RGB, 7, 7)
    res = @inferred draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = false)
    @test all(expected .== res) == true

    img = zeros(RGB, 7, 7)
    expected = ones(RGB, 7, 7)
    res = @inferred draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = false)
    @test all(expected .== res) == true

    # cases without boundary_value
    img = zeros(RGB,7,7)
    expected = ones(RGB,7,7)
    res = @inferred draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1)); closed = false)
    @test all(expected .== res) == true

    img = zeros(RGB,7,7)
    expected = ones(RGB,7,7)
    res = @inferred draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(1)); closed = false)
    @test all(expected .== res) == true

    # cases without both fill_value and boundary_value
    img = zeros(RGB,7,7)
    expected = ones(RGB,7,7)
    res = @inferred draw(img, verts, BoundaryFill(4, 4); closed = false)
    @test all(expected .== res) == true

    img = zeros(RGB,7,7)
    expected = ones(RGB,7,7)
    res = @inferred draw!(img, verts, BoundaryFill(4, 4); closed = false)
    @test all(expected .== res) == true

    # cases with no parameters
    img = zeros(RGB,7,7)
    expected = copy(img)
    expected .= RGB{N0f8}(1)
    res = @inferred draw(img, verts, BoundaryFill(); closed = false)
    @test all(expected .== res) == true

    img = zeros(RGB,7,7)
    expected = copy(img)
    expected .= RGB{N0f8}(1)
    res = @inferred draw!(img, verts, BoundaryFill(); closed = false)
    @test all(expected .== res) == true

    # cases for use of CartesianIndex
    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw(img, verts, BoundaryFill(CartesianIndex(4, 4); fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw!(img, verts, BoundaryFill(CartesianIndex(4, 4); fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    # Cases for use of Point
    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw(img, verts, BoundaryFill(Point(5, 4); fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    img = zeros(RGB, 7, 7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB{N0f8}(1)
    res = @inferred draw!(img, verts, BoundaryFill(Point(5, 4); fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
    @test all(expected .== res) == true

    # cases when verts are outside the domain
    img = zeros(RGB, 7, 7)
    verts = [CartesianIndex(2, 2), CartesianIndex(2, 0), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2, 2)]
    err = DomainError(CartesianIndex(2, 0), "Vertice $(CartesianIndex(2, 0)) outside the image array domain")
    @test_throws err draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)

    img = zeros(RGB, 7, 7)
    verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(0, 6), CartesianIndex(6, 2), CartesianIndex(2, 2)]
    err = DomainError(CartesianIndex(0, 6), "Vertice $(CartesianIndex(0, 6)) outside the image array domain")
    @test_throws err draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)

    # verts need to be reset to original
    verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2, 2)]

    # cases when initial point is outside image domain
    img = zeros(RGB, 7, 7)
    @test_throws BoundsError draw(img, verts, BoundaryFill(0, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)

    img = zeros(RGB, 7, 7)
    @test_throws BoundsError draw!(img, verts, BoundaryFill(4, 9; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)

    # cases without boundary_value and fill value = RGB(0, 1, 0)
    img = zeros(RGB,7,7)
    expected = copy(img)
    expected[2:6, 2:6] .= RGB(0, 1, 0)
    res = @inferred draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(0, 1, 0)); closed = true)
    @test all(expected .== res) == true

    img = zeros(RGB,7,7)
    expected = copy(img)
    expected[:,:] .= RGB(0, 1, 0)
    res = @inferred draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(0, 1, 0)); closed = false)
    @test all(expected .== res) == true

end


text: abstract type AbstractDataPoint end # For connecting Dataset module's AbstractDataPoint to main module


text:
@testset "ExponentialLikelihood" begin
    for args in ((), (exp,), (ExpLink(),))
        @test ExponentialLikelihood(args...) isa ExponentialLikelihood{ExpLink}
    end

    lik = ExponentialLikelihood()
    test_interface(lik, SqExponentialKernel(), rand(10))
end


text:
#######################################
## Conversion (assumes input isvalid)
#######################################

function natparam(g::GaussMP)::GaussNP
    P = pinv(2g.mu2.data-g.mu1*g.mu1')
    GaussianNatParam(P*g.mu1, -P)
end
function natparam(g::DGaussMP)::DGaussNP
    P = 1./(2g.mu2-g.mu1.^2)
    DiagGaussianNatParam(P.*g.mu1, -P)
end
function meanparam(g::GaussNP, correction::Float=1.0)::GaussMP
    # NOTE correction -> may want to use (N-P-2)/(N-1) if g.theta2 is noisy,
    # NOTE see litt for estimation of precision matrix.
    cov = -pinv(g.theta2.data)*correction
    mu1 = cov * g.theta1
    GaussianMeanParam(mu1, (mu1*mu1' + cov)/2)
end
function meanparam(g::DGaussNP)::DGaussMP
    cov = -1./g.theta2
    mu1 = cov.*g.theta1
    DiagGaussianMeanParam(mu1, (mu1.^2+cov)/2)
end

Base.full(g::DGaussNP) = GaussianNatParam(g.theta1, diagm(g.theta2))
Base.full(g::DGaussMP) = GaussianMeanParam(g.mu1, diagm(g.mu2))

##############################
## Operations extending base
##############################

Base.mean(g::GaussNP)  = -g.theta2\g.theta1
Base.mean(g::DGaussNP) = -g.theta1./g.theta2
Base.mean(g::GaussMP)  = g.mu1
Base.mean(g::DGaussMP) = g.mu1
Base.cov(g ::GaussNP)  = -pinv(g.theta2.data)
Base.cov(g ::DGaussNP) = -1./g.theta2
Base.cov(g ::GaussMP)  = 2g.mu2.data-g.mu1*g.mu1'
Base.cov(g ::DGaussMP) = 2g.mu2-g.mu1.^2
Base.var(g ::FGauss)   = diag(cov(g))
Base.var(g ::DGauss)   = cov(g)
Base.std(g ::Gauss)    = sqrt.(var(g))

Base.isvalid(g::GaussNP)  = isposdef(-g.theta2.data)
Base.isvalid(g::DGaussNP) = all(-g.theta2 .> 0)
Base.isvalid(g::GaussMP)  = isposdef(2g.mu2.data-g.mu1*g.mu1')
Base.isvalid(g::DGaussMP) = all(2g.mu2-g.mu1.^2 .> 0)

Base.ones(::Type{GaussNP}, d::Int) =
    GaussianNatParam(mean=zeros(d), cov=eye(d))
Base.ones(::Type{DGaussNP}, d::Int) =
    DiagGaussianNatParam(mean=zeros(d), cov=ones(d))
Base.ones(::Type{GaussMP}, d::Int) =
    GaussianMeanParam(mean=zeros(d), cov=eye(d))
Base.ones(::Type{DGaussMP}, d::Int) =
    DiagGaussianMeanParam(mean=zeros(d), cov=ones(d))
# NOTE zeros does not really make sense, just use a one with variance
# to infinity or to zero based on use case

Base.length(g::Union{GaussNP, DGaussNP}) = length(g.theta1)
Base.length(g::Union{GaussMP, DGaussMP}) = length(g.mu1)

Base.rand(g::GaussNP, n::Int=1) =
    broadcast(+, mean(g), chol(-g.theta2)\randn(length(g), n))
Base.rand(g::GaussMP, n::Int=1) =
    broadcast(+, mean(g), chol(cov(g))*randn(length(g), n))
Base.rand(g::Union{DGaussNP, DGaussMP}, n::Int=1) =
    broadcast(+, mean(g), broadcast(*, std(g), randn(length(g),n)))

Base.vec(g::GaussNP)  = [g.theta1;g.theta2[:]]
Base.vec(g::DGaussNP) = [g.theta1;g.theta2]
Base.vec(g::GaussMP)  = [g.mu1;g.mu2[:]]
Base.vec(g::DGaussMP) = [g.mu1;g.mu2]

###################
## Safe operators
###################

+(g1::GaussNP, g2::GaussNP) =
    GaussianNatParam(vec(g1)+vec(g2), length(g1))
+(g1::GaussMP, g2::GaussMP) =
    GaussianMeanParam(vec(g1)+vec(g2), length(g1))
+(g1::DGaussNP, g2::DGaussNP) =
    DiagGaussianNatParam(vec(g1)+vec(g2), length(g1))
+(g1::DGaussMP, g2::DGaussMP) =
    DiagGaussianMeanParam(vec(g1)+vec(g2), length(g1))

*(a::Real,  g::GaussNP)  = GaussianNatParam(a*vec(g), length(g))
*(a::Real,  g::DGaussNP) = DiagGaussianNatParam(a*vec(g), length(g))
*(a::Real,  g::GaussMP)  = GaussianMeanParam(a*vec(g), length(g))
*(a::Real,  g::DGaussMP) = DiagGaussianMeanParam(a*vec(g), length(g))

*(g::Gauss, a::Real)   = *(a,g)
/(a::Real,  g::Gauss)   = *(1.0/a,g)
/(g::Gauss, a::Real)   = /(a,g)

#####################
## Unsafe operators (NOTE result is not necessarily valid Gaussian)
#####################

-(g1::GaussNP, g2::GaussNP) =
    GaussianNatParam(vec(g1)-vec(g2), length(g1))
-(g1::GaussMP, g2::GaussMP) =
    GaussianMeanParam(vec(g1)-vec(g2), length(g1))
-(g1::DGaussNP, g2::DGaussNP) =
    DiagGaussianNatParam(vec(g1)-vec(g2), length(g1))
-(g1::DGaussMP, g2::DGaussMP) =
    DiagGaussianMeanParam(vec(g1)-vec(g2), length(g1))

#########################
## Comparison operators
#########################

Base.norm(g::Gauss) = norm(vec(g))

function Base.isapprox(g1::Gauss, g2::Gauss;
                       rtol::Real=sqrt(eps()), atol::Real=0)
    norm(g1-g2) <= atol + rtol*max(norm(g1), norm(g2))
end

########################
## Projection operator
########################

function project(g::GaussNP; minmu::Float=-Inf, maxmu::Float=Inf,
                 minvar::Float=0., maxvar::Float=Inf)::GaussNP
    (D, V)  = eig(-g.theta2)
    D, V    = real(D), real(V)
    Dthresh = max.(1.0/maxvar, min.(1.0/minvar, D))
    mthresh = max.(minmu, min.(maxmu, mean(g)))
    Pthresh = V*Diagonal(Dthresh)*V'
    GaussianNatParam(Pthresh*mthresh, -Pthresh)
end
function project(g::GaussMP; minmu::Float=-Inf, maxmu::Float=Inf,
                 minvar::Float=0., maxvar::Float=Inf)::GaussMP
    (D,V)   = eig(cov(g))
    D, V    = real(D), real(V)
    Dthresh = max.(minvar, min.(maxvar, D))
    mthresh = max.(minmu, min.(maxmu, mean(g)))
    GaussianMeanParam(mean=mthresh, cov=V*Diagonal(Dthresh)*V')
end
function project(g::DGaussNP; minmu=-Inf, maxmu=Inf,
                 minvar=0, maxvar=Inf)::DGaussNP
    mthresh = max.(minmu,  min.(maxmu,  mean(g)))
    pthresh = max.(1.0 / maxvar, min.(1.0 / minvar, -g.theta2))
    DiagGaussianNatParam(pthresh.*mthresh, -pthresh)
end
function project(g::DGaussMP; minmu=-Inf, maxmu=Inf,
                 minvar=0, maxvar=Inf)::DGaussMP
    mthresh = max.(minmu,  min.(maxmu,  mean(g)))
    vthresh = max.(minvar, min.(maxvar, var(g)))
    DiagGaussianMeanParam(mean=mthresh, cov=vthresh)
end


text:
@doc raw"""
	Spinor(α,β)

Spinor(α,β) with Cayley-Klein parameters α and β.

Based on "Introduction to the Shinnar-Le Roux algorithm", Patrick Le Roux (1995).

A spinor is a way to represent 3D rotations, the underlying representation is a 2 X 2 complex unitary matrix (``\alpha,\beta\in\mathbb{C}``):

```math
R=\left[\begin{array}{cc}
\alpha & -\beta^{*}\\
\beta & \alpha^{*}
\end{array}\right],
```
with ``|\alpha|^2+|\beta|^2 = 1``.

This later operates on the ``2\times2`` representation of ``(x,y,z)`` as follows ``V^{+} = R V R^{*}``.
"""
struct Spinor
	α::Complex
	β::Complex
end
Base.show(io::IO,s::Spinor) = print(io, "Spinor(α = ",round(s.α,digits=3),", β = ",round(s.β,digits=3),")")

"""
Spinor multiplication identity.
	(α1,β1)×(α2,β2) = (α1 α2 - β2⋆ β1 , β2 α1 + α2⋆ β1)
"""
*(s1::Spinor,s2::Spinor) = begin
	Spinor(s1.α*s2.α - conj(s2.β)*s1.β,
		   s1.α*s2.β + conj(s2.α)*s1.β)
end

#Rotation matrices
"Spinor clockwise rotation matrix with angle φ with respect to z-axis."
Rz(φ) = Spinor(exp(-im*φ/2), 0)
"Spinor clockwise rotation matrix with angle Θ with respect to y-axis."
Ry(θ) = Spinor(cos(θ/2), sin(θ/2))
"""
Spinor clockwise rotation matrix with angle Θ with respect to x-axis.

Rx(θ) = Rz(-π/2) Ry(θ) Rz(π/2)
"""
Rx(θ) = Spinor(cos(θ/2), -im*sin(θ/2))
"""
Spinor rotation matrix.

Rg(φ1,θ,φ2) = Rz(φ2) Ry(θ) Rz(φ1)
"""
Rg(φ1,θ,φ2) = Spinor(cos(θ/2)*exp(-im*(φ1+φ2)/2), sin(θ/2)*exp(-im*(φ1-φ2)/2))
"""
Spinor rotation matrix with angle θ with axis in the xy plane u=(cosφ,sinφ).

Rφ(φ,θ) = Rg(-φ,θ,φ) = Rz(φ) Ry(θ) Rz(-φ)
"""
Rφ(φ,θ) = Spinor(cos(θ/2),exp(im*φ)*sin(θ/2))

@doc raw"""
Pauly, J., Le Roux, P., Nishimura, D., & Macovski, A. (1991).
Parameter relations for the Shinnar-Le Roux selective excitation pulse design algorithm (NMR imaging).
IEEE Transactions on Medical Imaging, 10(1), 53–65. doi:10.1109/42.75611

Spinor rotation matrix.

Rotation of φ with respect to the axis of rotation n=(nx,ny,nz).

```math
\varphi=-\gamma\Delta t\sqrt{\left|B_{1}\right|^{2}+\left(\boldsymbol{G}\cdot\boldsymbol{x}\right)^{2}}=-\gamma\Delta t\left\Vert \boldsymbol{B}\right\Vert
```
```math
\boldsymbol{n}=\boldsymbol{B}/\left\Vert \boldsymbol{B}\right\Vert
```
"""
Q(φ, nxy, nz) = begin
	Spinor( cos(φ/2)-im*nz*sin(φ/2), -im*nxy*sin(φ/2) )
end
"""
It calculates |\\alpha|^2+|\\beta|^2 of the Cayley-Klein parameters
"""
abs(s::Spinor) = abs(s.α)^2 + abs(s.β)^2

"""RF Object"""
mutable struct RF
	A # Amplitud/Phase B1x + i B1y [T]
	T # Durations [s]
	Δf::Float64    # Frequency offset [Hz]
	delay::Float64 # Delay [s]
	function RF(A,T,Δf,delay)
		@argcheck all(T .>= 0) && delay >= 0 "RF timings must be positive."
		new(A, T, Δf, delay)
    end
	function RF(A,T,Δf)
		@argcheck all(T .>= 0) "RF timings must be positive."
		new(A, T, Δf, 0.)
    end
	function RF(A,T)
		@argcheck all(T .>= 0) >= 0 "RF timings must be positive."
		new(A, T, 0., 0.)
    end
end
#Properties
*(α::ComplexF64, x::RF) = RF(α*x.A,x.T,x.Δf,x.delay)
"Duration `T` [s] of RF datatype."
dur(x::RF) = sum(x.T)
"Duration `T` [s] of the RF array Array{RF,1}."
dur(x::Array{RF,1}) = sum(sum(x[i].T) for i=1:size(x,1))
"Duration `T` [s] of the RF array Array{RF,2}."
dur(x::Array{RF,2}) = maximum(sum([sum(x[i,j].T) for i=1:size(x,1),j=1:size(x,2)],dims=2))
"Generate an RF sequence with amplitudes sampled from a function."
RF_fun(f::Function,T::Real,N::Int64=300) = begin
	t = range(0,stop=T,length=N)
	A = f.(t)
	RFs = RF(A,T)
end
"Calculates the flip-angle α [deg] of an RF object. α = γ ∫ B1(τ) dτ"
get_flip_angle(x::RF) = begin
	A, NA, T, NT = x.A, length(x.A), x.T, length(x.T)
	dT = T / NA * NT
	α = round(360 * γ * abs(sum(A .* dT)), digits=3) #Pulseq
	α
end
"Calculates the RF center. This calculation includes the RF delay."
get_RF_center(x::RF) = begin
	A, NA, T, NT, delay = x.A, length(x.A), x.T, length(x.T), x.delay
	dT = T / NA * NT .* ones(NA)
	t = cumsum([0; dT])[1:end-1]
	i_center = argmax(abs.(A))
	t_center = t[i_center] + dT[i_center]/2
	t_center + delay
end

one(T::Spinor) = Spinor(1,0)
getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x,f)
getproperty(x::Matrix{RF}, f::Symbol) = begin
	if     f == :Bx
		real.(getproperty.(x,:A))
	elseif f == :By
		imag.(getproperty.(x,:A))
	elseif f == :Δf
		getproperty.(x,:Δf)
	elseif f == :T || f == :delay
		getproperty.(x[1,:],f)
	elseif f == :dur
		T, delay = x.T, x.delay
		ΔT = [sum(t) for t=T] .+ delay
		ΔT
	else
		getproperty.(x,f)
	end
end
#aux
Base.show(io::IO,x::RF) = begin
	r(x) = round.(x,digits=4)
	compact = get(io, :compact, false)
	if !compact
		wave = length(x.A) == 1 ? r(x.A*1e6) : "∿"
		print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"RF($(wave) uT, $(r(sum(x.T)*1e3)) ms, $(r(x.Δf)) Hz)")
	else
		wave = length(x.A) == 1 ? "⊓" : "∿"
		print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿")*"($(r((x.delay+sum(x.T))*1e3)) ms)")
	end
end


text:
immutable CurvePt
    pos::VecSE2 # global position and orientation
    s::Float64  # distance along the curve
    k::Float64  # curvature
    kd::Float64 # derivative of curvature

    CurvePt(pos::VecSE2, s::Float64, k::Float64=NaN, kd::Float64=NaN) = new(pos, s, k, kd)
end
Base.show(io::IO, pt::CurvePt) = @printf(io, "CurvePt({%.3f, %.3f, %.3f}, %.3f, %.3f, %.3f)", pt.pos.x, pt.pos.y, pt.pos.θ, pt.s, pt.k, pt.kd)

Vec.lerp(a::CurvePt, b::CurvePt, t::Float64) = CurvePt(lerp(a.pos, b.pos, t), a.s + (b.s - a.s)*t, a.k + (b.k - a.k)*t, a.kd + (b.kd - a.kd)*t)

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

typealias Curve Vector{CurvePt}

"""
    get_lerp_time_unclamped(A::VecE2, B::VecE2, Q::VecE2)
Get the interpolation scalar t for the point on the line AB closest to Q
This point is P = A + (B-A)*t
"""
function get_lerp_time_unclamped(A::VecE2, B::VecE2, Q::VecE2)

    a = Q - A
    b = B - A
    c = proj(a, b, VecE2)

    if b.x != 0.0
        t = c.x / b.x
    elseif b.y != 0.0
        t = c.y / b.y
    else
        t = 0.0 # no lerping to be done
    end

    t
end
get_lerp_time_unclamped(A::VecSE2, B::VecSE2, Q::VecSE2) = get_lerp_time_unclamped(convert(VecE2, A), convert(VecE2, B), convert(VecE2, Q))
get_lerp_time_unclamped(A::CurvePt, B::CurvePt, Q::VecSE2) = get_lerp_time_unclamped(convert(VecE2, A.pos), convert(VecE2, B.pos), convert(VecE2, Q))

"""
    get_lerp_time(A::VecE2, B::VecE2, Q::VecE2)
Get lerp time t∈[0,1] such that lerp(A, B) is as close as possible to Q
"""
get_lerp_time(A::VecE2, B::VecE2, Q::VecE2) = clamp(get_lerp_time_unclamped(A, B, Q), 0.0, 1.0)
get_lerp_time(A::CurvePt, B::CurvePt, Q::VecSE2) = get_lerp_time(convert(VecE2, A.pos), convert(VecE2, B.pos), convert(VecE2, Q))

immutable CurveIndex
    i::Int     # index in curve, ∈ [1:length(curve)-1]
    t::Float64 # ∈ [0,1] for linear interpolation
end
const CURVEINDEX_START = CurveIndex(1,0.0)
Base.show(io::IO, ind::CurveIndex) = @printf(io, "CurveIndex(%d, %.3f)", ind.i, ind.t)

curveindex_end(curve::Curve) = CurveIndex(length(curve)-1,1.0)

Base.getindex(curve::Curve, ind::CurveIndex) = lerp(curve[ind.i], curve[ind.i+1], ind.t)

function is_at_curve_end(ind::CurveIndex, curve::Curve)
    (ind.i == 1 && ind.t == 0.0) ||
    (ind.i == length(curve)-1 && ind.t == 1.0)
end

function index_closest_to_point(curve::Curve, target::AbstractVec)

    a = 1
    b = length(curve)
    c = div(a+b, 2)

    @assert(length(curve) ≥ b)

    sqdist_a = abs2(curve[a].pos - target)
    sqdist_b = abs2(curve[b].pos - target)
    sqdist_c = abs2(curve[c].pos - target)

    while true
        if b == a
            return a
        elseif b == a + 1
            return sqdist_b < sqdist_a ? b : a
        elseif c == a + 1 && c == b - 1
            if sqdist_a < sqdist_b && sqdist_a < sqdist_c
                return a
            elseif sqdist_b < sqdist_a && sqdist_b < sqdist_c
                return b
            else
                return c
            end
        end

        left = div(a+c, 2)
        sqdist_l = abs2(curve[left].pos - target)

        right = div(c+b, 2)
        sqdist_r = abs2(curve[right].pos - target)

        if sqdist_l < sqdist_r
            b = c
            sqdist_b = sqdist_c
            c = left
            sqdist_c = sqdist_l
        else
            a = c
            sqdist_a = sqdist_c
            c = right
            sqdist_c = sqdist_r
        end
    end

    error("index_closest_to_point reached unreachable statement")
end

"""
    get_curve_index(curve::Curve, s::Float64)
Return the CurveIndex for the closest s-location on the curve
"""
function get_curve_index(curve::Curve, s::Float64)

    if s ≤ 0.0
        return CURVEINDEX_START
    elseif s ≥ curve[end].s
        return curveindex_end(curve)
    end

    a = 1
    b = length(curve)

    fa = curve[a].s - s
    fb = curve[b].s - s

    n = 1
    while true
        if b == a+1
            extind = a + -fa/(fb-fa)
            ind = floor(Int, extind)
            t = rem(extind, 1.0)
            return CurveIndex(ind, t)
        end

        c = div(a+b, 2)
        fc = curve[c].s - s
        n += 1

        if sign(fc) == sign(fa)
            a, fa = c, fc
        else
            b, fb = c, fc
        end
    end

    error("get_curve_index failed for s=$s")
end

"""
    get_curve_index(ind::CurveIndex, curve::Curve, Δs::Float64)
Return the CurveIndex at ind's s position + Δs
"""
function get_curve_index(ind::CurveIndex, curve::Curve, Δs::Float64)

    L = length(curve)
    ind_lo, ind_hi = ind.i, ind.i+1

    s_lo = curve[ind_lo].s
    s_hi = curve[ind_hi].s
    s = lerp(s_lo, s_hi, ind.t)

    if Δs ≥ 0.0

        if s + Δs ≥ s_hi && ind_hi < L
            while s + Δs ≥ s_hi && ind_hi < L
                Δs -= (s_hi - s)
                s = s_hi
                ind_lo += 1
                ind_hi += 1
                s_lo = curve[ind_lo].s
                s_hi = curve[ind_hi].s
            end
        else
            Δs = s + Δs - s_lo
        end

        t = Δs/(s_hi - s_lo)
        CurveIndex(ind_lo, t)
    else
        while s + Δs < s_lo  && ind_lo > 1
            Δs += (s - s_lo)
            s = s_lo
            ind_lo -= 1
            ind_hi -= 1
            s_lo = curve[ind_lo].s
            s_hi = curve[ind_hi].s
        end

        Δs = s + Δs - s_lo
        t = Δs/(s_hi - s_lo)
        CurveIndex(ind_lo, t)
    end
end

"""
    CurveProjection
The result of a point projected to a Curve
"""
immutable CurveProjection
    ind::CurveIndex
    t::Float64 # lane offset
    ϕ::Float64 # lane-relative heading [rad]
end
Base.show(io::IO, curveproj::CurveProjection) = @printf(io, "CurveProjection({%d, %.3f}, %.3f, %.3f)", curveproj.ind.i, curveproj.ind.t, curveproj.t, curveproj.ϕ)
function get_curve_projection(posG::VecSE2, footpoint::VecSE2, ind::CurveIndex)
    F = inertial2body(posG, footpoint)
    CurveProjection(ind, F.y, F.θ)
end

"""
    Vec.proj(posG::VecSE2, curve::Curve)
Return a CurveProjection obtained by projecting posG onto the curve
"""
function Vec.proj(posG::VecSE2, curve::Curve)

    ind = index_closest_to_point(curve, posG)::Int

    curveind = CurveIndex(0,NaN)
    footpoint = VecSE2(NaN, NaN, NaN)

    if ind > 1 && ind < length(curve)
        t_lo = get_lerp_time(curve[ind-1], curve[ind],   posG)
        t_hi = get_lerp_time(curve[ind],   curve[ind+1], posG)

        p_lo = lerp(curve[ind-1].pos, curve[ind].pos,   t_lo)
        p_hi = lerp(curve[ind].pos,   curve[ind+1].pos, t_hi)

        d_lo = hypot(p_lo - posG)
        d_hi = hypot(p_hi - posG)

        if d_lo < d_hi
            footpoint = p_lo
            curveind = CurveIndex(ind-1, t_lo)
        else
            footpoint = p_hi
            curveind = CurveIndex(ind, t_hi)
        end
    elseif ind == 1
        t = get_lerp_time( curve[1], curve[2], posG )
        footpoint = lerp( curve[1].pos, curve[2].pos, t)
        curveind = CurveIndex(ind, t)
    else # ind == length(curve)
        t = get_lerp_time( curve[end-1], curve[end], posG )
        footpoint = lerp( curve[end-1].pos, curve[end].pos, t)
        curveind = CurveIndex(ind-1, t)
    end

    get_curve_projection(posG, footpoint, curveind)
end


text:
function nyquistPlot(Z; fignum=1, label="")
    s = subplot(111)
    title("Nyquist plot")
    xlabel("Re\$(Z) \\ [\\Omega]\$")
    ylabel("-Im\$(Z) \\ [\\Omega]\$")
    PyPlot.plot(real(Z), -imag(Z), label=label, "-x")
    if label!=""
      legend()
    end
    grid(true)
    s.set_aspect(1.0)
end

function get_RC_prms_from_Z(omega, Z)
  a = real(Z)
  b = imag(Z)
  return (
    (a*a + b*b)/a,
    - b/(omega*(a*a + b*b))
    )
end

function get_R_RC_prms_from_2Z(omegas, Zs)
  a1, b1 = real(Zs[1]), imag(Zs[1])
  a2, b2 = real(Zs[2]), imag(Zs[2])

  tol = 1.0e-4
  if abs(a1 - a2) < tol*(a1+a2)/2
    R_ohm = a1
    R = 0
    C = 0
  else

    A = a1-a2
    B = ((a1^2 + b1^2) - (a2^2 + b2^2))
    C = a2*(a1^2 + b1^2) - a1*(a2^2 + b2^2)

    D = B^2 - 4*A*C

    s1 = (-B + sqrt(D))/
            (2*A)

    R_ohm = -s1
    R, C = get_RC_prms_from_Z(omegas[1], Zs[1] - R_ohm)
    #get_RC_prms_from_Z(omegas[2], Zs[2] - R_ohm)
  end

  return R_ohm, R, C
end

function get_Z_from_R_RC(f_range, R_ohm, R, C)
  [R_ohm + R/(1 + im*2*pi*f*R*C) for f in f_range]
end

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
#####
## Functions for preprocessing data in TintiNet.jl
#

"""
    pad_aa_seq(v::Vector{String}, n::Integer, p::String)

    This function receives a `Vector{String}` `v` and returns a same-length `Vector{String}`
    with every element of `v` padded to length `n` with symbol `p`
"""
pad_aa_seq(v::Vector{String}, n::Integer, p::String) = [v; fill(p, max(n - length(v), 0))]

"""
    preprocess(batchofseqs::Vector{Vector{String}}, n::Integer, p::String)

    This function applies `pad_aa_seq` to every element of `batchofseqs` with arguments
    `n` and `p` (default 120 and \"-\")
"""
preprocess(batchofseqs::Vector{Vector{String}}, n = 128, p = "-") = [ pad_aa_seq([singleseq...], n, "-") for singleseq in batchofseqs ]

"""
    prepare_mask(batchofseqs::Vector{Vector{String}})

    This function returns a mask to be used on the masked loss function for training.
    The mask is an `Array{Float}` filled with ones in positions different from `unk` and zeros otherwise.
    This function makes use of `Transformers.Basic.getmask`
"""
function prepare_mask(batchofseqs::Vector{Vector{String}}, seqlen)

    mask = getmask(preprocess(batchofseqs, seqlen)) # Uses getmask() from Transformers.Basic

    for i in 1:length(batchofseqs)
        for j in 1:length(batchofseqs[i])
            if batchofseqs[i][j] == "-"
                mask[1, j, i] = 0.0
            end
        end
    end

    return(mask)

end

"""
    split_train_test(
        sequences::Vector{Vector{String}},
        structures::Vector{Vector{String}},
        maximum_length::Int,
        minimum_length=32;
        compute_complete_only=true,
        test_all_incomplete_obs=false,
        training_percentage=0.6)

    This function splits a collection of target inputs (`sequences`) and outputs (`structures`) into a test set and training set.
    Sequences/structures are filtered off if their length is smaller than `minimum_length`.
    Sequences/structures are eith length greater than `maximum_length` are trimmed.

    If `compute_complete_only` is `false` and `test_all_incomplete_obs` is `false`,
    all samples will be split according to the fraction `training_percentage`, independant of their completeness, by default.

    If `compute_complete_only` is `true` and `test_all_incomplete_obs` is `false`,
    samples with `unk` positions will be removed and the remaining samples will be split according to the fraction `training_percentage`.

    If `compute_complete_only` is `false` and `test_all_incomplete_obs` is `true`,
    samples with `unk` positions will be placed entirely on the test set, with remaining samples placed on train set, _regardless of_ `training_percentage`.

    If `compute_complete_only` is `true` and `test_all_incomplete_obs` is `true`,
    all samples will be placed on the training set. This is not an adequate use of split_train_test for practical purposes, but can be useful for debugging.
"""
function split_train_test(
            sequences::Vector{Vector{String}},
            structures::Vector{Vector{String}},
            maximum_length::Int,
            minimum_length=32;
            compute_complete_only=false,
            test_on_incomplete_obs=false,
            training_percentage=0.7)

    ## Step 1. Size filters

    sequences = [ el[1:minimum([length(el), maximum_length])] for el in sequences ]
    structures = [ el[1:minimum([length(el), maximum_length])] for el in structures ]

    equal_length_passed_idx = [ size(el) for el in sequences ] .== [ size(el) for el in structures ]
    max_length_passed_idx = [ size(el)[1] for el in sequences ] .<= maximum_length
    min_length_passed_idx = [ size(el)[1] for el in sequences ] .>= minimum_length

    filtered_selection = equal_length_passed_idx .& max_length_passed_idx .& min_length_passed_idx

    filtered_sequences = sequences[filtered_selection]
    filtered_structures = structures[filtered_selection]

    ## Step 2. Splitting between training and test sets.

    if (compute_complete_only & test_all_incomplete_obs) # This is just to handle misunderstading of function arguments

        @warn "Function called with compute_complete_only = true and train_all_complete_obs = true.\nTest set will contain zero observations."

        train_sequences = filtered_sequences
        test_sequences = Vector{Vector{String}}(undef, 0)
        train_structures = filtered_structures
        test_structures = Vector{Vector{String}}(undef, 0)

    elseif (compute_complete_only) # The complete sequences/structures selected in Step 1 will be split between training and test sets.

        n = length(filtered_sequences)
        idx = shuffle(1:n)
        train_idx = view(idx, 1:floor(Int, training_percentage*n))
        test_idx = view(idx, (floor(Int, training_percentage*n)+1):n)
        lengths = length.(filtered_sequences)

        train_sequences = filtered_sequences[train_idx]
        test_sequences = filtered_sequences[test_idx]
        train_structures = filtered_structures[train_idx]
        test_structures = filtered_structures[test_idx]

    else # The complete sequences/structures will be the training set, and the incomplete sequences/structures will be the testing set.

        train_idx = [ sum(el .== "-") for el in filtered_structures ] .== 0
        test_idx = .!train_idx

        train_sequences = filtered_sequences[train_idx]
        test_sequences = filtered_sequences[test_idx]
        train_structures = filtered_structures[train_idx]
        test_structures = filtered_structures[test_idx]

    end

    ## Step 3. Ordering the sets by length

    train_order = sortperm(length.(train_sequences))
    test_order = sortperm(length.(test_sequences))

    train_sequences_ordered = train_sequences[train_order]
    test_sequences_ordered = test_sequences[test_order]
    train_structures_ordered = train_structures[train_order]
    test_structures_ordered = test_structures[test_order]
    train_lengths_ordered = length.(train_sequences_ordered)
    test_lengths_ordered = length.(test_sequences_ordered)

    return(train_sequences_ordered, test_sequences_ordered, train_structures_ordered, test_structures_ordered, train_lengths_ordered, test_lengths_ordered)

end


text:
using MultiScaleArrays, OrdinaryDiffEq, DiffEqBase, StochasticDiffEq, Statistics
using Test

struct PlantSettings{T} x::T end
struct OrganParams{T} y::T end

struct Organ{B<:Number,P} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
    name::Symbol
    params::P
end

struct Plant{B,S,N<:Tuple{Vararg{<:Organ{<:Number}}}} <: AbstractMultiScaleArray{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
    settings::S
end

struct Community{B,N<:Tuple{Vararg{<:Plant{<:Number}}}} <: AbstractMultiScaleArray{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
end

mutable struct Scenario{B,N<:Tuple{Vararg{<:Community{<:Number}}}} <: AbstractMultiScaleArrayHead{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
end

organ1 = Organ([1.1,2.1,3.1], :Shoot, OrganParams(:grows_up))
organ2 = Organ([4.1,5.1,6.1], :Root, OrganParams("grows down"))
organ3 = Organ([1.2,2.2,3.2], :Shoot, OrganParams(true))
organ4 = Organ([4.2,5.2,6.2], :Root, OrganParams(1//3))
plant1 = construct(Plant, (deepcopy(organ1), deepcopy(organ2)), Float64[], PlantSettings(1))
plant2 = construct(Plant, (deepcopy(organ3), deepcopy(organ4)), Float64[], PlantSettings(1.0))
community = construct(Community, (deepcopy(plant1), deepcopy(plant2), ))
scenario = construct(Scenario, (deepcopy(community),))

@inferred getindex(organ1, 1)
@inferred getindex(plant1, 3)
@inferred getindex(community, 4)
@inferred getindex(scenario, 8)

@test scenario[1] == 1.1
@test scenario[2] == 2.1
@test scenario[3] == 3.1
@test scenario[4] == 4.1
@test scenario[5] == 5.1
@test scenario[6] == 6.1
@test scenario[7] == 1.2
@test scenario[8] == 2.2
@test scenario[9] == 3.2
@test scenario[10] == 4.2
@test scenario[11] == 5.2
@test scenario[12] == 6.2

@test getindices(scenario, 1) == 1:12
@test getindices(scenario, 1, 1) == 1:6
@test getindices(scenario, 1, 2) == 7:12
@test getindices(scenario, 1, 1, 1) == 1:3
@test getindices(scenario, 1, 1, 2) == 4:6
@test getindices(scenario, 1, 2, 1) == 7:9
@test getindices(scenario, 1, 2, 2) == 10:12

organ_ode = function (dorgan,organ,p,t)
    m = mean(organ)
    for i in eachindex(organ)
        dorgan[i] = -m*organ[i]
    end
end
f = function (dscenario,scenario,p,t)
    for (organ, y, z) in LevelIterIdx(scenario, 2)
        organ_ode(@view(dscenario[y:z]),organ,p,t)
    end
end
affect! = function (integrator)
    add_node!(integrator, integrator.u[1, 1, 1], 1, 1)
end

println("ODE with tuple nodes")

prob = ODEProblem(f, scenario, (0.0, 1.0))

sol = solve(prob, Tsit5())

@test length(sol[end]) == 12

organ1 .+ organ2
scenario .+ scenario


text:
# functions related to hyper-geometric distribution

import .RFunctions:
    hyperpdf,
    hyperlogpdf,
    hypercdf,
    hyperccdf,
    hyperlogcdf,
    hyperlogccdf,
    hyperinvcdf,
    hyperinvccdf,
    hyperinvlogcdf,
    hyperinvlogccdf


text:
function nmex(a::Matrix{Float64}, x::Vector{Float64};
    tol::Float64=1.0e-10, limit::Int64=1000, show=0)

  iters = 0
  n = size(a, 1)
  converged = false
  local big
  while true
    iters += 1
    x1 = a * x
    big = 0.0
    for i in 1:n
      abs(x1[i]) > abs(big) && (big = x1[i])
    end
    x1 = x1 / big
    (checkit(x1, x, tol) || (iters == limit)) && break
    x = x1
    iters <= show && println("$iters $big $x")
  end
  l2 = norm(x1)
  x1 = x1 / l2
  (iters, big, x1)
end


text:
binDig = 88;
set_bigfloat_precision(binDig);
# big or BigFloat can be used on 2.0 and integers because they are exactly representable, and big(pi) should also (about) have all digits.
relacc = BigFloat(2.0)^(-binDig*big(4)/5);

for wei = 1:4

if wei == 1
	function Q(z) big(7)/10*z^3+big(3)/2; end
	alpha = big(28)/10;
	maxN = 128;
	zt = Complex{BigFloat}(1,-1)/big(100);
        function betan(n) big(n)^(1/3)*(big(3)/2*big(+7)/10*gamma(big(3.5))/gamma(big(4))/sqrt(big(pi)))^(-1/big(3)); end;
	label = "Mon0p7ThirdDegImag"
	begIn = BigFloat(0.0);

elseif wei == 2
	function Q(z) BigFloat(1.0)*z^4 -BigFloat(3)*z^2; end
	alpha = BigFloat(0.0);
	zt = big(31)/100;
	maxN = 128*2+1
        function betan(n) big(n)^(1/2)*(big(2)/2*big(+1.0)*gamma(big(2.5))/gamma(big(3))/sqrt(big(pi)))^(-1/big(2)); end
	label = "Hermite_x^4"
	begIn = BigFloat(-Inf);

elseif wei == 3
	function Q(z) BigFloat(exp(z)); end
	alpha = BigFloat(-1/2);
	zt = BigFloat(6)/10;
	maxN = 128;
        function betan(n) big(log(4*big(pi)*n) - log(log(4*big(pi)*n))); end
	label = "Exp"
	begIn = BigFloat(0);

elseif wei == 4
	function Q(z) big(1)*z^6 -big(21)/10*z^5 + big(3)*z^3 -big(6)*z^2 +big(9); end
	alpha = big(11)/10;
	maxN = 128;
	zt = Complex{BigFloat}(-1,-2);
        function betan(n) big(n)^(1/6)*(big(6)/2*big(1)*gamma(big(6.5))/gamma(big(7))/sqrt(big(pi)))^(-1/big(6)); end;
	label = "genPoly"
	begIn = BigFloat(0);
end

printtoc = 30;
anPs = zeros(BigFloat,maxN+2,1);
bnm1Ps = zeros(BigFloat,maxN+2,1);
bnm1Ps[1] = sqrt(quadgk(x -> x^alpha*exp(-Q(x)), begIn, BigFloat(+Inf), reltol=relacc)[1]);

function pn(n,z)
        pnm2 = 0;
        if n < 0
            return 0*z;
	end
        pnm1 = 1/bnm1Ps[1];
        for j = 1:n
            pnew = ((z-anPs[j])*pnm1 -bnm1Ps[j]*pnm2 )/bnm1Ps[j+1];
            pnm2 = pnm1;
            pnm1 = pnew;
	end
        return pnm1;
end

function saveAll()
# open(string("/home/path/bnm1Ps", label, ".txt"), "w") do g # to specify another path
open(string("bnm1Ps", label, ".txt"), "w") do g
        write(g, string(bnm1Ps));
end
open(string("anPs", label, ".txt"), "w") do g
        write(g, string(anPs));
end
# writedlm(string("/bnm1Ps", label, ".txt"), bnm1Ps); # Also possible

myVal = zeros(Complex{Float64}, convert(Int64,floor(log(maxN)/log(2))),1);
for ni = 1:length(myVal)
        myVal[ni] = pn(2^ni, zt*betan(2^ni));
end

open(string("result", label, ".txt"), "w") do g
        write(g, string(label, ": evaluated polys are: ", myVal) );
end
if wei == 3
        for ni = 1:length(myVal)
                myVal[ni] = bnm1Ps[2^ni+1];
        end
        open(string("coefs", label, ".txt"), "w") do g
                write(g, string(label, ": Recurrence coefficients are: ", myVal) );
        end
end

end

start = time();
ticc = start;

for nt = 1:(maxN+1)
	ea = big(0);
        if wei == 2 # We know the integrand is odd so the integral is zero, to speed up the computations.
		anPs[nt] = big(0);
	else
		(anPs[nt], ea) = quadgk(x-> x^alpha*exp(-Q(x))*x*(pn(nt-1,x))^2, begIn, BigFloat(+Inf), reltol=relacc, abstol = 2^(-binDig/6*5)); # Add abstol for when values are zero (ex. Hermite with odd integrand)
	end
        (tb, eb) = quadgk(x -> x^alpha*exp(-Q(x))*x*pn(nt-1,x)*((x-anPs[nt])*pn(nt-1,x) -bnm1Ps[nt]*pn(nt-2,x)), begIn, BigFloat(+Inf), reltol=relacc, abstol = 2^(-binDig/6*5));
	bnm1Ps[nt+1] = sqrt(tb);
	if (abs(ea) > max(abs(anPs[nt]*2*relacc), 2^(-binDig/6*5) ) ) || (abs(eb) > max(abs(tb*2*relacc), 2^(-binDig/6*5) ) )
		warn(ea, ", ", anPs[nt], " Integrals not computed accurately: ", eb, " ", tb);
	end
        if (time() > ticc + printtoc) || (nt == 10)
            ticc = time();
#            print(wei, " =wei,nt= ", nt, " ", now(), " ", Float64(anPs[nt]), " ", Float64(bnm1Ps[nt+1]), " ", Libc.strftime((maxN^2-nt^2)/nt^2*(ticc-start)+ticc), "\n"); # Print time info that is easier to read if Libc and now() are available.
            saveAll();
            print(wei, " =wei,nt= ", nt, " ", ticc, " ", convert(Float64,anPs[nt]), " ", convert(Float64,bnm1Ps[nt+1]), " ", (maxN^2-nt^2)*(ticc-start)/nt^2, "\n");
	end
end

saveAll();
print("Ended weight ", wei, ".\n")

end # loop over nonstd weights

# std Hermite:
myVal= zeros(BigFloat, 9,1)
for ni = 1:9
        n = 2^ni;
	x = sqrt(388/big(100)*n);
        pnm2 = big(0);
        pnm1 = big(pi)^(-1/4);
        for j = 1:(2*n)
            pnew = 2*x*pnm1/sqrt(big(2)*j) -(j-1)*pnm2/sqrt(j*max(1,big(j-1)));
            pnm2 = pnm1;
            pnm1 = pnew;
	end
        myVal[ni] = pnm1;
end
open(string("resultStdNormHermite.txt"), "w") do g
        write(g, string("Std normalized Hermite: evaluated polys are: ", myVal) );

end


text:
module MeshPlotter

using PyPlot

include("plot.jl")

end


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# Define a special stride where x and z are updated.
function get_stride(model::Hopper2D, traj::ContactTraj)
    stride = zeros(SizedVector{model.dim.q})
    stride[1:2] = traj.q[end-1][1:2] - traj.q[1][1:2]
    return stride
end

# get hopper model
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
s_sim = get_simulation("hopper_2D", "stairs3_2D_lc", "stairs")

# MPC parameters
N_sample = 10
H_mpc = 10
κ_mpc = 2.0e-4
n_opts = NewtonOptions(
    r_tol = 3e-4,
    max_iter = 5)
mpc_opts = LinearizedMPCOptions(
    altitude_update = true,
    altitude_impact_threshold = 0.1,
    altitude_verbose = true)

################################################################################
# Stair climbing
################################################################################

# get stair trajectory
ref_traj_ = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_stair_ref.jld2"),
    load_type=:split_traj_alt)
ref_traj = deepcopy(ref_traj_)

visualize_robot!(vis, model, ref_traj)

# Horizon
H = ref_traj.H
h = ref_traj.h
h_sim = h / N_sample
H_sim = 240*N_sample

# Objective
obj = TrackingVelocityObjective(s.model, s.env, H_mpc,
    v = [Diagonal(1e-3 * [1e-2,1,1,10]) for t = 1:H_mpc],
    q = [[Diagonal(1e-0 * [1e1,1e-1,1,1])   for t = 1:H_mpc-5]; [Diagonal(1e-1 * [1,1e-1,1e1,0.1])   for t = 1:5]],
    u = [Diagonal(1e-0 * [1e0, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1e-100 * ones(s.model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1e-100 * ones(s.model.dim.c * friction_dim(s.env))) for t = 1:H_mpc])

# Policy
p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = n_opts,
    mpc_opts = mpc_opts)

# Get initial configurations
q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

# p = open_loop_policy(ref_traj.u, N_sample = N_sample)
# Simulate the stair climbing
sim_stair = ContactImplicitMPC.simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = ContactImplicitMPC.InteriorPointOptions(
        γ_reg = 0.0,
        undercut = Inf,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

@time status = simulate!(sim_stair)

# plot_surface!(vis, s_sim.env, n=400)
anim = visualize_robot!(vis, model, sim_stair.traj, sample=10, name=:Sim, α=1.0)

plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.(sim_stair.traj.q)...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[1:model.dim.u] for u in sim_stair.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
plot!(plt[3,1], hcat(Vector.([γ[1:model.dim.c] for γ in sim_stair.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

# ghost
ref_traj_full_ = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_stairs_3_v3.jld2"),
    load_type=:split_traj_alt)
ref_traj_full = deepcopy(ref_traj_full_)

idx = [10, 300, 500, 1100, 1300, 1900, 2100]
α = [0.2, 0.2, 0.4, 0.4, 0.6, 0.6, 0.6]
hopper_parkour_ghost!(vis, sim_stair, sim_stair.traj, ref_traj_full, idx = idx, α = α)

################################################################################
# Front flip
################################################################################

# get trajectory
ref_traj_ = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_tall_flip_ref.jld2"),
    load_type=:split_traj_alt)
ref_traj = deepcopy(ref_traj_)
# offset the trajectory
for t = 1:ref_traj.H+2
    ref_traj.q[t][1] += sim_stair.traj.q[end][1]
end

# time
H = ref_traj.H
h = ref_traj.h
h_sim = h / N_sample
H_sim = 64*N_sample

# Objective
obj = TrackingVelocityObjective(s.model, s.env, H_mpc,
    v = [Diagonal(1e-13 * [1e-2,1,1,10]) for t = 1:H_mpc],
    q = [[Diagonal(1e-10 * [1e1,1e1,1,1])   for t = 1:H_mpc-5]; [Diagonal(1e-11 * [1,1e1,1e1,0.1])   for t = 1:5]],
    u = [Diagonal(1e-0 * [1e0, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1e-100 * ones(s.model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1e-100 * ones(s.model.dim.c * friction_dim(s.env))) for t = 1:H_mpc])

# Policy
p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = n_opts,
    mpc_opts = mpc_opts,
    )

# Initial configurations
q0_sim = deepcopy(SVector{model.dim.q}(sim_stair.traj.q[end-1]))
q1_sim = deepcopy(SVector{model.dim.q}(sim_stair.traj.q[end]))

sim_flip = ContactImplicitMPC.simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = ContactImplicitMPC.InteriorPointOptions(
        γ_reg = 0.0,
        undercut = Inf,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

@time status = ContactImplicitMPC.simulate!(sim_flip)

# plot_surface!(vis, s.env, xlims = [-1, 3], ylims = [-1, 1])
anim = visualize_robot!(vis, model, sim_flip.traj, sample=10, name=:Sim, α=1.0)
anim = visualize_robot!(vis, model, ref_traj, anim=anim, name=:Ref, α=0.3)

################################################################################
# Full trajectory
################################################################################

ref_traj_full = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_stairs_flip_ref.jld2"),
    load_type=:split_traj_alt)

N_sample = 10
sim_traj_full = [sim_stair.traj.q[1:end-2]; sim_flip.traj.q]
anim = visualize_robot!(vis, model,
    [sim_traj_full[1:N_sample:end]..., [sim_traj_full[end] for i = 1:50]...],
    name = :Sim, α = 1.0, h = h_sim*N_sample)
anim = visualize_robot!(vis, model,
    [ref_traj_full.q..., [ref_traj_full.q[end] for i = 1:50]...],
    anim = anim, name = :Ref, α = 0.15, h = h_sim*N_sample)
stairs!(vis)
# settransform!(vis["/Cameras/default"],
#         compose(Translation(0.0, -95.0, -1.0), LinearMap(RotY(0.0 * π) * RotZ(-π / 2.0))))
# setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 20)

plot_lines!(vis, model, ref_traj_full.q, offset = -0.3)
plot_lines!(vis, model, sim_traj_full, offset = -0.5, size = 6)

filename = "hopper_parkour_ref"
MeshCat.convert_frames_to_video(
    "/home/simon/Downloads/$filename.tar",
    "/home/simon/Documents/video/$filename.mp4", overwrite=true)

convert_video_to_gif(
    "/home/simon/Documents/video/$filename.mp4",
    "/home/simon/Documents/video/$filename.gif", overwrite=true)


text:
export partition_sparse
using SparseArrays

function partition_sparse(colptr, nnz, sz, nparts)
    nnz_per_chunk = nnz/nparts
    dom = ArrayDomain(map(x->1:x, sz))

    nxt = nnz_per_chunk
    colstart = 1
    splits = UnitRange[]

    for i=1:nparts-1
       idxs = searchsorted(colptr, nxt)
       if length(idxs) == 0
           # the exact index was not found. latch on to the nearest column boundary
           idx = first(idxs)
           if abs(colptr[idx]-nxt) < abs(colptr[idx-1]-nxt)
               nextidx = idx
           else
               nextidx = idx-1
           end
       else
           nextidx = last(idxs)
       end
       push!(splits, colstart:nextidx)
       colstart=nextidx+1
       nxt=nxt + nnz_per_chunk
    end
    push!(splits, colstart:(length(colptr)-1))

    cumlength = cumsum(map(length, splits))
    subdomains = DomainBlocks((1,1), ([size(dom, 1)], cumlength))
    DomainSplit(dom, subdomains)
end

partition_sparse(S::SparseMatrixCSC, nparts) =
    partition_sparse(S.colptr, length(S.nzval), size(S), nparts)


text:
using Plots
using LaTeXStrings

function interpolate(x1, y1, x2, y2, x)
    return ( y1 .+ (y2 .- y1)/(x2 .- x1).*(x .- x1))
end

function FindCoordinates(Ξ::Array{Skeleton,1}, t)
    tt = [Ξ[i].t for i in 1:length(Ξ)]
    i = searchsortedfirst(tt, t)
    #i = findfirst(x -> x>t ,tt)
    return(Skeleton(interpolate(Ξ[i-1].t, Ξ[i-1].ξ, Ξ[i].t, Ξ[i].ξ, t), t))
end

#x1 = Skeleton(ones(2<<1- 1), 0.0)
#x2 = Skeleton(fill(0, 2<<1 - 1), 1.0)
#FindCoordinates([x1,x2], 0.5).ξ

function plotmixing(y::Array{Skeleton,1}, b, T::Float64, L::Int64, u::Float64, v::Float64, hist, trasform = x -> x)
    # if hist
    #     ind = fill(0,L)
    #     [ind[i] =Faber(i,i) for i in 1:L]
    #     #dt = range(0.0, stop=T, length=2<<(L) + 1)
    #     P = []
    #     hist = zeros(length(b), length(ind))
    #     for i in b
    #         ξ_interp = FindCoordinates(y, i).ξ
    #         hist[j,:] = ξ_interp[ind]
    #         j = j + 1
    #     end
    #     return hist
    # else
    p = Plots.plot()
    P = []
    dt = range(0.0, stop=T, length=2<<(L) + 1)
    N = length(b)
    j = 1
    for i in b
        ξ_interp = FindCoordinates(y, i).ξ
        dx = fs_expansion(ξ_interp, u, v, L, T)
         push!(P, trasform.(dx))
    end
    Plots.plot!(p, dt, P, color = :darkrainbow, line_z = (1:N)' , linewidth=0.1, alpha = 0.2, leg = false, colorbar = true)
    return p
end

# function plotmixing(y::Array{Skeleton,1}, b, T::Float64, L::Int64, u::Float64, v::Float64, hist, trasform = x -> x)
#     if hist
#         ind = fill(0,L)
#         [ind[i] =Faber(i,i) for i in 1:L]
#         #dt = range(0.0, stop=T, length=2<<(L) + 1)
#         P = []
#         hist = zeros(length(b), length(ind))
#         for i in b
#             ξ_interp = FindCoordinates(y, i).ξ
#             hist[j,:] = ξ_interp[ind]
#             j = j + 1
#         end
#         return hist
#     else
#         p = Plots.plot()
#         dt = range(0.0, stop=T, length=2<<(L) + 1)
#         j = 1
#         for i in b
#             ξ_interp = FindCoordinates(y, i).ξ
#             dx = fs_expansion(ξ_interp, u, v, L, T)
#              push!(P, trasform.(dx))
#         end
#         Plots.plot(dt, P,  line_z = j, linewidth=0.3, alpha = 0.2, leg = false, colorbar = true)
#         return p
#     end
# end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge